source/wham/src-NEWSC-NEWCORR/wham_calc1.F

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