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