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,'(6f15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 535             write (iout,'(6f15.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,1.0d0/(1.987d-3*beta_h(j,i)),
 736      &      PotEmin_all(j,i)
 737         enddo
 738       enddo
 739       write (iout,*) "potEmin_min",potEmin_min
 740 #endif
 741
 742 #ifdef MPI
 743       do t=0,tmax
 744         hfin_ent_p(t)=0.0d0
 745       enddo
 746 #else
 747       do t=0,tmax
 748         hfin_ent(t)=0.0d0
 749       enddo
 750 #endif
 751       write (iout,*) "--------------hist"
 752 #ifdef MPI
 753       do iparm=1,nParmSet
 754         do i=0,nGridT
 755           sumW_p(i,iparm)=0.0d0
 756           sumE_p(i,iparm)=0.0d0
 757           sumEbis_p(i,iparm)=0.0d0
 758           sumEsq_p(i,iparm)=0.0d0
 759           do j=1,nQ+2
 760             sumQ_p(j,i,iparm)=0.0d0
 761             sumQsq_p(j,i,iparm)=0.0d0
 762             sumEQ_p(j,i,iparm)=0.0d0
 763           enddo
 764         enddo
 765       enddo
 766       upindE_p=0
 767 #else
 768       do iparm=1,nParmSet
 769         do i=0,nGridT
 770           sumW(i,iparm)=0.0d0
 771           sumE(i,iparm)=0.0d0
 772           sumEbis(i,iparm)=0.0d0
 773           sumEsq(i,iparm)=0.0d0
 774           do j=1,nQ+2
 775             sumQ(j,i,iparm)=0.0d0
 776             sumQsq(j,i,iparm)=0.0d0
 777             sumEQ(j,i,iparm)=0.0d0
 778           enddo
 779         enddo
 780       enddo
 781       upindE=0
 782 #endif
 783 c 8/26/05 entropy distribution
 784 #ifdef MPI
 785       entmin_p=1.0d10
 786       entmax_p=-1.0d10
 787       do t=1,scount(me1)
 788 c        ent=-dlog(entfac(t))
 789         ent=entfac(t)
 790         if (ent.lt.entmin_p) entmin_p=ent
 791         if (ent.gt.entmax_p) entmax_p=ent
 792       enddo
 793       write (iout,*) "entmin",entmin_p," entmax",entmax_p
 794       call flush(iout)
 795       call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,
 796      &  WHAM_COMM,IERROR)
 797       call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,
 798      &  WHAM_COMM,IERROR)
 799       ientmax=entmax-entmin
 800       if (ientmax.gt.2000) ientmax=2000
 801       write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
 802       call flush(iout)
 803       do t=1,scount(me1)
 804 c        ient=-dlog(entfac(t))-entmin
 805         ient=entfac(t)-entmin
 806         if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
 807       enddo
 808       call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,
 809      &  MPI_SUM,WHAM_COMM,IERROR)
 810       if (me1.eq.Master) then
 811         write (iout,*) "Entropy histogram"
 812         do i=0,ientmax
 813           write(iout,'(f15.4,i10)') entmin+i,histent(i)
 814         enddo
 815       endif
 816 #else
 817       entmin=1.0d10
 818       entmax=-1.0d10
 819       do t=1,ntot(islice)
 820         ent=entfac(t)
 821         if (ent.lt.entmin) entmin=ent
 822         if (ent.gt.entmax) entmax=ent
 823       enddo
 824       ientmax=-dlog(entmax)-entmin
 825       if (ientmax.gt.2000) ientmax=2000
 826       do t=1,ntot(islice)
 827         ient=entfac(t)-entmin
 828         if (ient.le.2000) histent(ient)=histent(ient)+1
 829       enddo
 830       write (iout,*) "Entropy histogram"
 831       do i=0,ientmax
 832         write(iout,'(2f15.4)') entmin+i,histent(i)
 833       enddo
 834 #endif
 835
 836 #ifdef MPI
 837 c      write (iout,*) "me1",me1," scount",scount(me1)
 838
 839       do iparm=1,nParmSet
 840
 841 #ifdef MPI
 842         do ib=1,nT_h(iparm)
 843           do t=0,tmax
 844             hfin_p(t,ib)=0.0d0
 845           enddo
 846         enddo
 847         do i=1,maxindE
 848           histE_p(i)=0.0d0
 849         enddo
 850 #else
 851         do ib=1,nT_h(iparm)
 852           do t=0,tmax
 853             hfin(t,ib)=0.0d0
 854           enddo
 855         enddo
 856         do i=1,maxindE
 857           histE(i)=0.0d0
 858         enddo
 859 #endif
 860         do ib=1,nT_h(iparm)
 861           do i=0,MaxBinRms
 862             do j=0,MaxBinRgy
 863               hrmsrgy(j,i,ib)=0.0d0
 864 #ifdef MPI
 865               hrmsrgy_p(j,i,ib)=0.0d0
 866 #endif
 867             enddo
 868           enddo
 869         enddo
 870
 871         do t=1,scount(me1)
 872 #else
 873         do t=1,ntot(islice)
 874 #endif
 875           ind = ind_point(t)
 876 #ifdef MPI
 877           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 878 #else
 879           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 880 #endif
 881           call restore_parm(iparm)
 882           evdw=enetb(21,t,iparm)
 883           evdw_t=enetb(1,t,iparm)
 884 #ifdef SCP14
 885           evdw2_14=enetb(17,t,iparm)
 886           evdw2=enetb(2,t,iparm)+evdw2_14
 887 #else
 888           evdw2=enetb(2,t,iparm)
 889           evdw2_14=0.0d0
 890 #endif
 891 #ifdef SPLITELE
 892           ees=enetb(3,t,iparm)
 893           evdw1=enetb(16,t,iparm)
 894 #else
 895           ees=enetb(3,t,iparm)
 896           evdw1=0.0d0
 897 #endif
 898           ecorr=enetb(4,t,iparm)
 899           ecorr5=enetb(5,t,iparm)
 900           ecorr6=enetb(6,t,iparm)
 901           eel_loc=enetb(7,t,iparm)
 902           eello_turn3=enetb(8,t,iparm)
 903           eello_turn4=enetb(9,t,iparm)
 904           eturn6=enetb(10,t,iparm)
 905           ebe=enetb(11,t,iparm)
 906           escloc=enetb(12,t,iparm)
 907           etors=enetb(13,t,iparm)
 908           etors_d=enetb(14,t,iparm)
 909           ehpb=enetb(15,t,iparm)
 910           estr=enetb(18,t,iparm)
 911           esccor=enetb(19,t,iparm)
 912           edihcnstr=enetb(20,t,iparm)
 913           do k=0,nGridT
 914             betaT=startGridT+k*delta_T
 915             temper=betaT
 916 c            fT=T0/betaT
 917 c            ft=2*T0/(T0+betaT)
 918             if (rescale_mode.eq.1) then
 919               quot=betaT/T0
 920               quotl=1.0d0
 921               kfacl=1.0d0
 922               do l=1,5
 923                 quotl1=quotl
 924                 quotl=quotl*quot
 925                 kfacl=kfacl*kfac
 926                 denom=kfacl-1.0d0+quotl
 927                 fT(l)=kfacl/denom
 928                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
 929                 ftbis(l)=l*kfacl*quotl1*
 930      &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
 931               enddo
 932 #if defined(FUNCTH)
 933               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 934      &                  320.0d0
 935               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 936              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 937      &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 938 #elif defined(FUNCT)
 939               fT(6)=betaT/T0
 940               ftprim(6)=1.0d0/T0
 941               ftbis(6)=0.0d0
 942 #else
 943               fT(6)=1.0d0
 944               ftprim(6)=0.0d0
 945               ftbis(6)=0.0d0
 946 #endif
 947             else if (rescale_mode.eq.2) then
 948               quot=betaT/T0
 949               quotl=1.0d0
 950               do l=1,5
 951                 quotl1=quotl
 952                 quotl=quotl*quot
 953                 eplus=dexp(quotl)
 954                 eminus=dexp(-quotl)
 955                 logfac=1.0d0/dlog(eplus+eminus)
 956                 tanhT=(eplus-eminus)/(eplus+eminus)
 957                 fT(l)=1.12692801104297249644d0*logfac
 958                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
 959                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
 960      &          2*l*quotl1/T0*logfac*
 961      &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
 962      &          +ftprim(l)*tanhT)
 963               enddo
 964 #if defined(FUNCTH)
 965               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
 966      &                 320.0d0
 967               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
 968              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
 969      &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
 970 #elif defined(FUNCT)
 971               fT(6)=betaT/T0
 972               ftprim(6)=1.0d0/T0
 973               ftbis(6)=0.0d0
 974 #else
 975               fT(6)=1.0d0
 976               ftprim(6)=0.0d0
 977               ftbis(6)=0.0d0
 978 #endif
 979             else if (rescale_mode.eq.0) then
 980               do l=1,5
 981                 fT(l)=1.0d0
 982                 ftprim(l)=0.0d0
 983               enddo
 984             else
 985               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
 986      &          rescale_mode
 987               call flush(iout)
 988               return1
 989             endif
 990 c            write (iout,*) "ftprim",ftprim
 991 c            write (iout,*) "ftbis",ftbis
 992             betaT=1.0d0/(1.987D-3*betaT)
 993             if (betaT.ge.beta_h(1,iparm)) then
 994               potEmin=potEmin_all(1,iparm)+
 995      &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
 996      &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
 997      &         (1.0/betaT-1.0/beta_h(1,iparm))
 998 #ifdef DEBUG
 999               write(iout,*) "first",temper,potEmin
