source/wham/src/wham_calc1.F

   1       subroutine WHAM_CALC(islice,*)
   2 ! Weighed Histogram Analysis Method (WHAM) code
   3 ! Written by A. Liwo based on the work of Kumar et al.,
   4 ! J.Comput.Chem., 13, 1011 (1992)
   5 !
   6 ! 2/1/05 Multiple temperatures allowed.
   7 ! 2/2/05 Free energies calculated directly from data points
   8 !  acc. to Eq. (21) of Kumar et al.; final histograms also
   9 !  constructed based on this equation.
  10 ! 2/12/05 Multiple parameter sets included
  11 !
  12 ! 2/2/05 Parallel version
  13       implicit none
  14       include "DIMENSIONS"
  15       include "DIMENSIONS.ZSCOPT"
  16       include "DIMENSIONS.FREE"
  17       integer MaxBinRms,MaxBinRgy
  18       parameter (MaxBinRms=100,MaxBinRgy=100)
  19       integer MaxHdim
  20 c      parameter (MaxHdim=200000)
  21       parameter (MaxHdim=200)
  22       integer maxinde
  23       parameter (maxinde=200)
  24 #ifdef MPI
  25       include "mpif.h"
  26       include "COMMON.MPI"
  27       integer ierror,errcode,status(MPI_STATUS_SIZE)
  28 #endif
  29       include "COMMON.CONTROL"
  30       include "COMMON.IOUNITS"
  31       include "COMMON.FREE"
  32       include "COMMON.ENERGIES"
  33       include "COMMON.FFIELD"
  34       include "COMMON.SBRIDGE"
  35       include "COMMON.PROT"
  36       include "COMMON.ENEPS"
  37       integer MaxPoint,MaxPointProc
  38       parameter (MaxPoint=MaxStr,
  39      & MaxPointProc=MaxStr_Proc)
  40       double precision finorm_max,potfac,entmin,entmax,expfac,vf
  41       double precision entfac_min,entfac_min_t
  42       parameter (finorm_max=1.0d0)
  43       integer islice
  44       integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln
  45       integer start,end,iharm,ib,iib,nbin1,nbin,nbin_rms,nbin_rgy,
  46      &  nbin_rmsrgy,liczba,iparm,nFi,indrgy,indrms
  47       integer htot(0:MaxHdim),histent(0:2000)
  48       double precision v(MaxPointProc,MaxR,MaxT_h,Max_Parm)
  49       double precision energia(0:max_ene)
  50 #ifdef MPI
  51       integer tmax_t,upindE_p
  52       double precision fi_p(MaxR,MaxT_h,Max_Parm),
  53      &  fi_p_min(MaxR,MaxT_h,Max_Parm)
  54       double precision sumW_p(0:Max_GridT,Max_Parm),
  55      & sumE_p(0:Max_GridT,Max_Parm),sumEsq_p(0:Max_GridT,Max_Parm),
  56      & sumQ_p(MaxQ1,0:Max_GridT,Max_Parm),
  57      & sumQsq_p(MaxQ1,0:Max_GridT,Max_Parm),
  58      & sumEQ_p(MaxQ1,0:Max_GridT,Max_Parm),
  59      & sumEprim_p(MaxQ1,0:Max_GridT,Max_Parm),
  60      & sumEbis_p(0:Max_GridT,Max_Parm)
  61       double precision hfin_p(0:MaxHdim,maxT_h),
  62      & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
  63      & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
  64       double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
  65       double precision potEmin_t_all(maxT_h,Max_Parm),entmin_p,entmax_p
  66       integer histent_p(0:2000)
  67       logical lprint /.true./
  68 #endif
  69       double precision rgymin,rmsmin,rgymax,rmsmax
  70       double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm),
  71      & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm),
  72      & sumQsq(MaxQ1,0:NGridT,Max_Parm),sumEQ(MaxQ1,0:NGridT,Max_Parm),
  73      & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
  74      & weight,econstr
  75       double precision fi(MaxR,maxT_h,Max_Parm),
  76      & fi_min(MaxR,maxT_h,Max_Parm),
  77      & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
  78      & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
  79      & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
  80      & potEmin_all(maxT_h,Max_Parm),potEmin,potEmin_min,ent,
  81      & hfin_ent(0:MaxHdim),vmax,aux
  82       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
  83      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,
  84      &  eplus,eminus,logfac,tanhT,tt
  85       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
  86      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
  87      &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
  88      &  ehomology_constr,edfadis,edfator,edfanei,edfabet
  89
  90
  91       integer ind_point(maxpoint),upindE,indE
  92       character*16 plik
  93       character*1 licz1
  94       character*2 licz2
  95       character*3 licz3
  96       character*128 nazwa
  97       integer ilen
  98       external ilen
  99       write(licz2,'(bz,i2.2)') islice
 100       nbin1 = 1.0d0/delta
 101       write (iout,'(//80(1h-)/"Solving WHAM equations for slice",
 102      &  i2/80(1h-)//)') islice
 103       write (iout,*) "delta",delta," nbin1",nbin1
 104       write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
 105       call flush(iout)
 106       dmin=0.0d0
 107       tmax=0
 108       do i=1,nParmset
 109         do j=1,nT_h(i)
 110           potEmin_all(j,i)=1.0d10
 111         enddo
 112       enddo
 113       rgymin=1.0d10
 114       rmsmin=1.0d10
 115       rgymax=0.0d0
 116       rmsmax=0.0d0
 117       do t=0,MaxN
 118         htot(t)=0
 119       enddo
 120 #ifdef MPI
 121       do i=1,scount(me1)
 122 #else
 123       do i=1,ntot(islice)
 124 #endif
 125         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
 126         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
 127         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
 128         if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
 129         ind_point(i)=0
 130         do j=nQ,1,-1
 131           ind=(q(j,i)-dmin+1.0d-8)/delta
 132           if (j.eq.1) then
 133             ind_point(i)=ind_point(i)+ind
 134           else
 135             ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
 136           endif
 137           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
 138             write (iout,*) "Error - index exceeds range for point",i,
 139      &      " q=",q(j,i)," ind",ind_point(i)
 140 #ifdef MPI
 141             write (iout,*) "Processor",me1
 142             call flush(iout)
 143             call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
 144 #endif
 145             stop
 146           endif
 147         enddo ! j
 148         if (ind_point(i).gt.tmax) tmax=ind_point(i)
 149         htot(ind_point(i))=htot(ind_point(i))+1
 150 #ifdef DEBUG
 151         write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),
 152      &   " htot",htot(ind_point(i))
 153         call flush(iout)
 154 #endif
 155       enddo ! i
 156       call flush(iout)
 157
 158       nbin=nbin1**nQ-1
 159       write (iout,'(a)') "Numbers of counts in Q bins"
 160       do t=0,tmax
 161         if (htot(t).gt.0) then
 162         write (iout,'(i15,$)') t
 163         liczba=t
 164         do j=1,nQ
 165           jj = mod(liczba,nbin1)
 166           liczba=liczba/nbin1
 167           write (iout,'(i5,$)') jj
 168         enddo
 169         write (iout,'(i8)') htot(t)
 170         endif
 171       enddo
 172       do iparm=1,nParmSet
 173       write (iout,'(a,i3)') "Number of data points for parameter set",
 174      & iparm
 175       write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),
 176      & ib=1,nT_h(iparm))
 177       write (iout,'(i8)') stot(islice)
 178       write (iout,'(a)')
 179       enddo
 180       call flush(iout)
 181
 182 #ifdef MPI
 183       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
 184      &  WHAM_COMM,IERROR)
 185       tmax=tmax_t
 186       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
 187      &  MPI_MIN,WHAM_COMM,IERROR)
 188       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
 189      &  MPI_MAX,WHAM_COMM,IERROR)
 190       call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,
 191      &  MPI_MIN,WHAM_COMM,IERROR)
 192       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
 193      &  MPI_MAX,WHAM_COMM,IERROR)
 194       rgymin=rgymin_t
 195       rgymax=rgymax_t
 196       rmsmin=rmsmin_t
 197       rmsmax=rmsmax_t
 198 #endif
 199       rmsmin=deltrms*dint(rmsmin/deltrms)
 200       rmsmax=deltrms*dint(rmsmax/deltrms)
 201       rgymin=deltrms*dint(rgymin/deltrgy)
 202       rgymax=deltrms*dint(rgymax/deltrgy)
 203       nbin_rms=(rmsmax-rmsmin)/deltrms
 204       nbin_rgy=(rgymax-rgymin)/deltrgy
 205       write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,
 206      & " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
 207       nFi=0
 208       do i=1,nParmSet
 209         do j=1,nT_h(i)
 210           nFi=nFi+nR(j,i)
 211         enddo
 212       enddo
 213       write (iout,*) "nFi",nFi
 214 ! Compute the Boltzmann factor corresponing to restrain potentials in different
 215 ! simulations.
 216 #ifdef MPI
 217       do i=1,scount(me1)
 218 #else
 219       do i=1,ntot(islice)
 220 #endif
 221 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 222         do iparm=1,nParmSet
 223 #ifdef DEBUG
 224           write (iout,'(2i5,22f8.2)') i,iparm,
 225      &     (enetb(k,i,iparm),k=1,22)
 226 #endif
 227           call restore_parm(iparm)
 228 #ifdef DEBUG
 229           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 230      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 231      &      wtor_d,wsccor,wbond
 232 #endif
 233           do ib=1,nT_h(iparm)
 234             if (rescale_mode.eq.1) then
 235               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 236               quotl=1.0d0
 237               kfacl=1.0d0
 238               do l=1,5
 239                 quotl1=quotl
 240                 quotl=quotl*quot
 241                 kfacl=kfacl*kfac
 242                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 243               enddo
 244 #if defined(FUNCTH)
 245               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 246               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 247 #elif defined(FUNCT)
 248               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 249 #else
 250               ft(6)=1.0d0
 251 #endif
 252             else if (rescale_mode.eq.2) then
 253               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 254               quotl=1.0d0
 255               do l=1,5
 256                 quotl=quotl*quot
 257                 fT(l)=1.12692801104297249644d0/
 258      &             dlog(dexp(quotl)+dexp(-quotl))
 259               enddo
 260 #if defined(FUNCTH)
 261               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 262               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 263 #elif defined(FUNCT)
 264               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 265 #else
 266               ft(6)=1.0d0
 267 #endif
 268 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 269             else if (rescale_mode.eq.0) then
 270               do l=1,6
 271                 fT(l)=1.0d0
 272               enddo
 273             else
 274               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 275      &          rescale_mode
 276               call flush(iout)
 277               return1
 278             endif
 279             evdw=enetb(1,i,iparm)
 280             evdw_t=enetb(21,i,iparm)
 281 #ifdef SCP14
 282             evdw2_14=enetb(17,i,iparm)
 283             evdw2=enetb(2,i,iparm)+evdw2_14
 284 #else
 285             evdw2=enetb(2,i,iparm)
 286             evdw2_14=0.0d0
 287 #endif
 288 #ifdef SPLITELE
 289             ees=enetb(3,i,iparm)
 290             evdw1=enetb(16,i,iparm)
 291 #else
 292             ees=enetb(3,i,iparm)
 293             evdw1=0.0d0
 294 #endif
 295             ecorr=enetb(4,i,iparm)
 296             ecorr5=enetb(5,i,iparm)
 297             ecorr6=enetb(6,i,iparm)
 298             eel_loc=enetb(7,i,iparm)
 299             eello_turn3=enetb(8,i,iparm)
 300             eello_turn4=enetb(9,i,iparm)
 301             eturn6=enetb(10,i,iparm)
 302             ebe=enetb(11,i,iparm)
 303             escloc=enetb(12,i,iparm)
 304             etors=enetb(13,i,iparm)
 305             etors_d=enetb(14,i,iparm)
 306             ehpb=enetb(15,i,iparm)
 307             estr=enetb(18,i,iparm)
 308             esccor=enetb(19,i,iparm)
 309             edihcnstr=enetb(20,i,iparm)
 310             ehomology_constr=enetb(22,i,iparm)
 311             edfadis=enetb(23,i,iparm)
 312             edfator=enetb(24,i,iparm)
 313             edfanei=enetb(25,i,iparm)
 314             edfabet=enetb(26,i,iparm)
 315 #ifdef DEBUG
 316             write (iout,'(3i5,6f5.2,16f12.3)') i,ib,iparm,(ft(l),l=1,6),
 317      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
 318      &       etors,etors_d,eello_turn3,eello_turn4,esccor,
 319      &       ehomology_constr,edfadis,edfator,edfanei,edfabet
 320 #endif
 321
 322 #ifdef SPLITELE
 323             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 324      &      +wvdwpp*evdw1
 325      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 326      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 327      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 328      &      +ft(2)*wturn3*eello_turn3
 329      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 330      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 331      &      +wbond*estr+wdfa_dist*edfadis
 332      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 333 #else
 334             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 335      &      +ft(1)*welec*(ees+evdw1)
 336      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 337      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 338      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 339      &      +ft(2)*wturn3*eello_turn3
 340      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 341      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 342      &      +wbond*estr+wdfa_dist*edfadis
 343      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 344 #endif
 345 #ifdef DEBUG
 346             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
 347      &        etot
 348 #endif
 349 #ifdef DEBUG
 350             if (iparm.eq.1 .and. ib.eq.1) then
 351             write (iout,*)"Conformation",i
 352             energia(0)=etot
 353             do k=1,max_ene
 354               energia(k)=enetb(k,i,iparm)
 355             enddo
 356             call enerprint(energia(0),fT)
 357             endif
 358 #endif
 359 #ifdef DEBUG
 360             write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
 361 #endif
 362             if (homol_nset.gt.1) then
 363
 364             do kk=1,nR(ib,iparm)
 365               Econstr=waga_homology(kk)*ehomology_constr
 366               v(i,kk,ib,iparm)=
 367      &          -beta_h(ib,iparm)*(etot+Econstr)
 368 #ifdef DEBUG
 369               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
 370      &         etot,Econstr,v(i,kk,ib,iparm)
 371 #endif
 372             enddo ! kk
 373
 374             else
 375
 376             etot=etot+ehomology_constr
 377             do kk=1,nR(ib,iparm)
 378               Econstr=0.0d0
 379               do j=1,nQ
 380                 dd = q(j,i)
 381                 Econstr=Econstr+Kh(j,kk,ib,iparm)
 382      &           *(dd-q0(j,kk,ib,iparm))**2
 383               enddo
 384               v(i,kk,ib,iparm)=
 385      &          -beta_h(ib,iparm)*(etot+Econstr)
 386 #ifdef DEBUG
 387               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
 388      &         etot,v(i,kk,ib,iparm)
 389 #endif
 390             enddo ! kk
 391
 392             endif
 393           enddo   ! ib
 394         enddo     ! iparm
 395       enddo       ! i
 396 ! Simple iteration to calculate free energies corresponding to all simulation
 397 ! runs.
 398       do iter=1,maxit
 399
 400 ! Compute new free-energy values corresponding to the righ-hand side of the
 401 ! equation and their derivatives.
 402         write (iout,*) "------------------------fi"
 403         entfac_min=1.0d10
 404 #ifdef MPI
 405         do t=1,scount(me1)
 406 #else
 407         do t=1,ntot(islice)
 408 #endif
 409           vmax=-1.0d+20
 410           do i=1,nParmSet
 411             do k=1,nT_h(i)
 412               do l=1,nR(k,i)
 413                 vf=v(t,l,k,i)+f(l,k,i)
 414                 if (vf.gt.vmax) vmax=vf
 415               enddo
 416             enddo
 417           enddo
 418           denom=0.0d0
 419           do i=1,nParmSet
 420             do k=1,nT_h(i)
 421               do l=1,nR(k,i)
 422                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 423                 if (aux.gt.-200.0d0)
 424      &          denom=denom+snk(l,k,i,islice)*dexp(aux)
 425               enddo
 426             enddo
 427           enddo
 428           entfac(t)=-dlog(denom)-vmax
 429           if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 430 #ifdef DEBUG
 431           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 432 #endif
 433         enddo
 434 c#ifdef MPI
 435 c        write (iout,*) "entfac_min before AllReduce",entfac_min
 436 c        call MPI_AllReduce(entfac_min,entfac_min_t,1,
 437 c     &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 438 c        entfac_min=entfac_min_t
 439 c        write (iout,*) "entfac_min after AllReduce",entfac_min
 440 c#endif
 441 c#ifdef MPI
 442 c        do t=1,scount(me)
 443 c          entfac(t)=entfac(t)-entfac_min
 444 c        enddo
 445 c#else
 446 c        do t=1,ntot(islice)
 447 c          entfac(t)=entfac(t)-entfac_min
 448 c        enddo
 449 c#endif
 450         do iparm=1,nParmSet
 451           do iib=1,nT_h(iparm)
 452             do ii=1,nR(iib,iparm)
 453 #ifdef MPI
 454               fi_p_min(ii,iib,iparm)=-1.0d10
 455               do t=1,scount(me)
 456                 aux=v(t,ii,iib,iparm)+entfac(t)
 457                 if (aux.gt.fi_p_min(ii,iib,iparm))
 458      &                                   fi_p_min(ii,iib,iparm)=aux
 459               enddo
 460 #else
 461               do t=1,ntot(islice)
 462                 aux=v(t,ii,iib,iparm)+entfac(t)
 463                 if (aux.gt.fi_min(ii,iib,iparm))
 464      &                                     fi_min(ii,iib,iparm)=aux
 465               enddo
 466 #endif
 467             enddo ! ii
 468           enddo ! iib
 469         enddo ! iparm
 470 #ifdef MPI
 471 #ifdef DEBUG
 472         write (iout,*) "fi_min before AllReduce"
 473         do i=1,nParmSet
 474           do j=1,nT_h(i)
 475             write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
 476           enddo
 477         enddo
 478 #endif
 479         call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
 480      &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 481 #ifdef DEBUG
 482         write (iout,*) "fi_min after AllReduce"
 483         do i=1,nParmSet
 484           do j=1,nT_h(i)
 485             write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
 486           enddo
 487         enddo
 488 #endif
 489 #endif
 490         do iparm=1,nParmSet
 491           do iib=1,nT_h(iparm)
 492             do ii=1,nR(iib,iparm)
 493 #ifdef MPI
 494               fi_p(ii,iib,iparm)=0.0d0
 495               do t=1,scount(me)
 496                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
 497      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 498 #ifdef DEBUG
 499               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
 500      &         v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
 501      &         fi_p(ii,iib,iparm)
 502 #endif
 503               enddo
 504 #else
 505               fi(ii,iib,iparm)=0.0d0
 506               do t=1,ntot(islice)
 507                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
 508      &           +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 509               enddo
 510 #endif
 511             enddo ! ii
 512           enddo ! iib
 513         enddo ! iparm
 514
 515 #ifdef MPI
 516 #ifdef DEBUG
 517         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 518         do iparm=1,nParmSet
 519           do ib=1,nT_h(iparm)
 520             write (iout,*) "iparm",iparm," ib",ib
 521             write (iout,*) "beta=",beta_h(ib,iparm)
 522             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 523           enddo
 524         enddo
 525 #endif
 526 #ifdef DEBUG
 527         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
 528      &   maxR*MaxT_h*nParmSet
 529         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
 530      &   " WHAM_COMM",WHAM_COMM
 531 #endif
 532         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
 533      &   MPI_DOUBLE_PRECISION,
 534      &   MPI_SUM,Master,WHAM_COMM,IERROR)
 535 #ifdef DEBUG
 536         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 537         do iparm=1,nParmSet
 538           write (iout,*) "iparm",iparm
 539           do ib=1,nT_h(iparm)
 540             write (iout,*) "beta=",beta_h(ib,iparm)
 541             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 542           enddo
 543         enddo
 544 #endif
 545         if (me1.eq.Master) then
 546 #endif
 547         avefi=0.0d0
 548         do iparm=1,nParmSet
 549           do ib=1,nT_h(iparm)
 550             do i=1,nR(ib,iparm)
 551               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
 552               avefi=avefi+fi(i,ib,iparm)
 553             enddo
 554           enddo
 555         enddo
 556         avefi=avefi/nFi
 557         do iparm=1,nParmSet
 558           write (iout,*) "Parameter set",iparm
 559           do ib =1,nT_h(iparm)
 560             write (iout,*) "beta=",beta_h(ib,iparm)
 561             do i=1,nR(ib,iparm)
 562               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 563             enddo
 564             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 565             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 566           enddo
 567         enddo
 568
 569 ! Compute the norm of free-energy increments.
 570         finorm=0.0d0
 571         do iparm=1,nParmSet
 572           do ib=1,nT_h(iparm)
 573             do i=1,nR(ib,iparm)
 574               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 575               f(i,ib,iparm)=fi(i,ib,iparm)
 576             enddo
 577           enddo
 578         enddo
 579
 580         write (iout,*) 'Iteration',iter,' finorm',finorm
 581
 582 #ifdef MPI
 583         endif
 584         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,
 585      &   MPI_DOUBLE_PRECISION,Master,
 586      &   WHAM_COMM,IERROR)
 587         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,
 588      &   WHAM_COMM,IERROR)
 589 #endif
 590 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 591         if (finorm.lt.fimin) then
 592           write (iout,*) 'Iteration converged'
 593           goto 20
 594         endif
 595
 596       enddo ! iter
 597
 598    20 continue
 599 ! Now, put together the histograms from all simulations, in order to get the
 600 ! unbiased total histogram.
 601
 602 C Determine the minimum free energies
 603 #ifdef MPI
 604       do i=1,scount(me1)
 605 #else
 606       do i=1,ntot(islice)
 607 #endif
 608 c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 609         do iparm=1,nParmSet
 610 #ifdef DEBUG
 611           write (iout,'(2i5,21f8.2)') i,iparm,
 612      &     (enetb(k,i,iparm),k=1,22)
 613 #endif
 614           call restore_parm(iparm)
 615 #ifdef DEBUG
 616           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
 617      &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
 618      &      wtor_d,wsccor,wbond
 619 #endif
 620           do ib=1,nT_h(iparm)
 621             if (rescale_mode.eq.1) then
 622               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 623               quotl=1.0d0
 624               kfacl=1.0d0
 625               do l=1,5
 626                 quotl1=quotl
 627                 quotl=quotl*quot
 628                 kfacl=kfacl*kfac
 629                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 630               enddo
 631 #if defined(FUNCTH)
 632               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 633               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 634 #elif defined(FUNCT)
 635               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 636 #else
 637               ft(6)=1.0d0
 638 #endif
 639             else if (rescale_mode.eq.2) then
 640               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 641               quotl=1.0d0
 642               do l=1,5
 643                 quotl=quotl*quot
 644                 fT(l)=1.12692801104297249644d0/
 645      &             dlog(dexp(quotl)+dexp(-quotl))
 646               enddo
 647 #if defined(FUNCTH)
 648               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 649               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 650 #elif defined(FUNCT)
 651               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 652 #else
 653               ft(6)=1.0d0
 654 #endif
 655 c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 656             else if (rescale_mode.eq.0) then
 657               do l=1,6
 658                 fT(l)=1.0d0
 659               enddo
 660             else
 661               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 662      &          rescale_mode
 663               call flush(iout)
 664               return1
 665             endif
 666             evdw=enetb(1,i,iparm)
 667             evdw_t=enetb(21,i,iparm)
 668 #ifdef SCP14
 669             evdw2_14=enetb(17,i,iparm)
 670             evdw2=enetb(2,i,iparm)+evdw2_14
 671 #else
 672             evdw2=enetb(2,i,iparm)
 673             evdw2_14=0.0d0
 674 #endif
 675 #ifdef SPLITELE
 676             ees=enetb(3,i,iparm)
 677             evdw1=enetb(16,i,iparm)
 678 #else
 679             ees=enetb(3,i,iparm)
 680             evdw1=0.0d0
 681 #endif
 682             ecorr=enetb(4,i,iparm)
 683             ecorr5=enetb(5,i,iparm)
 684             ecorr6=enetb(6,i,iparm)
 685             eel_loc=enetb(7,i,iparm)
 686             eello_turn3=enetb(8,i,iparm)
 687             eello_turn4=enetb(9,i,iparm)
 688             eturn6=enetb(10,i,iparm)
 689             ebe=enetb(11,i,iparm)
 690             escloc=enetb(12,i,iparm)
 691             etors=enetb(13,i,iparm)
 692             etors_d=enetb(14,i,iparm)
 693             ehpb=enetb(15,i,iparm)
 694             estr=enetb(18,i,iparm)
 695             esccor=enetb(19,i,iparm)
 696             edihcnstr=enetb(20,i,iparm)
 697             ehomology_constr=enetb(22,i,iparm)
 698             if (homol_nset.gt.1)
 699      &       ehomology_constr=waga_homology(homol_nset)*ehomology_constr
 700             edfadis=enetb(23,i,iparm)
 701             edfator=enetb(24,i,iparm)
 702             edfanei=enetb(25,i,iparm)
 703             edfabet=enetb(26,i,iparm)
 704 #ifdef DEBUG
 705             write (iout,'(3i5,6f5.2,17f12.3)') i,ib,iparm,(ft(l),l=1,6),
 706      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
 707      &       etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr,
 708      &       ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
 709      &       wdfa_nei*edfanei+wdfa_beta*edfabet
 710 #endif
 711
 712 #ifdef SPLITELE
 713             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
 714      &      +wvdwpp*evdw1
 715      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 716      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 717      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 718      &      +ft(2)*wturn3*eello_turn3
 719      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
 720      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 721      &      +wbond*estr+ehomology_constr+wdfa_dist*edfadis
 722      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 723 #else
 724             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
 725      &      +ft(1)*welec*(ees+evdw1)
 726      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
 727      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
 728      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
 729      &      +ft(2)*wturn3*eello_turn3
 730      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
 731      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
 732      &      +wbond*estr+ehomology_constr+wdfa_dist*edfadis
 733      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 734 #endif
 735 #ifdef DEBUG
 736             write (iout,*) "i",i," ib",ib,
 737      &      " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
 738      &      " entfac",entfac(i)," ecorr",etot-entfac(i)/beta_h(ib,iparm)
 739             etot=etot-entfac(i)/beta_h(ib,iparm)
 740 #endif
 741             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 742 #ifdef DEBUG
 743             write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
 744 #endif
 745           enddo   ! ib
 746         enddo     ! iparm
 747       enddo       ! i
 748 #ifdef DEBUG
 749       write (iout,*) "The potEmin array before reduction"
 750       do i=1,nParmSet
 751         write (iout,*) "Parameter set",i
 752         do j=1,nT_h(i)
 753           write (iout,*) j,PotEmin_all(j,i)
 754         enddo
 755       enddo
 756       write (iout,*) "potEmin_min",potEmin_min
 757 #endif
 758 #ifdef MPI
 759 C Determine the minimum energes for all parameter sets and temperatures
 760       call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
 761      &  maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
 762       do i=1,nParmSet
 763         do j=1,nT_h(i)
 764           potEmin_all(j,i)=potEmin_t_all(j,i)
 765         enddo
 766       enddo
 767 #endif
 768       potEmin_min=potEmin_all(1,1)
 769       do i=1,nParmSet
 770         do j=1,nT_h(i)
 771           if (potEmin_all(j,i).lt.potEmin_min)
 772      &             potEmin_min=potEmin_all(j,i)
 773         enddo
 774       enddo
 775 #ifdef DEBUG
 776       write (iout,*) "The potEmin array"
 777       do i=1,nParmSet
 778         write (iout,*) "Parameter set",i
 779         do j=1,nT_h(i)
 780           write (iout,*) j,PotEmin_all(j,i)
 781         enddo
 782       enddo
 783       write (iout,*) "potEmin_min",potEmin_min
 784 #endif
 785
 786 #ifdef MPI
 787       do t=0,tmax
 788         hfin_ent_p(t)=0.0d0
 789       enddo
 790 #else
 791       do t=0,tmax
 792         hfin_ent(t)=0.0d0
 793       enddo
 794 #endif
 795       write (iout,*) "--------------hist"
 796 #ifdef MPI
 797       do iparm=1,nParmSet
 798         do i=0,nGridT
 799           sumW_p(i,iparm)=0.0d0
 800           sumE_p(i,iparm)=0.0d0
 801           sumEbis_p(i,iparm)=0.0d0
 802           sumEsq_p(i,iparm)=0.0d0
 803           do j=1,nQ+2
 804             sumQ_p(j,i,iparm)=0.0d0
 805             sumQsq_p(j,i,iparm)=0.0d0
 806             sumEQ_p(j,i,iparm)=0.0d0
 807           enddo
 808         enddo
 809       enddo
 810       upindE_p=0
 811 #else
 812       do iparm=1,nParmSet
 813         do i=0,nGridT
 814           sumW(i,iparm)=0.0d0
 815           sumE(i,iparm)=0.0d0
 816           sumEbis(i,iparm)=0.0d0
 817           sumEsq(i,iparm)=0.0d0
 818           do j=1,nQ+2
 819             sumQ(j,i,iparm)=0.0d0
 820             sumQsq(j,i,iparm)=0.0d0
 821             sumEQ(j,i,iparm)=0.0d0
 822           enddo
 823         enddo
 824       enddo
 825       upindE=0
 826 #endif
 827 c 8/26/05 entropy distribution
 828 #ifdef MPI
 829       entmin_p=1.0d10
 830       entmax_p=-1.0d10
 831       do t=1,scount(me1)
 832 c        ent=-dlog(entfac(t))
 833         ent=entfac(t)
 834         if (ent.lt.entmin_p) entmin_p=ent
 835         if (ent.gt.entmax_p) entmax_p=ent
 836       enddo
 837       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 838       call flush(iout)
 839       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 840      &  WHAM_COMM,IERROR)
 841       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 842      &  WHAM_COMM,IERROR)
 843       ientmax=entmax-entmin
 844       if (ientmax.gt.2000) ientmax=2000
 845       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 846       call flush(iout)
 847       do t=1,scount(me1)
 848 c        ient=-dlog(entfac(t))-entmin
 849         ient=entfac(t)-entmin
 850         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 851       enddo
 852       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 853      &  MPI_SUM,WHAM_COMM,IERROR)
 854       if (me1.eq.Master) then
 855         write (iout,*) "Entropy histogram"
 856         do i=0,ientmax
 857           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 858         enddo
 859       endif
 860 #else
 861       entmin=1.0d10
 862       entmax=-1.0d10
 863       do t=1,ntot(islice)
 864         ent=entfac(t)
 865         if (ent.lt.entmin) entmin=ent
 866         if (ent.gt.entmax) entmax=ent
 867       enddo
 868       ientmax=-dlog(entmax)-entmin
 869       if (ientmax.gt.2000) ientmax=2000
 870       do t=1,ntot(islice)
 871         ient=entfac(t)-entmin
 872         if (ient.le.2000) histent(ient)=histent(ient)+1
 873       enddo
 874       write (iout,*) "Entropy histogram"
 875       do i=0,ientmax
 876         write(iout,'(2f15.4)') entmin+i,histent(i)
 877       enddo
 878 #endif
 879
 880 #ifdef MPI
 881 c      write (iout,*) "me1",me1," scount",scount(me1)
 882
 883       do iparm=1,nParmSet
 884
 885 #ifdef MPI
 886         do ib=1,nT_h(iparm)
 887           do t=0,tmax
 888             hfin_p(t,ib)=0.0d0
 889           enddo
 890         enddo
 891         do i=1,maxindE
 892           histE_p(i)=0.0d0
 893         enddo
 894 #else
 895         do ib=1,nT_h(iparm)
 896           do t=0,tmax
 897             hfin(t,ib)=0.0d0
 898           enddo
 899         enddo
 900         do i=1,maxindE
 901           histE(i)=0.0d0
 902         enddo
 903 #endif
 904         do ib=1,nT_h(iparm)
 905           do i=0,MaxBinRms
 906             do j=0,MaxBinRgy
 907               hrmsrgy(j,i,ib)=0.0d0
 908 #ifdef MPI
 909               hrmsrgy_p(j,i,ib)=0.0d0
 910 #endif
 911             enddo
 912           enddo
 913         enddo
 914
 915         do t=1,scount(me1)
 916 #else
 917         do t=1,ntot(islice)
 918 #endif
 919           ind = ind_point(t)
 920 #ifdef MPI
 921           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 922 #else
 923           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 924 #endif
 925           call restore_parm(iparm)
 926           evdw=enetb(21,t,iparm)
 927           evdw_t=enetb(1,t,iparm)
 928 #ifdef SCP14
 929           evdw2_14=enetb(17,t,iparm)
 930           evdw2=enetb(2,t,iparm)+evdw2_14
 931 #else
 932           evdw2=enetb(2,t,iparm)
 933           evdw2_14=0.0d0
 934 #endif
 935 #ifdef SPLITELE
 936           ees=enetb(3,t,iparm)
 937           evdw1=enetb(16,t,iparm)
 938 #else
 939           ees=enetb(3,t,iparm)
 940           evdw1=0.0d0
 941 #endif
 942           ecorr=enetb(4,t,iparm)
 943           ecorr5=enetb(5,t,iparm)
 944           ecorr6=enetb(6,t,iparm)
 945           eel_loc=enetb(7,t,iparm)
 946           eello_turn3=enetb(8,t,iparm)
 947           eello_turn4=enetb(9,t,iparm)
 948           eturn6=enetb(10,t,iparm)
 949           ebe=enetb(11,t,iparm)
 950           escloc=enetb(12,t,iparm)
 951           etors=enetb(13,t,iparm)
 952           etors_d=enetb(14,t,iparm)
 953           ehpb=enetb(15,t,iparm)
 954           estr=enetb(18,t,iparm)
 955           esccor=enetb(19,t,iparm)
 956           edihcnstr=enetb(20,t,iparm)
 957           ehomology_constr=enetb(22,t,iparm)
 958           if (homol_nset.gt.1)
 959      &       ehomology_constr=waga_homology(homol_nset)*ehomology_constr
 960           edfadis=enetb(23,t,iparm)
 961           edfator=enetb(24,t,iparm)
 962           edfanei=enetb(25,t,iparm)
 963           edfabet=enetb(26,t,iparm)
 964           do k=0,nGridT
 965             betaT=startGridT+k*delta_T
 966             temper=betaT
 967 c            fT=T0/betaT
 968 c            ft=2*T0/(T0+betaT)
 969             if (rescale_mode.eq.1) then
 970               quot=betaT/T0
 971               quotl=1.0d0
 972               kfacl=1.0d0
 973               do l=1,5
 974                 quotl1=quotl
 975                 quotl=quotl*quot
 976                 kfacl=kfacl*kfac
 977                 denom=kfacl-1.0d0+quotl
 978                 fT(l)=kfacl/denom
 979                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 980                 ftbis(l)=l*kfacl*quotl1*
 981      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 982               enddo
 983 #if defined(FUNCTH)
 984               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 985      &                  320.0d0
 986               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 987              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 988      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 989 #elif defined(FUNCT)
 990               fT(6)=betaT/T0
 991               ftprim(6)=1.0d0/T0
 992               ftbis(6)=0.0d0
 993 #else
 994               fT(6)=1.0d0
 995               ftprim(6)=0.0d0
 996               ftbis(6)=0.0d0
 997 #endif
 998             else if (rescale_mode.eq.2) then
 999               quot=betaT/T0
