source/wham/wham_calc.F90

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