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