1000 #endif
1001             else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
1002               potEmin=potEmin_all(nT_h(iparm),iparm)+
1003      &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
1004      &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
1005      &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
1006 #ifdef DEBUG
1007               write (iout,*) "last",temper,potEmin
1008 #endif
1009             else
1010               do l=1,nT_h(iparm)-1
1011                 if (betaT.le.beta_h(l,iparm) .and.
1012      &              betaT.gt.beta_h(l+1,iparm)) then
1013                   potEmin=potEmin_all(l,iparm)
1014 #ifdef DEBUG
1015                   write (iout,*) "l",l,
1016      &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1017      &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1018 #endif
1019                   exit
1020                 endif
1021               enddo
1022             endif
1023 #ifdef DEBUG
1024             write (iout,*) "k",k," potEmin",potEmin
1025 #endif
1026 #ifdef SPLITELE
1027             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1028      &      +wvdwpp*evdw1
1029      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1030      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1031      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1032      &      +ft(2)*wturn3*eello_turn3
1033      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1034      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1035      &      +wbond*estr
1036             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1037      &            +ftprim(1)*wtor*etors+
1038      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1039      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1040      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1041      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1042      &            ftprim(1)*wsccor*esccor
1043             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1044      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1045      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1046      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1047      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1048      &            ftbis(1)*wsccor*esccor
1049 #else
1050             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1051      &      +ft(1)*welec*(ees+evdw1)
1052      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1053      &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1054      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1055      &      +ft(2)*wturn3*eello_turn3
1056      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1057      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1058      &      +wbond*estr
1059             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1060      &           +ftprim(1)*wtor*etors+
1061      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1062      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1063      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1064      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1065      &            ftprim(1)*wsccor*esccor
1066             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1067      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1068      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1069      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1070      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1071      &            ftprim(1)*wsccor*esccor
1072 #endif
1073             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1074 #ifdef DEBUG
1075             write (iout,*) "iparm",iparm," t",t," temper",temper,
1076      &        " etot",etot," entfac",entfac(t),
1077      &        " efree",etot-entfac(t)/betaT," potEmin",potEmin,
1078      &        " boltz",-betaT*(etot-potEmin)+entfac(t),
1079      &        " weight",weight," ebis",ebis
1080 #endif
1081             etot=etot-temper*eprim
1082 #ifdef MPI
1083             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1084             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1085             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1086             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1087             do j=1,nQ+2
1088               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1089               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1090               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1091      &         +etot*q(j,t)*weight
1092             enddo
1093 #else
1094             sumW(k,iparm)=sumW(k,iparm)+weight
1095             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1096             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1097             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1098             do j=1,nQ+2
1099               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1100               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1101               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1102      &         +etot*q(j,t)*weight
1103             enddo
1104 #endif
1105           enddo
1106           indE = aint(potE(t,iparm)-aint(potEmin))
1107           if (indE.ge.0 .and. indE.le.maxinde) then
1108             if (indE.gt.upindE_p) upindE_p=indE
1109             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1110           endif
1111 #ifdef MPI
1112           do ib=1,nT_h(iparm)
1113             potEmin=potEmin_all(ib,iparm)
1114             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1115             hfin_p(ind,ib)=hfin_p(ind,ib)+
1116      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1117              if (rmsrgymap) then
1118                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1119                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1120                hrmsrgy_p(indrgy,indrms,ib)=
1121      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1122              endif
1123           enddo
1124 #else
1125           do ib=1,nT_h(iparm)
1126             potEmin=potEmin_all(ib,iparm)
1127             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1128             hfin(ind,ib)=hfin(ind,ib)+
1129      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1130              if (rmsrgymap) then
1131                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1132                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1133                hrmsrgy(indrgy,indrms,ib)=
1134      &           hrmsrgy(indrgy,indrms,ib)+expfac
1135              endif
1136           enddo
1137 #endif
1138         enddo ! t
1139         do ib=1,nT_h(iparm)
1140           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1141      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1142           if (rmsrgymap) then
1143           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1144      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1145      &       WHAM_COMM,IERROR)
1146           endif
1147         enddo
1148         call MPI_Reduce(upindE_p,upindE,1,
1149      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1150         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1151      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1152
1153         if (me1.eq.master) then
1154
1155         if (histout) then
1156
1157         write (iout,'(6x,$)')
1158         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1159      &   ib=1,nT_h(iparm))
1160         write (iout,*)
1161
1162         write (iout,'(/a)') 'Final histograms'
1163         if (histfile) then
1164           if (nslice.eq.1) then
1165             if (separate_parset) then
1166               write(licz3,"(bz,i3.3)") myparm
1167               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1168             else
1169               histname=prefix(:ilen(prefix))//'.hist'
1170             endif
1171           else
1172             if (separate_parset) then
1173               write(licz3,"(bz,i3.3)") myparm
1174               histname=prefix(:ilen(prefix))//'_par'//licz3//
1175      &         '_slice_'//licz2//'.hist'
1176             else
1177               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1178             endif
1179           endif
1180 #if defined(AIX) || defined(PGI)
1181           open (ihist,file=histname,position='append')
1182 #else
1183           open (ihist,file=histname,access='append')
1184 #endif
1185         endif
1186
1187         do t=0,tmax
1188           liczba=t
1189           sumH=0.0d0
1190           do ib=1,nT_h(iparm)
1191             sumH=sumH+hfin(t,ib)
1192           enddo
1193           if (sumH.gt.0.0d0) then
1194             do j=1,nQ
1195               jj = mod(liczba,nbin1)
1196               liczba=liczba/nbin1
1197               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1198               if (histfile)
1199      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1200             enddo
1201             do ib=1,nT_h(iparm)
1202               write (iout,'(e20.10,$)') hfin(t,ib)
1203               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1204             enddo
1205             write (iout,'(i5)') iparm
1206             if (histfile) write (ihist,'(i5)') iparm
1207           endif
1208         enddo
1209
1210         endif
1211
1212         if (entfile) then
1213           if (nslice.eq.1) then
1214             if (separate_parset) then
1215               write(licz3,"(bz,i3.3)") myparm
1216               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1217             else
1218               histname=prefix(:ilen(prefix))//'.ent'
1219             endif
1220           else
1221             if (separate_parset) then
1222               write(licz3,"(bz,i3.3)") myparm
1223               histname=prefix(:ilen(prefix))//'par_'//licz3//
1224      &           '_slice_'//licz2//'.ent'
1225             else
1226               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1227             endif
1228           endif
1229 #if defined(AIX) || defined(PGI)
1230           open (ihist,file=histname,position='append')
1231 #else
1232           open (ihist,file=histname,access='append')
1233 #endif
1234           write (ihist,'(a)') "# Microcanonical entropy"
1235           do i=0,upindE
1236             write (ihist,'(f8.0,$)') dint(potEmin)+i
1237             if (histE(i).gt.0.0e0) then
1238               write (ihist,'(f15.5,$)') dlog(histE(i))
1239             else
1240               write (ihist,'(f15.5,$)') 0.0d0
1241             endif
1242           enddo
1243           write (ihist,*)
1244           close(ihist)
1245         endif
1246         write (iout,*) "Microcanonical entropy"
1247         do i=0,upindE
1248           write (iout,'(f8.0,$)') dint(potEmin)+i
1249           if (histE(i).gt.0.0e0) then
1250             write (iout,'(f15.5,$)') dlog(histE(i))
1251           else
1252             write (iout,'(f15.5,$)') 0.0d0
1253           endif
1254           write (iout,*)
1255         enddo
1256         if (rmsrgymap) then
1257           if (nslice.eq.1) then
1258             if (separate_parset) then
1259               write(licz3,"(bz,i3.3)") myparm
1260               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1261             else
1262               histname=prefix(:ilen(prefix))//'.rmsrgy'
1263             endif
1264           else
1265             if (separate_parset) then
1266               write(licz3,"(bz,i3.3)") myparm
1267               histname=prefix(:ilen(prefix))//'_par'//licz3//
1268      &         '_slice_'//licz2//'.rmsrgy'
1269             else
1270              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1271             endif
1272           endif
1273 #if defined(AIX) || defined(PGI)
1274           open (ihist,file=histname,position='append')
1275 #else
1276           open (ihist,file=histname,access='append')
1277 #endif
1278           do i=0,nbin_rms
1279             do j=0,nbin_rgy
1280               write(ihist,'(2f8.2,$)')
1281      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1282               do ib=1,nT_h(iparm)
1283                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1284                   write(ihist,'(e14.5,$)')
1285      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1286      &              +potEmin
1287                 else
1288                   write(ihist,'(e14.5,$)') 1.0d6
1289                 endif
1290               enddo
1291               write (ihist,'(i2)') iparm
1292             enddo
1293           enddo
1294           close(ihist)
1295         endif
1296         endif
1297       enddo ! iparm
1298 #ifdef MPI
1299       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1300      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1301       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1302      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1303       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1304      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1305       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1306      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1307       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1308      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1309       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1310      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1311      &   WHAM_COMM,IERROR)
1312       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1313      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1314      &   WHAM_COMM,IERROR)
1315       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1316      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1317      &   WHAM_COMM,IERROR)
1318       if (me.eq.master) then
1319 #endif
1320       write (iout,'(/a)') 'Thermal characteristics of folding'
1321       if (nslice.eq.1) then
1322         nazwa=prefix
1323       else
1324         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1325       endif
1326       iln=ilen(nazwa)
1327       if (nparmset.eq.1 .and. .not.separate_parset) then
1328         nazwa=nazwa(:iln)//".thermal"
1329       else if (nparmset.eq.1 .and. separate_parset) then
1330         write(licz3,"(bz,i3.3)") myparm
1331         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1332       endif
1333       do iparm=1,nParmSet
1334       if (nparmset.gt.1) then
1335         write(licz3,"(bz,i3.3)") iparm
1336         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1337       endif
1338       open(34,file=nazwa)
1339       if (separate_parset) then
1340         write (iout,'(a,i3)') "Parameter set",myparm
1341       else
1342         write (iout,'(a,i3)') "Parameter set",iparm
1343       endif
1344       do i=0,NGridT
1345         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1346         if (betaT.ge.beta_h(1,iparm)) then
1347           potEmin=potEmin_all(1,iparm)
1348         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1349           potEmin=potEmin_all(nT_h(iparm),iparm)
1350         else
1351           do l=1,nT_h(iparm)-1
1352             if (betaT.le.beta_h(l,iparm) .and.
1353      &          betaT.gt.beta_h(l+1,iparm)) then
1354               potEmin=potEmin_all(l,iparm)
1355               exit
1356             endif
1357           enddo
1358         endif
1359
1360 c        write (iout,*) "i",i," potEmin",potEmin
1361
1362         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1363         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1364      &    sumW(i,iparm)
1365         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1366      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1367         do j=1,nQ+2
1368           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1369           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1370      &     -sumQ(j,i,iparm)**2
1371           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1372      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1373         enddo
1374         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1375      &   (startGridT+i*delta_T))+potEmin
1376         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1377      &   sumW(i,iparm),sumE(i,iparm)
1378         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1379         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1380      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1381         write (iout,*)
1382         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1383      &   sumW(i,iparm),sumE(i,iparm)
1384         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1385         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1386      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1387         write (34,*)
1388       enddo
1389       close(34)
1390       enddo
1391       if (histout) then
1392       do t=0,tmax
1393         if (hfin_ent(t).gt.0.0d0) then
1394           liczba=t
1395           jj = mod(liczba,nbin1)
1396           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1397      &     hfin_ent(t)
1398           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1399      &      dmin+(jj+0.5d0)*delta,
1400      &     hfin_ent(t)
1401         endif
1402       enddo
1403       if (histfile) close(ihist)
1404       endif
1405
1406 #ifdef ZSCORE
1407 ! Write data for zscore
1408       if (nslice.eq.1) then
1409         zscname=prefix(:ilen(prefix))//".zsc"
1410       else
1411         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1412       endif
1413 #if defined(AIX) || defined(PGI)
1414       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1415 #else
1416       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1417 #endif
1418       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1419       do iparm=1,nParmSet
1420         write (izsc,'("NT=",i1)') nT_h(iparm)
1421       do ib=1,nT_h(iparm)
1422         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1423      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1424         jj = min0(nR(ib,iparm),7)
1425         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1426         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1427         write (izsc,'("&")')
1428         if (nR(ib,iparm).gt.7) then
1429           do ii=8,nR(ib,iparm),9
1430             jj = min0(nR(ib,iparm),ii+8)
1431             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1432             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1433             write (izsc,'("&")')
1434           enddo
1435         endif
1436         write (izsc,'("FI=",$)')
1437         jj=min0(nR(ib,iparm),7)
1438         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1439         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1440         write (izsc,'("&")')
1441         if (nR(ib,iparm).gt.7) then
1442           do ii=8,nR(ib,iparm),9
1443             jj = min0(nR(ib,iparm),ii+8)
1444             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1445             if (jj.eq.nR(ib,iparm)) then
1446               write (izsc,*)
1447             else
1448               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1449               write (izsc,'(t80,"&")')
1450             endif
1451           enddo
1452         endif
1453         do i=1,nR(ib,iparm)
1454           write (izsc,'("KH=",$)')
1455           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1456           write (izsc,'(" Q0=",$)')
1457           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1458           write (izsc,*)
1459         enddo
1460       enddo
1461       enddo
1462       close(izsc)
1463 #endif
1464 #ifdef MPI
1465       endif
1466 #endif
1467
1468       return
1469
1470       end