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