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
  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+nss*ebr+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+nss*ebr+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,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 536             write (iout,'(8f10.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+nss*ebr+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+nss*ebr+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,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 c              write(iout,*) "first",temper,potEmin
 996             else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
 997               potEmin=potEmin_all(nT_h(iparm),iparm)
 998 c              write (iout,*) "last",temper,potEmin
 999             else
1000               do l=1,nT_h(iparm)-1
1001                 if (betaT.le.beta_h(l,iparm) .and.
1002      &              betaT.gt.beta_h(l+1,iparm)) then
1003                   potEmin=potEmin_all(l,iparm)
1004 c                  write (iout,*) "l",l,
1005 c     &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
1006 c     &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
1007                   exit
1008                 endif
1009               enddo
1010             endif
1011 c            write (iout,*) ib," PotEmin",potEmin
1012 #ifdef SPLITELE
1013             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
1014      &      +wvdwpp*evdw1
1015      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1016      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1017      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1018      &      +ft(2)*wturn3*eello_turn3
1019      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
1020      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1021      &      +wbond*estr
1022             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
1023      &            +ftprim(1)*wtor*etors+
1024      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1025      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1026      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1027      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1028      &            ftprim(1)*wsccor*esccor
1029             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
1030      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1031      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1032      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1033      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1034      &            ftbis(1)*wsccor*esccor
1035 #else
1036             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
1037      &      +ft(1)*welec*(ees+evdw1)
1038      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
1039      &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1040      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
1041      &      +ft(2)*wturn3*eello_turn3
1042      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
1043      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
1044      &      +wbond*estr
1045             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
1046      &           +ftprim(1)*wtor*etors+
1047      &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1048      &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1049      &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
1050      &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1051      &            ftprim(1)*wsccor*esccor
1052             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
1053      &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1054      &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1055      &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
1056      &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1057      &            ftprim(1)*wsccor*esccor
1058 #endif
1059             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1060 #ifdef DEBUG
1061             write (iout,*) "iparm",iparm," t",t," temper",temper,
1062      &        " etot",etot," entfac",entfac(t),
1063      &        " efree",etot-entfac(t)/betaT," potEmin",potEmin,
1064      &        " boltz",-betaT*(etot-potEmin)+entfac(t),
1065      &        " weight",weight," ebis",ebis
1066 #endif
1067             etot=etot-temper*eprim
1068 #ifdef MPI
1069             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1070             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1071             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1072             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1073             do j=1,nQ+2
1074               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1075               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1076               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm)
1077      &         +etot*q(j,t)*weight
1078             enddo
1079 #else
1080             sumW(k,iparm)=sumW(k,iparm)+weight
1081             sumE(k,iparm)=sumE(k,iparm)+etot*weight
1082             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1083             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1084             do j=1,nQ+2
1085               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1086               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1087               sumEQ(j,k,iparm)=sumEQ(j,k,iparm)
1088      &         +etot*q(j,t)*weight
1089             enddo
1090 #endif
1091           enddo
1092           indE = aint(potE(t,iparm)-aint(potEmin))
1093           if (indE.ge.0 .and. indE.le.maxinde) then
1094             if (indE.gt.upindE_p) upindE_p=indE
1095             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1096           endif
1097 #ifdef MPI
1098           do ib=1,nT_h(iparm)
1099             potEmin=potEmin_all(ib,iparm)
1100             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1101             hfin_p(ind,ib)=hfin_p(ind,ib)+
1102      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1103              if (rmsrgymap) then
1104                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1105                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1106                hrmsrgy_p(indrgy,indrms,ib)=
1107      &           hrmsrgy_p(indrgy,indrms,ib)+expfac
1108              endif
1109           enddo
1110 #else
1111           do ib=1,nT_h(iparm)
1112             potEmin=potEmin_all(ib,iparm)
1113             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1114             hfin(ind,ib)=hfin(ind,ib)+
1115      &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1116              if (rmsrgymap) then
1117                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1118                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1119                hrmsrgy(indrgy,indrms,ib)=
1120      &           hrmsrgy(indrgy,indrms,ib)+expfac
1121              endif
1122           enddo
1123 #endif
1124         enddo ! t
1125         do ib=1,nT_h(iparm)
1126           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
1127      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1128           if (rmsrgymap) then
1129           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),
1130      &   (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1131      &       WHAM_COMM,IERROR)
1132           endif
1133         enddo
1134         call MPI_Reduce(upindE_p,upindE,1,
1135      &    MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
1136         call MPI_Reduce(histE_p(0),histE(0),maxindE,
1137      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1138
1139         if (me1.eq.master) then
1140
1141         if (histout) then
1142
1143         write (iout,'(6x,$)')
1144         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),
1145      &   ib=1,nT_h(iparm))
1146         write (iout,*)
1147
1148         write (iout,'(/a)') 'Final histograms'
1149         if (histfile) then
1150           if (nslice.eq.1) then
1151             if (separate_parset) then
1152               write(licz3,"(bz,i3.3)") myparm
1153               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
1154             else
1155               histname=prefix(:ilen(prefix))//'.hist'
1156             endif
1157           else
1158             if (separate_parset) then
1159               write(licz3,"(bz,i3.3)") myparm
1160               histname=prefix(:ilen(prefix))//'_par'//licz3//
1161      &         '_slice_'//licz2//'.hist'
1162             else
1163               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
1164             endif
1165           endif
1166 #if defined(AIX) || defined(PGI)
1167           open (ihist,file=histname,position='append')
1168 #else
1169           open (ihist,file=histname,access='append')
1170 #endif
1171         endif
1172
1173         do t=0,tmax
1174           liczba=t
1175           sumH=0.0d0
1176           do ib=1,nT_h(iparm)
1177             sumH=sumH+hfin(t,ib)
1178           enddo
1179           if (sumH.gt.0.0d0) then
1180             do j=1,nQ
1181               jj = mod(liczba,nbin1)
1182               liczba=liczba/nbin1
1183               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1184               if (histfile)
1185      &           write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
1186             enddo
1187             do ib=1,nT_h(iparm)
1188               write (iout,'(e20.10,$)') hfin(t,ib)
1189               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
1190             enddo
1191             write (iout,'(i5)') iparm
1192             if (histfile) write (ihist,'(i5)') iparm
1193           endif
1194         enddo
1195
1196         endif
1197
1198         if (entfile) then
1199           if (nslice.eq.1) then
1200             if (separate_parset) then
1201               write(licz3,"(bz,i3.3)") myparm
1202               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
1203             else
1204               histname=prefix(:ilen(prefix))//'.ent'
1205             endif
1206           else
1207             if (separate_parset) then
1208               write(licz3,"(bz,i3.3)") myparm
1209               histname=prefix(:ilen(prefix))//'par_'//licz3//
1210      &           '_slice_'//licz2//'.ent'
1211             else
1212               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
1213             endif
1214           endif
1215 #if defined(AIX) || defined(PGI)
1216           open (ihist,file=histname,position='append')
1217 #else
1218           open (ihist,file=histname,access='append')
1219 #endif
1220           write (ihist,'(a)') "# Microcanonical entropy"
1221           do i=0,upindE
1222             write (ihist,'(f8.0,$)') dint(potEmin)+i
1223             if (histE(i).gt.0.0e0) then
1224               write (ihist,'(f15.5,$)') dlog(histE(i))
1225             else
1226               write (ihist,'(f15.5,$)') 0.0d0
1227             endif
1228           enddo
1229           write (ihist,*)
1230           close(ihist)
1231         endif
1232         write (iout,*) "Microcanonical entropy"
1233         do i=0,upindE
1234           write (iout,'(f8.0,$)') dint(potEmin)+i
1235           if (histE(i).gt.0.0e0) then
1236             write (iout,'(f15.5,$)') dlog(histE(i))
1237           else
1238             write (iout,'(f15.5,$)') 0.0d0
1239           endif
1240           write (iout,*)
1241         enddo
1242         if (rmsrgymap) then
1243           if (nslice.eq.1) then
1244             if (separate_parset) then
1245               write(licz3,"(bz,i3.3)") myparm
1246               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
1247             else
1248               histname=prefix(:ilen(prefix))//'.rmsrgy'
1249             endif
1250           else
1251             if (separate_parset) then
1252               write(licz3,"(bz,i3.3)") myparm
1253               histname=prefix(:ilen(prefix))//'_par'//licz3//
1254      &         '_slice_'//licz2//'.rmsrgy'
1255             else
1256              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
1257             endif
1258           endif
1259 #if defined(AIX) || defined(PGI)
1260           open (ihist,file=histname,position='append')
1261 #else
1262           open (ihist,file=histname,access='append')
1263 #endif
1264           do i=0,nbin_rms
1265             do j=0,nbin_rgy
1266               write(ihist,'(2f8.2,$)')
1267      &          rgymin+deltrgy*j,rmsmin+deltrms*i
1268               do ib=1,nT_h(iparm)
1269                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
1270                   write(ihist,'(e14.5,$)')
1271      &              -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm)
1272      &              +potEmin
1273                 else
1274                   write(ihist,'(e14.5,$)') 1.0d6
1275                 endif
1276               enddo
1277               write (ihist,'(i2)') iparm
1278             enddo
1279           enddo
1280           close(ihist)
1281         endif
1282         endif
1283       enddo ! iparm
1284 #ifdef MPI
1285       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,
1286      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1287       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,
1288      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1289       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,
1290      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1291       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,
1292      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1293       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,
1294      &   MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
1295       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),
1296      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1297      &   WHAM_COMM,IERROR)
1298       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),
1299      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1300      &   WHAM_COMM,IERROR)
1301       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),
1302      &   MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,
1303      &   WHAM_COMM,IERROR)
1304       if (me.eq.master) then
1305 #endif
1306       write (iout,'(/a)') 'Thermal characteristics of folding'
1307       if (nslice.eq.1) then
1308         nazwa=prefix
1309       else
1310         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
1311       endif
1312       iln=ilen(nazwa)
1313       if (nparmset.eq.1 .and. .not.separate_parset) then
1314         nazwa=nazwa(:iln)//".thermal"
1315       else if (nparmset.eq.1 .and. separate_parset) then
1316         write(licz3,"(bz,i3.3)") myparm
1317         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1318       endif
1319       do iparm=1,nParmSet
1320       if (nparmset.gt.1) then
1321         write(licz3,"(bz,i3.3)") iparm
1322         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
1323       endif
1324       open(34,file=nazwa)
1325       if (separate_parset) then
1326         write (iout,'(a,i3)') "Parameter set",myparm
1327       else
1328         write (iout,'(a,i3)') "Parameter set",iparm
1329       endif
1330       do i=0,NGridT
1331         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
1332         if (betaT.ge.beta_h(1,iparm)) then
1333           potEmin=potEmin_all(1,iparm)
1334         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
1335           potEmin=potEmin_all(nT_h(iparm),iparm)
1336         else
1337           do l=1,nT_h(iparm)-1
1338             if (betaT.le.beta_h(l,iparm) .and.
1339      &          betaT.gt.beta_h(l+1,iparm)) then
1340               potEmin=potEmin_all(l,iparm)
1341               exit
1342             endif
1343           enddo
1344         endif
1345         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
1346         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
1347      &    sumW(i,iparm)
1348         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm)
1349      &    -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
1350         do j=1,nQ+2
1351           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
1352           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm)
1353      &     -sumQ(j,i,iparm)**2
1354           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
1355      &     -sumQ(j,i,iparm)*sumE(i,iparm)
1356         enddo
1357         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
1358      &   (startGridT+i*delta_T))+potEmin
1359         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1360      &   sumW(i,iparm),sumE(i,iparm)
1361         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1362         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1363      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1364         write (iout,*)
1365         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
1366      &   sumW(i,iparm),sumE(i,iparm)
1367         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
1368         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),
1369      &   (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
1370         write (34,*)
1371       enddo
1372       close(34)
1373       enddo
1374       if (histout) then
1375       do t=0,tmax
1376         if (hfin_ent(t).gt.0.0d0) then
1377           liczba=t
1378           jj = mod(liczba,nbin1)
1379           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,
1380      &     hfin_ent(t)
1381           if (histfile) write (ihist,'(f6.3,e20.10," ent")')
1382      &      dmin+(jj+0.5d0)*delta,
1383      &     hfin_ent(t)
1384         endif
1385       enddo
1386       if (histfile) close(ihist)
1387       endif
1388
1389 #ifdef ZSCORE
1390 ! Write data for zscore
1391       if (nslice.eq.1) then
1392         zscname=prefix(:ilen(prefix))//".zsc"
1393       else
1394         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
1395       endif
1396 #if defined(AIX) || defined(PGI)
1397       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
1398 #else
1399       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
1400 #endif
1401       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
1402       do iparm=1,nParmSet
1403         write (izsc,'("NT=",i1)') nT_h(iparm)
1404       do ib=1,nT_h(iparm)
1405         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)')
1406      &    1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
1407         jj = min0(nR(ib,iparm),7)
1408         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
1409         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
1410         write (izsc,'("&")')
1411         if (nR(ib,iparm).gt.7) then
1412           do ii=8,nR(ib,iparm),9
1413             jj = min0(nR(ib,iparm),ii+8)
1414             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
1415             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
1416             write (izsc,'("&")')
1417           enddo
1418         endif
1419         write (izsc,'("FI=",$)')
1420         jj=min0(nR(ib,iparm),7)
1421         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
1422         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
1423         write (izsc,'("&")')
1424         if (nR(ib,iparm).gt.7) then
1425           do ii=8,nR(ib,iparm),9
1426             jj = min0(nR(ib,iparm),ii+8)
1427             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
1428             if (jj.eq.nR(ib,iparm)) then
1429               write (izsc,*)
1430             else
1431               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
1432               write (izsc,'(t80,"&")')
1433             endif
1434           enddo
1435         endif
1436         do i=1,nR(ib,iparm)
1437           write (izsc,'("KH=",$)')
1438           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
1439           write (izsc,'(" Q0=",$)')
1440           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
1441           write (izsc,*)
1442         enddo
1443       enddo
1444       enddo
1445       close(izsc)
1446 #endif
1447 #ifdef MPI
1448       endif
1449 #endif
1450
1451       return
1452
1453       end