1000               quotl=1.0d0
1001               do l=1,5
1002                 quotl1=quotl
1003                 quotl=quotl*quot
1004                 eplus=dexp(quotl)
1005                 eminus=dexp(-quotl)
1006                 logfac=1.0d0/dlog(eplus+eminus)
1007                 tanhT=(eplus-eminus)/(eplus+eminus)
1008                 fT(l)=1.12692801104297249644d0*logfac
1009                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1010                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
1011      &          2*l*quotl1/T0*logfac*
1012      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
1013      &          +ftprim(l)*tanhT)
1014               enddo
1015 #if defined(FUNCTH)
1016               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
1017      &                 320.0d0
1018               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1019              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
1020      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1021 #elif defined(FUNCT)
1022               fT(6)=betaT/T0
1023               ftprim(6)=1.0d0/T0
1024               ftbis(6)=0.0d0
1025 #else
1026               fT(6)=1.0d0
1027               ftprim(6)=0.0d0
1028               ftbis(6)=0.0d0
1029 #endif
1030             else if (rescale_mode.eq.0) then
1031               do l=1,5
1032                 fT(l)=1.0d0
1033                 ftprim(l)=0.0d0
1034               enddo
1035             else
1036               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
1037      &          rescale_mode
1038               call flush(iout)
1039               return1
1040             endif
1041 c            write (iout,*) "ftprim",ftprim
1042 c            write (iout,*) "ftbis",ftbis
1043             betaT=1.0d0/(1.987D-3*betaT)
1044             if (betaT.ge.beta_h(1,iparm)) then
1045               potEmin=potEmin_all(1,iparm)
1046 c              write(iout,*) "first",temper,potEmin
1047             else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1048               potEmin=potEmin_all(nT_h(iparm),iparm)
1049 c              write (iout,*) "last",temper,potEmin
1050             else
1051               do l=1,nT_h(iparm)-1
1052                 if (betaT.le.beta_h(l,iparm) .and.
1053      &              betaT.gt.beta_h(l+1,iparm)) then
1054                   potEmin=potEmin_all(l,iparm)
1055 c                  write (iout,*) "l",l,
1056 c     &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1057 c     &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1058                   exit
1059                 endif
1060               enddo
1061             endif
1062 c            write (iout,*) ib," PotEmin",potEmin
1063 #ifdef SPLITELE
1064             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1065      &      +wvdwpp*evdw1
1066      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1067      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1068      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1069      &      +ft(2)*wturn3*eello_turn3
1070      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1071      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1072      &      +wbond*estr+ehomology_constr
1073             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1074      &            +ftprim(1)*wtor*etors+
1075      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1076      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1077      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1078      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1079      &            ftprim(1)*wsccor*esccor
1080             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1081      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1082      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1083      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1084      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1085      &            ftbis(1)*wsccor*esccor
1086 #else
1087             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1088      &      +ft(1)*welec*(ees+evdw1)
1089      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1090      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1091      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1092      &      +ft(2)*wturn3*eello_turn3
1093      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1094      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1095      &      +wbond*estr+ehomology_constr
1096             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1097      &           +ftprim(1)*wtor*etors+
1098      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1099      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1100      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1101      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1102      &            ftprim(1)*wsccor*esccor
1103             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1104      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1105      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1106      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1107      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1108      &            ftprim(1)*wsccor*esccor
1109 #endif
1110             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1111 #ifdef DEBUG
1112             write (iout,'(3i5,6f5.2,17f12.3)') i,ib,iparm,(ft(l),l=1,6),
1113      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
1114      &       etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr,
1115      &       ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
1116      &       wdfa_nei*edfanei+wdfa_beta*edfabet
1117             write (iout,*) "iparm",iparm," t",t," temper",temper,
1118      &        " etot",etot," entfac",entfac(t),
1119      &        " efree",etot-entfac(t)/betaT," potEmin",potEmin,
1120      &        " boltz",-betaT*(etot-potEmin)+entfac(t),
1121      &        " weight",weight," ebis",ebis
1122 #endif
1123             etot=etot-temper*eprim
1124 #ifdef MPI
1125             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1126             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1127             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1128             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1129             do j=1,nQ+2
1130               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1131               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1132               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1133      &         +etot*q(j,t)*weight
1134             enddo
1135 #else
1136             sumW(k,iparm)=sumW(k,iparm)+weight
1137             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1138             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1139             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1140             do j=1,nQ+2
1141               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1142               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1143               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1144      &         +etot*q(j,t)*weight
1145             enddo
1146 #endif
1147           enddo
1148           indE = aint(potE(t,iparm)-aint(potEmin))
1149           if (indE.ge.0 .and. indE.le.maxinde) then
1150             if (indE.gt.upindE_p) upindE_p=indE
1151             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1152           endif
1153 #ifdef MPI
1154           do ib=1,nT_h(iparm)
1155             potEmin=potEmin_all(ib,iparm)
1156             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1157             hfin_p(ind,ib)=hfin_p(ind,ib)+
1158      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1159              if (rmsrgymap) then
1160                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1161                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1162                hrmsrgy_p(indrgy,indrms,ib)=
1163      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1164              endif
1165           enddo
1166 #else
1167           do ib=1,nT_h(iparm)
1168             potEmin=potEmin_all(ib,iparm)
1169             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1170             hfin(ind,ib)=hfin(ind,ib)+
1171      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1172              if (rmsrgymap) then
1173                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1174                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1175                hrmsrgy(indrgy,indrms,ib)=
1176      &           hrmsrgy(indrgy,indrms,ib)+expfac
1177              endif
1178           enddo
1179 #endif
1180         enddo ! t
1181         do ib=1,nT_h(iparm)
1182           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1183      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1184           if (rmsrgymap) then
1185           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1186      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1187      &       WHAM_COMM,IERROR)
1188           endif
1189         enddo
1190         call MPI_Reduce(upindE_p,upindE,1,
1191      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1192         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1193      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1194
1195         if (me1.eq.master) then
1196
1197         if (histout) then
1198
1199         write (iout,'(6x,$)')
1200         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1201      &   ib=1,nT_h(iparm))
1202         write (iout,*)
1203
1204         write (iout,'(/a)') 'Final histograms'
1205         if (histfile) then
1206           if (nslice.eq.1) then
1207             if (separate_parset) then
1208               write(licz3,"(bz,i3.3)") myparm
1209               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1210             else
1211               histname=prefix(:ilen(prefix))//'.hist'
1212             endif
1213           else
1214             if (separate_parset) then
1215               write(licz3,"(bz,i3.3)") myparm
1216               histname=prefix(:ilen(prefix))//'_par'//licz3//
1217      &         '_slice_'//licz2//'.hist'
1218             else
1219               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1220             endif
1221           endif
1222 #if defined(AIX) || defined(PGI)
1223           open (ihist,file=histname,position='append')
1224 #else
1225           open (ihist,file=histname,access='append')
1226 #endif
1227         endif
1228
1229         do t=0,tmax
1230           liczba=t
1231           sumH=0.0d0
1232           do ib=1,nT_h(iparm)
1233             sumH=sumH+hfin(t,ib)
1234           enddo
1235           if (sumH.gt.0.0d0) then
1236             do j=1,nQ
1237               jj = mod(liczba,nbin1)
1238               liczba=liczba/nbin1
1239               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1240               if (histfile)
1241      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1242             enddo
1243             do ib=1,nT_h(iparm)
1244               write (iout,'(e20.10,$)') hfin(t,ib)
1245               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1246             enddo
1247             write (iout,'(i5)') iparm
1248             if (histfile) write (ihist,'(i5)') iparm
1249           endif
1250         enddo
1251
1252         endif
1253
1254         if (entfile) then
1255           if (nslice.eq.1) then
1256             if (separate_parset) then
1257               write(licz3,"(bz,i3.3)") myparm
1258               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1259             else
1260               histname=prefix(:ilen(prefix))//'.ent'
1261             endif
1262           else
1263             if (separate_parset) then
1264               write(licz3,"(bz,i3.3)") myparm
1265               histname=prefix(:ilen(prefix))//'par_'//licz3//
1266      &           '_slice_'//licz2//'.ent'
1267             else
1268               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1269             endif
1270           endif
1271 #if defined(AIX) || defined(PGI)
1272           open (ihist,file=histname,position='append')
1273 #else
1274           open (ihist,file=histname,access='append')
1275 #endif
1276           write (ihist,'(a)') "# Microcanonical entropy"
1277           do i=0,upindE
1278             write (ihist,'(f8.0,$)') dint(potEmin)+i
1279             if (histE(i).gt.0.0e0) then
1280               write (ihist,'(f15.5,$)') dlog(histE(i))
1281             else
1282               write (ihist,'(f15.5,$)') 0.0d0
1283             endif
1284           enddo
1285           write (ihist,*)
1286           close(ihist)
1287         endif
1288         write (iout,*) "Microcanonical entropy"
1289         do i=0,upindE
1290           write (iout,'(f8.0,$)') dint(potEmin)+i
1291           if (histE(i).gt.0.0e0) then
1292             write (iout,'(f15.5,$)') dlog(histE(i))
1293           else
1294             write (iout,'(f15.5,$)') 0.0d0
1295           endif
1296           write (iout,*)
1297         enddo
1298         if (rmsrgymap) then
1299           if (nslice.eq.1) then
1300             if (separate_parset) then
1301               write(licz3,"(bz,i3.3)") myparm
1302               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1303             else
1304               histname=prefix(:ilen(prefix))//'.rmsrgy'
1305             endif
1306           else
1307             if (separate_parset) then
1308               write(licz3,"(bz,i3.3)") myparm
1309               histname=prefix(:ilen(prefix))//'_par'//licz3//
1310      &         '_slice_'//licz2//'.rmsrgy'
1311             else
1312              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1313             endif
1314           endif
1315 #if defined(AIX) || defined(PGI)
1316           open (ihist,file=histname,position='append')
1317 #else
1318           open (ihist,file=histname,access='append')
1319 #endif
1320           do i=0,nbin_rms
1321             do j=0,nbin_rgy
1322               write(ihist,'(2f8.2,$)')
1323      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1324               do ib=1,nT_h(iparm)
1325                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1326                   write(ihist,'(e14.5,$)')
1327      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1328      &              +potEmin
1329                 else
1330                   write(ihist,'(e14.5,$)') 1.0d6
1331                 endif
1332               enddo
1333               write (ihist,'(i2)') iparm
1334             enddo
1335           enddo
1336           close(ihist)
1337         endif
1338         endif
1339       enddo ! iparm
1340 #ifdef MPI
1341       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1342      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1343       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1344      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1345       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1346      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1347       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1348      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1349       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1350      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1351       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1352      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1353      &   WHAM_COMM,IERROR)
1354       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1355      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1356      &   WHAM_COMM,IERROR)
1357       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1358      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1359      &   WHAM_COMM,IERROR)
1360       if (me.eq.master) then
1361 #endif
1362       write (iout,'(/a)') 'Thermal characteristics of folding'
1363       if (nslice.eq.1) then
1364         nazwa=prefix
1365       else
1366         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1367       endif
1368       iln=ilen(nazwa)
1369       if (nparmset.eq.1 .and. .not.separate_parset) then
1370         nazwa=nazwa(:iln)//".thermal"
1371       else if (nparmset.eq.1 .and. separate_parset) then
1372         write(licz3,"(bz,i3.3)") myparm
1373         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1374       endif
1375       do iparm=1,nParmSet
1376       if (nparmset.gt.1) then
1377         write(licz3,"(bz,i3.3)") iparm
1378         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1379       endif
1380       open(34,file=nazwa)
1381       if (separate_parset) then
1382         write (iout,'(a,i3)') "Parameter set",myparm
1383       else
1384         write (iout,'(a,i3)') "Parameter set",iparm
1385       endif
1386       do i=0,NGridT
1387         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1388         if (betaT.ge.beta_h(1,iparm)) then
1389           potEmin=potEmin_all(1,iparm)
1390         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1391           potEmin=potEmin_all(nT_h(iparm),iparm)
1392         else
1393           do l=1,nT_h(iparm)-1
1394             if (betaT.le.beta_h(l,iparm) .and.
1395      &          betaT.gt.beta_h(l+1,iparm)) then
1396               potEmin=potEmin_all(l,iparm)
1397               exit
1398             endif
1399           enddo
1400         endif
1401         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1402         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1403      &    sumW(i,iparm)
1404         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1405      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1406         do j=1,nQ+2
1407           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1408           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1409      &     -sumQ(j,i,iparm)**2
1410           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1411      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1412         enddo
1413         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1414      &   (startGridT+i*delta_T))+potEmin
1415         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1416      &   sumW(i,iparm),sumE(i,iparm)
1417         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1418         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1419      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1420         write (iout,*)
1421         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1422      &   sumW(i,iparm),sumE(i,iparm)
1423         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1424         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1425      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1426         write (34,*)
1427       enddo
1428       close(34)
1429       enddo
1430       if (histout) then
1431       do t=0,tmax
1432         if (hfin_ent(t).gt.0.0d0) then
1433           liczba=t
1434           jj = mod(liczba,nbin1)
1435           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1436      &     hfin_ent(t)
1437           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1438      &      dmin+(jj+0.5d0)*delta,
1439      &     hfin_ent(t)
1440         endif
1441       enddo
1442       if (histfile) close(ihist)
1443       endif
1444
1445 #ifdef ZSCORE
1446 ! Write data for zscore
1447       if (nslice.eq.1) then
1448         zscname=prefix(:ilen(prefix))//".zsc"
1449       else
1450         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1451       endif
1452 #if defined(AIX) || defined(PGI)
1453       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1454 #else
1455       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1456 #endif
1457       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1458       do iparm=1,nParmSet
1459         write (izsc,'("NT=",i1)') nT_h(iparm)
1460       do ib=1,nT_h(iparm)
1461         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1462      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1463         jj = min0(nR(ib,iparm),7)
1464         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1465         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1466         write (izsc,'("&")')
1467         if (nR(ib,iparm).gt.7) then
1468           do ii=8,nR(ib,iparm),9
1469             jj = min0(nR(ib,iparm),ii+8)
1470             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1471             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1472             write (izsc,'("&")')
1473           enddo
1474         endif
1475         write (izsc,'("FI=",$)')
1476         jj=min0(nR(ib,iparm),7)
1477         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1478         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1479         write (izsc,'("&")')
1480         if (nR(ib,iparm).gt.7) then
1481           do ii=8,nR(ib,iparm),9
1482             jj = min0(nR(ib,iparm),ii+8)
1483             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1484             if (jj.eq.nR(ib,iparm)) then
1485               write (izsc,*)
1486             else
1487               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1488               write (izsc,'(t80,"&")')
1489             endif
1490           enddo
1491         endif
1492         do i=1,nR(ib,iparm)
1493           write (izsc,'("KH=",$)')
1494           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1495           write (izsc,'(" Q0=",$)')
1496           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1497           write (izsc,*)
1498         enddo
1499       enddo
1500       enddo
1501       close(izsc)
1502 #endif
1503 #ifdef MPI
1504       endif
1505 #endif
1506
1507       return
1508
1509       end