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