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