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