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), &
  75                      fi_p_min(MaxR,MaxT_h,nParmSet) !(MaxR,MaxT_h,Max_Parm)
  76       real(kind=8),dimension(0:nGridT,nParmSet) :: sumW_p,sumE_p,&
  77               sumEbis_p,sumEsq_p !(0:nGridT,Max_Parm)
  78       real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ_p,&
  79               sumQsq_p,sumEQ_p,sumEprim_p !(MaxQ1,0:nGridT,Max_Parm)
  80       real(kind=8) :: hfin_p(0:MaxHdim,maxT_h),&
  81               hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,&
  82               hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h),weimax_(0:ngridT)
  83       real(kind=8) :: rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
  84       real(kind=8) :: potEmin_t,potEmin_t_all(maxT_h,Max_Parm)!,entmin_p,entmax_p
  85 !      integer :: histent_p(0:2000)
  86       logical :: lprint=.true.
  87 #endif
  88       real(kind=8) :: delta_T=1.0d0,iientmax
  89       real(kind=8) :: rgymin,rmsmin,rgymax,rmsmax
  90       real(kind=8),dimension(0:nGridT,nParmSet) :: sumW,sumE,&
  91               sumEsq,sumEprim,sumEbis !(0:NGridT,Max_Parm)
  92       real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ,&
  93               sumQsq,sumEQ !(MaxQ1,0:NGridT,Max_Parm)
  94       real(kind=8) :: betaT,weight,econstr
  95       real(kind=8) :: fi(MaxR,MaxT_h,nParmSet),& !(MaxR,maxT_h,Max_Parm)
  96               ddW,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,&
  97               hfin(0:MaxHdim,maxT_h),histE(0:maxindE),&
  98               hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),&
  99               potEmin,ent,&
 100               hfin_ent(0:MaxHdim),vmax,aux,fi_min(MaxR,maxT_h,nParmSet), &
 101               potEmin_all(maxT_h,Max_Parm),potEmin_min,entfac_min,&
 102               weimax(0:nGridT,Max_Parm)
 103       real(kind=8) :: fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,&
 104         eprim,ebis,temper,kfac=2.4d0,T0=300.0d0,startGridT=200.0d0,&
 105         eplus,eminus,logfac,tanhT,tt
 106       real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
 107         escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
 108         eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, &
 109         ecationcation,ecation_prot, evdwpp,eespp ,evdwpsb,eelpsb, &
 110         evdwsb, eelsb, estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
 111         ecorr_nucl, ecorr3_nucl,epeppho, escpho, epepbase,escbase,ecation_nucl,&
 112         elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang,eliptran,&
 113         ehomology_constr
 114
 115
 116       integer :: ind_point(maxpoint),upindE,indE
 117       character(len=16) :: plik
 118       character(len=1) :: licz1
 119       character(len=2) :: licz2
 120       character(len=3) :: licz3
 121       character(len=128) :: nazwa
 122 !      integer ilen
 123 !      external ilen
 124 !el      ientmax=0
 125 !el      ent=0.0d0
 126       write(licz2,'(bz,i2.2)') islice
 127       nbin1 = 1.0d0/delta
 128       write (iout,'(//80(1h-)/"Solving WHAM equations for slice",&
 129         i2/80(1h-)//)') islice
 130       write (iout,*) "delta",delta," nbin1",nbin1
 131       write (iout,*) "MaxN",MaxN," MaxQ",MaxQ," MaHdim",MaxHdim
 132       call flush(iout)
 133       dmin=0.0d0
 134       tmax=0
 135       potEmin=1.0d10
 136       do i=1,nParmset
 137         do j=1,nT_h(i)
 138           potEmin_all(j,i)=1.0d15
 139         enddo
 140       enddo
 141       entfac_min=1.0d10
 142       rgymin=1.0d10
 143       rmsmin=1.0d10
 144       rgymax=0.0d0
 145       rmsmax=0.0d0
 146       do t=0,MaxN
 147         htot(t)=0
 148       enddo
 149 #ifdef MPI
 150       do i=1,scount(me1)
 151 #else
 152       do i=1,ntot(islice)
 153 #endif
 154         do j=1,nParmSet
 155           if (potE(i,j).le.potEmin) potEmin=potE(i,j)
 156         enddo
 157         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
 158         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
 159         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
 160         if (q(nQ+2,i).gt.rgymax) rgymax=q(nQ+2,i)
 161         ind_point(i)=0
 162         do j=nQ,1,-1
 163           ind=(q(j,i)-dmin+1.0d-8)/delta
 164           if (j.eq.1) then
 165             ind_point(i)=ind_point(i)+ind
 166           else
 167             ind_point(i)=ind_point(i)+nbin1**(j-1)*ind
 168           endif
 169 !          write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
 170 !     &      ind_point(i)
 171           call flush(iout)
 172           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
 173             write (iout,*) "Error - index exceeds range for point",i,&
 174             " q=",q(j,i)," ind",ind_point(i)
 175 #ifdef MPI
 176             write (iout,*) "Processor",me1
 177             call flush(iout)
 178             call MPI_Abort(MPI_COMM_WORLD, Ierror, Errcode )
 179 #endif
 180             stop
 181           endif
 182         enddo ! j
 183         if (ind_point(i).gt.tmax) tmax=ind_point(i)
 184         htot(ind_point(i))=htot(ind_point(i))+1
 185 #ifdef DEBUG
 186         write (iout,*) "i",i,"q",(q(j,i),j=1,nQ)," ind",ind_point(i),&
 187          " htot",htot(ind_point(i))
 188         call flush(iout)
 189 #endif
 190       enddo ! i
 191       call flush(iout)
 192
 193       nbin=nbin1**nQ-1
 194       write (iout,'(a)') "Numbers of counts in Q bins"
 195       do t=0,tmax
 196         if (htot(t).gt.0) then
 197         write (iout,'(i15,$)') t
 198         liczbaW=t
 199         do j=1,nQ
 200           jj = mod(liczbaW,nbin1)
 201           liczbaW=liczbaW/nbin1
 202           write (iout,'(i5,$)') jj
 203         enddo
 204         write (iout,'(i8)') htot(t)
 205         endif
 206       enddo
 207       do iparm=1,nParmSet
 208       write (iout,'(a,i3)') "Number of data points for parameter set",&
 209        iparm
 210       write (iout,'(i7,$)') ((snk(m,ib,iparm,islice),m=1,nr(ib,iparm)),&
 211        ib=1,nT_h(iparm))
 212       write (iout,'(i8)') stot(islice)
 213       write (iout,'(a)')
 214       enddo
 215       call flush(iout)
 216
 217 #ifdef MPI
 218       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,&
 219         WHAM_COMM,IERROR)
 220       tmax=tmax_t
 221 !      call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,&
 222 !        MPI_MIN,WHAM_COMM,IERROR) !????
 223       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,&
 224         MPI_MIN,WHAM_COMM,IERROR)
 225       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,&
 226         MPI_MAX,WHAM_COMM,IERROR)
 227       call MPI_AllReduce(rgymin,rgymin_t,1,MPI_DOUBLE_PRECISION,&
 228         MPI_MIN,WHAM_COMM,IERROR)
 229       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,&
 230         MPI_MAX,WHAM_COMM,IERROR)
 231 !      potEmin=potEmin_t !/2 try now??
 232       rgymin=rgymin_t
 233       rgymax=rgymax_t
 234       rmsmin=rmsmin_t
 235       rmsmax=rmsmax_t
 236       write (iout,*) "potEmin",potEmin
 237 #endif
 238       rmsmin=deltrms*dint(rmsmin/deltrms)
 239       rmsmax=deltrms*dint(rmsmax/deltrms)
 240       rgymin=deltrms*dint(rgymin/deltrgy)
 241       rgymax=deltrms*dint(rgymax/deltrgy)
 242       nbin_rms=(rmsmax-rmsmin)/deltrms
 243       nbin_rgy=(rgymax-rgymin)/deltrgy
 244       write (iout,*) "rmsmin",rmsmin," rmsmax",rmsmax," rgymin",rgymin,&
 245        " rgymax",rgymax," nbin_rms",nbin_rms," nbin_rgy",nbin_rgy
 246       nFi=0
 247       do i=1,nParmSet
 248         do j=1,nT_h(i)
 249           nFi=nFi+nR(j,i)
 250         enddo
 251       enddo
 252       write (iout,*) "nFi",nFi
 253 ! Compute the Boltzmann factor corresponing to restrain potentials in different
 254 ! simulations.
 255 #ifdef MPI
 256       do i=1,scount(me1)
 257 #else
 258       do i=1,ntot(islice)
 259 #endif
 260 !        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 261         do iparm=1,nParmSet
 262 #ifdef DEBUG
 263           write (iout,'(2i5,21f8.2)') i,iparm,&
 264            (enetb(k,i,iparm),k=1,21)
 265 #endif
 266           call restore_parm(iparm)
 267 #ifdef DEBUG
 268           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
 269             wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
 270             wtor_d,wsccor,wbond,wcatcat
 271 #endif
 272           do ib=1,nT_h(iparm)
 273 !el old rascale weights
 274 !
 275 !            if (rescale_modeW.eq.1) then
 276 !              quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 277 !              quotl=1.0d0
 278 !              kfacl=1.0d0
 279 !              do l=1,5
 280 !                quotl1=quotl
 281 !                quotl=quotl*quot
 282 !                kfacl=kfacl*kfac
 283 !                fT(l)=kfacl/(kfacl-1.0d0+quotl)
 284 !              enddo
 285 !#if defined(FUNCTH)
 286 !              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 287 !              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 288 !#elif defined(FUNCT)
 289 !              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 290 !#else
 291 !              ft(6)=1.0d0
 292 !#endif
 293 !            else if (rescale_modeW.eq.2) then
 294 !              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 295 !              quotl=1.0d0
 296 !              do l=1,5
 297 !                quotl=quotl*quot
 298 !                fT(l)=1.12692801104297249644d0/ &
 299 !                   dlog(dexp(quotl)+dexp(-quotl))
 300 !              enddo
 301 !#if defined(FUNCTH)
 302 !              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 303 !              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 304 !#elif defined(FUNCT)
 305 !              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 306 !#else
 307 !              ft(6)=1.0d0
 308 !#endif
 309 !              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 310 !            else if (rescale_modeW.eq.0) then
 311 !              do l=1,6
 312 !                fT(l)=1.0d0
 313 !              enddo
 314 !            else
 315 !              write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
 316 !                rescale_modeW
 317 !              call flush(iout)
 318 !              return 1
 319 !            endif
 320 ! el end old rescale weights
 321             call rescale_weights(1.0d0/(beta_h(ib,iparm)*1.987D-3))
 322
 323 !            call etot(enetb(0,i,iparm))
 324             evdw=enetb(1,i,iparm)
 325 !            evdw_t=enetb(21,i,iparm)
 326             evdw_t=enetb(20,i,iparm)
 327 #ifdef SCP14
 328 !            evdw2_14=enetb(17,i,iparm)
 329             evdw2_14=enetb(18,i,iparm)
 330             evdw2=enetb(2,i,iparm)+evdw2_14
 331 #else
 332             evdw2=enetb(2,i,iparm)
 333             evdw2_14=0.0d0
 334 #endif
 335 #ifdef SPLITELE
 336             ees=enetb(3,i,iparm)
 337             evdw1=enetb(16,i,iparm)
 338 #else
 339             ees=enetb(3,i,iparm)
 340             evdw1=0.0d0
 341 #endif
 342             ecorr=enetb(4,i,iparm)
 343             ecorr5=enetb(5,i,iparm)
 344             ecorr6=enetb(6,i,iparm)
 345             eel_loc=enetb(7,i,iparm)
 346             eello_turn3=enetb(8,i,iparm)
 347             eello_turn4=enetb(9,i,iparm)
 348             eello_turn6=enetb(10,i,iparm)
 349             ebe=enetb(11,i,iparm)
 350             escloc=enetb(12,i,iparm)
 351             etors=enetb(13,i,iparm)
 352             etors_d=enetb(14,i,iparm)
 353             ehpb=enetb(15,i,iparm)
 354 !            estr=enetb(18,i,iparm)
 355             estr=enetb(17,i,iparm)
 356 !            esccor=enetb(19,i,iparm)
 357             esccor=enetb(21,i,iparm)
 358 !            edihcnstr=enetb(20,i,iparm)
 359             edihcnstr=enetb(19,i,iparm)
 360           ecationcation=enetb(41,i,iparm)
 361           ecation_prot=enetb(42,i,iparm)
 362             evdwpp =      enetb(26,i,iparm)
 363             eespp =      enetb(27,i,iparm)
 364             evdwpsb =      enetb(28,i,iparm)
 365             eelpsb =      enetb(29,i,iparm)
 366             evdwsb =      enetb(30,i,iparm)
 367             eelsb =      enetb(31,i,iparm)
 368             estr_nucl =      enetb(32,i,iparm)
 369             ebe_nucl =      enetb(33,i,iparm)
 370             esbloc =      enetb(34,i,iparm)
 371             etors_nucl =      enetb(35,i,iparm)
 372             etors_d_nucl =      enetb(36,i,iparm)
 373             ecorr_nucl =      enetb(37,i,iparm)
 374             ecorr3_nucl =      enetb(38,i,iparm)
 375             epeppho=   enetb(49,i,iparm)
 376             escpho=    enetb(48,i,iparm)
 377             epepbase=  enetb(47,i,iparm)
 378             escbase=   enetb(46,i,iparm)
 379             ecation_nucl= enetb(50,i,iparm)
 380             elipbond=enetb(52,i,iparm)
 381             elipang=enetb(53,i,iparm)
 382             eliplj=enetb(54,i,iparm)
 383             elipelec=enetb(55,i,iparm)
 384             ecat_prottran=enetb(56,i,iparm)
 385             ecation_protang=enetb(57,i,iparm)
 386             eliptran=enetb(22,i,iparm)
 387             ehomology_constr=enetb(51,i,iparm)
 388 #ifdef DEBUG
 389             write (iout,'(3i5,6f5.2,15f12.3)') i,ib,iparm,(ft(l),l=1,6),&
 390              evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
 391              etors,etors_d,eello_turn3,eello_turn4,esccor,ecationcation,&
 392              ehomology_constr
 393 #endif
 394
 395 !#ifdef SPLITELE
 396 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
 397 !            +wvdwpp*evdw1 &
 398 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 399 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 400 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 401 !            +ft(2)*wturn3*eello_turn3 &
 402 !            +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
 403 !            +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 404 !            +wbond*estr
 405 !#else
 406 !            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
 407 !            +ft(1)*welec*(ees+evdw1) &
 408 !            +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 409 !            +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 410 !            +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 411 !            +ft(2)*wturn3*eello_turn3 &
 412 !            +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
 413 !            +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 414 !            +wbond*estr
 415 !#endif
 416
 417 #ifdef SPLITELE
 418             etot=wsc*evdw+wscp*evdw2+welec*ees &
 419             +wvdwpp*evdw1 &
 420             +wang*ebe+wtor*etors+wscloc*escloc &
 421             +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
 422             +wcorr6*ecorr6+wturn4*eello_turn4 &
 423             +wturn3*eello_turn3 &
 424             +wturn6*eello_turn6+wel_loc*eel_loc &
 425             +edihcnstr+wtor_d*etors_d+wsccor*esccor &
 426             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
 427             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 428             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 429             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 430             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
 431             +wscbase*escbase&
 432             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
 433             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
 434             +wliptran*eliptran
 435
 436
 437 #else
 438             etot=wsc*evdw+wscp*evdw2 &
 439             +welec*(ees+evdw1) &
 440             +wang*ebe+wtor*etors+wscloc*escloc &
 441             +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
 442             +wcorr6*ecorr6+wturn4*eello_turn4 &
 443             +wturn3*eello_turn3 &
 444             +wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
 445             +wtor_d*etors_d+wsccor*esccor &
 446             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
 447             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 448             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 449             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 450             +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
 451             +wscbase*escbase&
 452             +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
 453             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
 454             +wliptran*eliptran
 455
 456
 457
 458
 459 #endif
 460
 461 #ifdef DEBUG
 462             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),&
 463               etot,potEmin
 464 #endif
 465 #ifdef DEBUG
 466             if (iparm.eq.1 .and. ib.eq.1) then
 467             write (iout,*)"Conformation",i
 468             energia(0)=etot
 469             do k=1,max_ene
 470               energia(k)=enetb(k,i,iparm)
 471             enddo
 472 !            call enerprint(energia(0),fT)
 473             call enerprint(energia(0))
 474             endif
 475 #endif
 476             write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
 477             write (iout,*) "waga_homology", waga_homology
 478 #ifdef DEBUG
 479             write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
 480 #endif
 481             if (homol_nset.gt.1) then
 482
 483               do kk=1,nR(ib,iparm)
 484                 Econstr=waga_homology(kk)*ehomology_constr
 485                 v(i,kk,ib,iparm)= &
 486                  -beta_h(ib,iparm)*(etot+Econstr)
 487 #ifdef DEBUG
 488                 write (iout,'(4i5,4e15.5)') i,kk,ib,iparm, &
 489                 etot,Econstr,v(i,kk,ib,iparm)
 490 #endif
 491               enddo ! kk
 492
 493             else
 494
 495               etot=etot+ehomology_constr
 496
 497
 498             do kk=1,nR(ib,iparm)
 499               Econstr=0.0d0
 500               do j=1,nQ
 501                 ddW = q(j,i)
 502                 Econstr=Econstr+Kh(j,kk,ib,iparm) &
 503                  *(ddW-q0(j,kk,ib,iparm))**2
 504               enddo
 505               v(i,kk,ib,iparm)= &
 506                 -beta_h(ib,iparm)*(etot-potEmin+Econstr)
 507 #ifdef DEBUG
 508               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,&
 509                etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 510 #endif
 511             enddo ! kk
 512            endif
 513           enddo   ! ib
 514         enddo     ! iparm
 515       enddo       ! i
 516 ! Simple iteration to calculate free energies corresponding to all simulation
 517 ! runs.
 518       do iter=1,maxit
 519
 520 ! Compute new free-energy values corresponding to the righ-hand side of the
 521 ! equation and their derivatives.
 522         write (iout,*) "------------------------fi"
 523 #ifdef MPI
 524         do t=1,scount(me1)
 525 #else
 526         do t=1,ntot(islice)
 527 #endif
 528           vmax=-1.0d+20
 529           do i=1,nParmSet
 530             do k=1,nT_h(i)
 531               do l=1,nR(k,i)
 532                 vf=v(t,l,k,i)+f(l,k,i)
 533                 if (vf.gt.vmax) vmax=vf
 534               enddo
 535             enddo
 536           enddo
 537           denom=0.0d0
 538           do i=1,nParmSet
 539             do k=1,nT_h(i)
 540               do l=1,nR(k,i)
 541                 aux=f(l,k,i)+v(t,l,k,i)-vmax
 542                 if (aux.gt.-200.0d0) &
 543                 denom=denom+snk(l,k,i,islice)*dexp(aux)
 544               enddo
 545             enddo
 546           enddo
 547           entfac(t)=-dlog(denom)-vmax
 548           if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 549 #ifdef DEBUG
 550           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 551 #endif
 552         enddo
 553         do iparm=1,nParmSet
 554           do iib=1,nT_h(iparm)
 555             do ii=1,nR(iib,iparm)
 556 #ifdef MPI
 557               fi_p_min(ii,iib,iparm)=-1.0d5
 558               do t=1,scount(me)
 559 !                fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
 560 !                  +dexp(v(t,ii,iib,iparm)+entfac(t))
 561                 aux=v(t,ii,iib,iparm)+entfac(t)
 562                 if (aux.gt.fi_p_min(ii,iib,iparm))   fi_p_min(ii,iib,iparm)=aux
 563
 564 !#define DEBUG
 565 #ifdef DEBUG
 566               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,&
 567                v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm),fi_p_min(ii,iib,iparm)
 568 #endif
 569 !#undef DEBUG
 570               enddo
 571 #else
 572 !              fi(ii,iib,iparm)=0.0d0
 573               do t=1,ntot(islice)
 574 !                fi(ii,iib,iparm)=fi(ii,iib,iparm) &
 575 !                  +dexp(v(t,ii,iib,iparm)+entfac(t))
 576                 aux=v(t,ii,iib,iparm)+entfac(t)
 577                 if (aux.gt.fi_min(ii,iib,iparm))
 578      &                                     fi_min(ii,iib,iparm)=aux
 579
 580               enddo
 581 #endif
 582             enddo ! ii
 583           enddo ! iib
 584         enddo ! iparm
 585
 586 #ifdef MPI
 587 #ifdef DEBUG
 588         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 589         do iparm=1,nParmSet
 590           do ib=1,nT_h(nparmset)
 591             write (iout,*) "iparm",iparm," ib",ib
 592             write (iout,*) "beta=",beta_h(ib,iparm)
 593             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 594             write (iout,'(8e15.5)') (fi_p_min(i,ib,iparm),i=1,nR(ib,iparm))
 595           enddo
 596         enddo
 597 #endif
 598         call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet, &
 599               MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
 600
 601 #ifdef DEBUG
 602         write (iout,*) "fi_min after AllReduce"
 603         do i=1,nParmSet
 604           do j=1,nT_h(i)
 605             write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
 606           enddo
 607         enddo
 608 #endif
 609 !#undef DEBUG
 610 #endif
 611         do iparm=1,nParmSet
 612           do iib=1,nT_h(iparm)
 613             do ii=1,nR(iib,iparm)
 614 #ifdef MPI
 615               fi_p(ii,iib,iparm)=0.0d0
 616               do t=1,scount(me)
 617                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
 618                 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 619 #ifdef DEBUG
 620               write (iout,'(4i5,4e15.5)') t,ii,iib,iparm, &
 621               v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm), &
 622               fi_p(ii,iib,iparm)
 623 #endif
 624               enddo
 625 #else
 626               fi(ii,iib,iparm)=0.0d0
 627               do t=1,ntot(islice)
 628                 fi(ii,iib,iparm)=fi(ii,iib,iparm) &
 629                 +dexp(v(t,ii,iib,iparm)+entfac(t)-fi_min(ii,iib,iparm))
 630               enddo
 631 #endif
 632             enddo ! ii
 633           enddo ! iib
 634         enddo ! iparm
 635
 636 #ifdef MPI
 637 #ifdef DEBUG
 638         write (iout,*) "fi before MPI_Reduce me",me,' master',master
 639         do iparm=1,nParmSet
 640           do ib=1,nT_h(nparmset)
 641             write (iout,*) "iparm",iparm," ib",ib
 642             write (iout,*) "beta=",beta_h(ib,iparm)
 643             write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
 644           enddo
 645         enddo
 646 #endif
 647 #ifdef DEBUG
 648         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet, &
 649         maxR*MaxT_h*nParmSet
 650         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD, &
 651         " WHAM_COMM",WHAM_COMM
 652 #endif
 653
 654         write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,&
 655          maxR*MaxT_h*nParmSet
 656         write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
 657          " WHAM_COMM",WHAM_COMM
 658         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,&
 659          MPI_DOUBLE_PRECISION,&
 660          MPI_SUM,Master,WHAM_COMM,IERROR)
 661 #ifdef DEBUG
 662         write (iout,*) "fi after MPI_Reduce nparmset",nparmset
 663         do iparm=1,nParmSet
 664           write (iout,*) "iparm",iparm
 665           do ib=1,nT_h(iparm)
 666             write (iout,*) "beta=",beta_h(ib,iparm)
 667             write (iout,'(8e15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 668           enddo
 669         enddo
 670 #endif
 671         if (me1.eq.Master) then
 672 #endif
 673         avefi=0.0d0
 674         do iparm=1,nParmSet
 675           do ib=1,nT_h(iparm)
 676             do i=1,nR(ib,iparm)
 677               fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm)
 678               avefi=avefi+fi(i,ib,iparm)
 679             enddo
 680           enddo
 681         enddo
 682         avefi=avefi/nFi
 683         do iparm=1,nParmSet
 684           write (iout,*) "Parameter set",iparm
 685           do ib =1,nT_h(iparm)
 686             write (iout,*) "beta=",beta_h(ib,iparm)
 687             do i=1,nR(ib,iparm)
 688               fi(i,ib,iparm)=fi(i,ib,iparm)-avefi
 689             enddo
 690             write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm))
 691             write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm))
 692           enddo
 693         enddo
 694
 695 ! Compute the norm of free-energy increments.
 696         finorm=0.0d0
 697         do iparm=1,nParmSet
 698           do ib=1,nT_h(iparm)
 699             do i=1,nR(ib,iparm)
 700               finorm=finorm+dabs(fi(i,ib,iparm)-f(i,ib,iparm))
 701               f(i,ib,iparm)=fi(i,ib,iparm)
 702             enddo
 703           enddo
 704         enddo
 705
 706         write (iout,*) 'Iteration',iter,' finorm',finorm
 707
 708 #ifdef MPI
 709         endif
 710         call MPI_Bcast(f(1,1,1),MaxR*MaxT_h*nParmSet,&
 711          MPI_DOUBLE_PRECISION,Master,&
 712          WHAM_COMM,IERROR)
 713         call MPI_Bcast(finorm,1,MPI_DOUBLE_PRECISION,Master,&
 714          WHAM_COMM,IERROR)
 715 #endif
 716 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 717         if (finorm.lt.fimin) then
 718           write (iout,*) 'Iteration converged'
 719           goto 20
 720         endif
 721
 722       enddo ! iter
 723
 724    20 continue
 725 ! Now, put together the histograms from all simulations, in order to get the
 726 ! unbiased total histogram.
 727 !C Determine the minimum free energies
 728 #ifdef MPI
 729       do i=1,scount(me1)
 730 #else
 731       do i=1,ntot(islice)
 732 #endif
 733 !c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
 734         do iparm=1,nParmSet
 735 #ifdef DEBUG
 736           write (iout,'(2i5,21f8.2)') i,iparm, &
 737           (enetb(k,i,iparm),k=1,26)
 738 #endif
 739           call restore_parm(iparm)
 740 #ifdef DEBUG
 741           write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc, &
 742            wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, &
 743            wtor_d,wsccor,wbond
 744 #endif
 745           do ib=1,nT_h(iparm)
 746             if (rescale_modeW.eq.1) then
 747               quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 748               quotl=1.0d0
 749               kfacl=1.0d0
 750               do l=1,5
 751                 quotl1=quotl
 752                 quotl=quotl*quot
 753                 kfacl=kfacl*kfac
 754                 fT(l)=kfacl/(kfacl-1.0d0+quotl)
 755               enddo
 756 #if defined(FUNCTH)
 757               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 758               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 759 #elif defined(FUNCT)
 760               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 761 #else
 762               ft(6)=1.0d0
 763 #endif
 764             else if (rescale_modeW.eq.2) then
 765               quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
 766               quotl=1.0d0
 767               do l=1,5
 768                 quotl=quotl*quot
 769                 fT(l)=1.12692801104297249644d0/ &
 770                   dlog(dexp(quotl)+dexp(-quotl))
 771                enddo
 772 #if defined(FUNCTH)
 773               tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
 774               ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
 775 #elif defined(FUNCT)
 776               ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
 777 #else
 778               ft(6)=1.0d0
 779 #endif
 780 !c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 781             else if (rescale_modeW.eq.0) then
 782               do l=1,6
 783                 fT(l)=1.0d0
 784               enddo
 785             else
 786               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
 787                rescale_mode
 788               call flush(iout)
 789               return 1
 790             endif
 791             evdw=enetb(1,i,iparm)
 792 !            evdw_t=enetb(21,i,iparm)
 793             evdw_t=enetb(20,i,iparm)
 794 #ifdef SCP14
 795 !            evdw2_14=enetb(17,i,iparm)
 796             evdw2_14=enetb(18,i,iparm)
 797             evdw2=enetb(2,i,iparm)+evdw2_14
 798 #else
 799             evdw2=enetb(2,i,iparm)
 800             evdw2_14=0.0d0
 801 #endif
 802 #ifdef SPLITELE
 803             ees=enetb(3,i,iparm)
 804             evdw1=enetb(16,i,iparm)
 805 #else
 806             ees=enetb(3,i,iparm)
 807             evdw1=0.0d0
 808 #endif
 809             ecorr=enetb(4,i,iparm)
 810             ecorr5=enetb(5,i,iparm)
 811             ecorr6=enetb(6,i,iparm)
 812             eel_loc=enetb(7,i,iparm)
 813             eello_turn3=enetb(8,i,iparm)
 814             eello_turn4=enetb(9,i,iparm)
 815             eello_turn6=enetb(10,i,iparm)
 816             ebe=enetb(11,i,iparm)
 817             escloc=enetb(12,i,iparm)
 818             etors=enetb(13,i,iparm)
 819             etors_d=enetb(14,i,iparm)
 820             ehpb=enetb(15,i,iparm)
 821 !            estr=enetb(18,i,iparm)
 822             estr=enetb(17,i,iparm)
 823 !            esccor=enetb(19,i,iparm)
 824             esccor=enetb(21,i,iparm)
 825 !            edihcnstr=enetb(20,i,iparm)
 826             edihcnstr=enetb(19,i,iparm)
 827 !            ehomology_constr=enetb(22,i,iparm)
 828 !            esaxs=enetb(26,i,iparm)
 829           ecationcation=enetb(41,i,iparm)
 830           ecation_prot=enetb(42,i,iparm)
 831             evdwpp =      enetb(26,i,iparm)
 832             eespp =      enetb(27,i,iparm)
 833             evdwpsb =      enetb(28,i,iparm)
 834             eelpsb =      enetb(29,i,iparm)
 835             evdwsb =      enetb(30,i,iparm)
 836             eelsb =      enetb(31,i,iparm)
 837             estr_nucl =      enetb(32,i,iparm)
 838             ebe_nucl =      enetb(33,i,iparm)
 839             esbloc =      enetb(34,i,iparm)
 840             etors_nucl =      enetb(35,i,iparm)
 841             etors_d_nucl =      enetb(36,i,iparm)
 842             ecorr_nucl =      enetb(37,i,iparm)
 843             ecorr3_nucl =      enetb(38,i,iparm)
 844             epeppho=   enetb(49,i,iparm)
 845             escpho=    enetb(48,i,iparm)
 846             epepbase=  enetb(47,i,iparm)
 847             escbase=   enetb(46,i,iparm)
 848             ecation_nucl=enetb(50,i,iparm)
 849             elipbond=enetb(52,i,iparm)
 850             elipang=enetb(53,i,iparm)
 851             eliplj=enetb(54,i,iparm)
 852             elipelec=enetb(55,i,iparm)
 853             ecat_prottran=enetb(56,i,iparm)
 854             ecation_protang=enetb(57,i,iparm)
 855             eliptran=enetb(22,i,iparm)
 856             ehomology_constr=enetb(51,i,iparm)
 857 #ifdef SPLITELE
 858             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
 859             +wvdwpp*evdw1 &
 860             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 861             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 862             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 863             +ft(2)*wturn3*eello_turn3 &
 864             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
 865             +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 866             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
 867             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 868             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 869             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 870             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
 871             +wcorr3_nucl*ecorr3_nucl&
 872             +wscbase*escbase&
 873             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
 874             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
 875             +wliptran*eliptran
 876
 877
 878 #else
 879             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
 880             +ft(1)*welec*(ees+evdw1) &
 881             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
 882             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
 883             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
 884             +ft(2)*wturn3*eello_turn3 &
 885             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
 886             +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 887             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
 888             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
 889             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
 890             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
 891             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
 892             +wcorr3_nucl*ecorr3_nucl&
 893             +wscbase*escbase&
 894             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
 895             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
 896             +wliptran*eliptran
 897
 898
 899
 900 #endif
 901             write (iout,*) "WTF,",etot,potEmin_all(ib,iparm),entfac(i)/beta_h(ib,iparm)
 902             etot=etot-entfac(i)/beta_h(ib,iparm)
 903             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
 904
 905           enddo   ! ib
 906         enddo     ! iparm
 907       enddo       ! i
 908 #ifdef DEBUG
 909       write (iout,*) "The potEmin array before reduction"
 910       do i=1,nParmSet
 911         write (iout,*) "Parameter set",i
 912         do j=1,nT_h(i)
 913           write (iout,*) j,PotEmin_all(j,i)
 914         enddo
 915       enddo
 916       write (iout,*) "potEmin_min",potEmin_min
 917 #endif
 918 #ifdef MPI
 919 !C Determine the minimum energes for all parameter sets and temperatures
 920       call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1), &
 921        maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
 922       do i=1,nParmSet
 923         do j=1,nT_h(i)
 924           potEmin_all(j,i)=potEmin_t_all(j,i)
 925         enddo
 926       enddo
 927 #endif
 928       potEmin_min=potEmin_all(1,1)
 929       do i=1,nParmSet
 930         do j=1,nT_h(i)
 931           if (potEmin_all(j,i).lt.potEmin_min) &
 932                   potEmin_min=potEmin_all(j,i)
 933         enddo
 934       enddo
 935 #ifdef DEBUG
 936       write (iout,*) "The potEmin array"
 937       do i=1,nParmSet
 938         write (iout,*) "Parameter set",i
 939         do j=1,nT_h(i)
 940           write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)), &
 941            PotEmin_all(j,i)
 942         enddo
 943       enddo
 944       write (iout,*) "potEmin_min",potEmin_min
 945 #endif
 946
 947 #ifdef MPI
 948       do t=0,tmax
 949         hfin_ent_p(t)=0.0d0
 950       enddo
 951 #else
 952       do t=0,tmax
 953         hfin_ent(t)=0.0d0
 954       enddo
 955 #endif
 956       write (iout,*) "--------------hist"
 957 #ifdef MPI
 958       do iparm=1,nParmSet
 959         do i=0,nGridT
 960           sumW_p(i,iparm)=0.0d0
 961           sumE_p(i,iparm)=0.0d0
 962           sumEbis_p(i,iparm)=0.0d0
 963           sumEsq_p(i,iparm)=0.0d0
 964           do j=1,nQ+2
 965             sumQ_p(j,i,iparm)=0.0d0
 966             sumQsq_p(j,i,iparm)=0.0d0
 967             sumEQ_p(j,i,iparm)=0.0d0
 968           enddo
 969         enddo
 970       enddo
 971       upindE_p=0
 972 #else
 973       do iparm=1,nParmSet
 974         do i=0,nGridT
 975           sumW(i,iparm)=0.0d0
 976           sumE(i,iparm)=0.0d0
 977           sumEbis(i,iparm)=0.0d0
 978           sumEsq(i,iparm)=0.0d0
 979           do j=1,nQ+2
 980             sumQ(j,i,iparm)=0.0d0
 981             sumQsq(j,i,iparm)=0.0d0
 982             sumEQ(j,i,iparm)=0.0d0
 983           enddo
 984         enddo
 985       enddo
 986       upindE=0
 987 #endif
 988 ! 8/26/05 entropy distribution
 989 !#ifdef MPI
 990 !      entmin_p=1.0d10
 991 !      entmax_p=-1.0d10
 992 !      do t=1,scount(me1)
 993 !!        ent=-dlog(entfac(t))
 994 !        ent=entfac(t)
 995 !        if (ent.lt.entmin_p) entmin_p=ent
 996 !        if (ent.gt.entmax_p) entmax_p=ent
 997 !      enddo
 998 !      write (iout,*) "entmin",entmin_p," entmax",entmax_p
 999 !!      write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
1000 !      call flush(iout)
1001 !      call MPI_Allreduce(entmin_p,entmin,1,MPI_DOUBLE_PRECISION,MPI_MIN,&
1002 !        WHAM_COMM,IERROR)
1003 !      call MPI_Allreduce(entmax_p,entmax,1,MPI_DOUBLE_PRECISION,MPI_MAX,&
1004 !        WHAM_COMM,IERROR)
1005 !      write (iout,*) "entmin",entmin," entmax",entmax
1006 !      write (iout,*) "entmin_p",entmin_p," entmax_p",entmax_p
1007 !      ientmax=entmax-entmin
1008 !iientmax=entmax-entmin !el
1009 !write (iout,*) "ientmax",ientmax,entmax,entmin
1010 !write (iout,*) "iientmax",iientmax
1011 !      if (ientmax.gt.2000) ientmax=2000
1012 !      write (iout,*) "entmin",entmin," entmax",entmax," ientmax",ientmax
1013 !      call flush(iout)
1014 !      do t=1,scount(me1)
1015 !!        ient=-dlog(entfac(t))-entmin
1016 !        ient=entfac(t)-entmin
1017 !        if (ient.le.2000) histent_p(ient)=histent_p(ient)+1
1018 !      enddo
1019 !      call MPI_Allreduce(histent_p(0),histent(0),ientmax+1,MPI_INTEGER,&
1020 !        MPI_SUM,WHAM_COMM,IERROR)
1021 !      if (me1.eq.Master) then
1022 !        write (iout,*) "Entropy histogram"
1023 !        do i=0,ientmax
1024 !          write(iout,'(f15.4,i10)') entmin+i,histent(i)
1025 !        enddo
1026 !      endif
1027 !#else
1028 !      entmin=1.0d10
1029 !      entmax=-1.0d10
1030 !      do t=1,ntot(islice)
1031 !        ent=entfac(t)
1032 !        if (ent.lt.entmin) entmin=ent
1033 !        if (ent.gt.entmax) entmax=ent
1034 !      enddo
1035 !      ientmax=-dlog(entmax)-entmin
1036 !      if (ientmax.gt.2000) ientmax=2000
1037 !      do t=1,ntot(islice)
1038 !        ient=entfac(t)-entmin
1039 !        if (ient.le.2000) histent(ient)=histent(ient)+1
1040 !      enddo
1041 !      write (iout,*) "Entropy histogram"
1042 !      do i=0,ientmax
1043 !        write(iout,'(2f15.4)') entmin+i,histent(i)
1044 !      enddo
1045 !#endif
1046
1047 #ifdef MPI
1048       write (iout,*) "me1",me1," scount",scount(me1) !d
1049
1050       do iparm=1,nParmSet
1051
1052 #ifdef MPI
1053         do ib=1,nT_h(iparm)
1054           do t=0,tmax
1055             hfin_p(t,ib)=0.0d0
1056           enddo
1057         enddo
1058         do i=1,maxindE
1059           histE_p(i)=0.0d0
1060         enddo
1061 #else
1062         do ib=1,nT_h(iparm)
1063           do t=0,tmax
1064             hfin(t,ib)=0.0d0
1065           enddo
1066         enddo
1067         do i=1,maxindE
1068           histE(i)=0.0d0
1069         enddo
1070 #endif
1071         do ib=1,nT_h(iparm)
1072           do i=0,MaxBinRms
1073             do j=0,MaxBinRgy
1074               hrmsrgy(j,i,ib)=0.0d0
1075 #ifdef MPI
1076               hrmsrgy_p(j,i,ib)=0.0d0
1077 #endif
1078             enddo
1079           enddo
1080         enddo
1081
1082         do t=1,scount(me1)
1083 #else
1084         do t=1,ntot(islice)
1085 #endif
1086           ind = ind_point(t)
1087 #ifdef MPI
1088           hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
1089 #else
1090           hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
1091 #endif
1092 !          write (iout,'(2i5,20f8.2)') "debug",t,t,(enetb(k,t,iparm),k=1,21)
1093           call restore_parm(iparm)
1094 !          evdw=enetb(21,t,iparm)
1095           evdw=enetb(20,t,iparm)
1096           evdw_t=enetb(1,t,iparm)
1097 #ifdef SCP14
1098 !          evdw2_14=enetb(17,t,iparm)
1099           evdw2_14=enetb(18,t,iparm)
1100           evdw2=enetb(2,t,iparm)+evdw2_14
1101 #else
1102           evdw2=enetb(2,t,iparm)
1103           evdw2_14=0.0d0
1104 #endif
1105 #ifdef SPLITELE
1106           ees=enetb(3,t,iparm)
1107           evdw1=enetb(16,t,iparm)
1108 #else
1109           ees=enetb(3,t,iparm)
1110           evdw1=0.0d0
1111 #endif
1112           ecorr=enetb(4,t,iparm)
1113           ecorr5=enetb(5,t,iparm)
1114           ecorr6=enetb(6,t,iparm)
1115           eel_loc=enetb(7,t,iparm)
1116           eello_turn3=enetb(8,t,iparm)
1117           eello_turn4=enetb(9,t,iparm)
1118           eello_turn6=enetb(10,t,iparm)
1119           ebe=enetb(11,t,iparm)
1120           escloc=enetb(12,t,iparm)
1121           etors=enetb(13,t,iparm)
1122           etors_d=enetb(14,t,iparm)
1123           ehpb=enetb(15,t,iparm)
1124 !          estr=enetb(18,t,iparm)
1125           estr=enetb(17,t,iparm)
1126 !          esccor=enetb(19,t,iparm)
1127           esccor=enetb(21,t,iparm)
1128 !          edihcnstr=enetb(20,t,iparm)
1129           edihcnstr=enetb(19,t,iparm)
1130           edihcnstr=0.0d0
1131           ecationcation=enetb(41,t,iparm)
1132           ecation_prot=enetb(42,t,iparm)
1133             evdwpp =      enetb(26,t,iparm)
1134             eespp =      enetb(27,t,iparm)
1135             evdwpsb =      enetb(28,t,iparm)
1136             eelpsb =      enetb(29,t,iparm)
1137             evdwsb =      enetb(30,t,iparm)
1138             eelsb =      enetb(31,t,iparm)
1139             estr_nucl =      enetb(32,t,iparm)
1140             ebe_nucl =      enetb(33,t,iparm)
1141             esbloc =      enetb(34,t,iparm)
1142             etors_nucl =      enetb(35,t,iparm)
1143             etors_d_nucl =      enetb(36,t,iparm)
1144             ecorr_nucl =      enetb(37,t,iparm)
1145             ecorr3_nucl =      enetb(38,t,iparm)
1146             epeppho=   enetb(49,t,iparm)
1147             escpho=    enetb(48,t,iparm)
1148             epepbase=  enetb(47,t,iparm)
1149             escbase=   enetb(46,t,iparm)
1150             ecation_nucl=enetb(50,t,iparm)
1151             elipbond=enetb(52,t,iparm)
1152             elipang=enetb(53,t,iparm)
1153             eliplj=enetb(54,t,iparm)
1154             elipelec=enetb(55,t,iparm)
1155             ecat_prottran=enetb(56,t,iparm)
1156             ecation_protang=enetb(57,t,iparm)
1157             eliptran=enetb(22,t,iparm)
1158             ehomology_constr=enetb(51,t,iparm)
1159
1160           do k=0,nGridT
1161             betaT=startGridT+k*delta_T
1162             temper=betaT
1163 !write(iout,*)"kkkkkkkk",betaT,startGridT,k,delta_T
1164 !d            fT=T0/betaT
1165 !d            ft=2*T0/(T0+betaT)
1166             if (rescale_modeW.eq.1) then
1167               quot=betaT/T0
1168               quotl=1.0d0
1169               kfacl=1.0d0
1170               do l=1,5
1171                 quotl1=quotl
1172                 quotl=quotl*quot
1173                 kfacl=kfacl*kfac
1174                 denom=kfacl-1.0d0+quotl
1175                 fT(l)=kfacl/denom
1176                 ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1177                 ftbis(l)=l*kfacl*quotl1* &
1178                  (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1179               enddo
1180 #if defined(FUNCTH)
1181               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1182                         320.0d0
1183               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1184              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1185                     /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1186 #elif defined(FUNCT)
1187               fT(6)=betaT/T0
1188               ftprim(6)=1.0d0/T0
1189               ftbis(6)=0.0d0
1190 #else
1191               fT(6)=1.0d0
1192               ftprim(6)=0.0d0
1193               ftbis(6)=0.0d0
1194 #endif
1195             else if (rescale_modeW.eq.2) then
1196               quot=betaT/T0
1197               quotl=1.0d0
1198               do l=1,5
1199                 quotl1=quotl
1200                 quotl=quotl*quot
1201                 eplus=dexp(quotl)
1202                 eminus=dexp(-quotl)
1203                 logfac=1.0d0/dlog(eplus+eminus)
1204                 tanhT=(eplus-eminus)/(eplus+eminus)
1205                 fT(l)=1.12692801104297249644d0*logfac
1206                 ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1207                 ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1208                 2*l*quotl1/T0*logfac* &
1209                 (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1210                 +ftprim(l)*tanhT)
1211               enddo
1212 #if defined(FUNCTH)
1213               ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1214                        320.0d0
1215               ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1216              ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1217                      /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1218 #elif defined(FUNCT)
1219               fT(6)=betaT/T0
1220               ftprim(6)=1.0d0/T0
1221               ftbis(6)=0.0d0
1222 #else
1223               fT(6)=1.0d0
1224               ftprim(6)=0.0d0
1225               ftbis(6)=0.0d0
1226 #endif
1227             else if (rescale_modeW.eq.0) then
1228               do l=1,5
1229                 fT(l)=1.0d0
1230                 ftprim(l)=0.0d0
1231               enddo
1232             else
1233               write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",&
1234                 rescale_modeW
1235               call flush(iout)
1236               return 1
1237             endif
1238 !            write (iout,*) "ftprim",ftprim
1239 !            write (iout,*) "ftbis",ftbis
1240             betaT=1.0d0/(1.987D-3*betaT)
1241 #ifdef SPLITELE
1242             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1243             +wvdwpp*evdw1 &
1244             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1245             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1246             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1247             +ft(2)*wturn3*eello_turn3 &
1248             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1249             +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1250             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1251             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1252             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1253             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1254             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1255             +wcorr3_nucl*ecorr3_nucl&
1256             +wscbase*escbase&
1257             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1258             + wcatnucl*ecation_nucl&
1259             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1260             +wliptran*eliptran
1261
1262
1263
1264             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees &
1265                  +ftprim(1)*wtor*etors+ &
1266                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1267                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1268                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1269                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1270                  ftprim(1)*wsccor*esccor +ftprim(1)*wtor_nucl*etors_nucl&
1271                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1272             ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+ &
1273                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1274                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1275                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1276                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1277                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1278                  +wtor_d_nucl*ftbis(2)*etors_d_nucl+ftbis(1)*wpepbase*epepbase
1279 #else
1280             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1281             +ft(1)*welec*(ees+evdw1) &
1282             +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1283             +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1284             +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1285             +ft(2)*wturn3*eello_turn3 &
1286             +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1287             +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1288             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
1289             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1290             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1291             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1292             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1293             +wcorr3_nucl*ecorr3_nucl&
1294             +wscbase*escbase&
1295             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
1296             + wcatnucl*ecation_nucl&
1297             +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
1298             +wliptran*eliptran
1299
1300
1301
1302             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
1303                 +ftprim(1)*wtor*etors+ &
1304                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1305                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1306                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1307                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1308                  ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
1309                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1310             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
1311                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1312                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1313                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1314                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1315                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1316                  +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
1317
1318 #endif
1319 !            weight=dexp(-betaT*(etot-potEmin)+entfac(t))
1320 !#ifdef DEBUG
1321 !            write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
1322 !              " etot",etot," entfac",entfac(t),&
1323 !              " weight",weight," ebis",ebis
1324 !#endif
1325 !            etot=etot-temper*eprim
1326 !#ifdef MPI
1327 !            sumW_p(k,iparm)=sumW_p(k,iparm)+weight
1328 !            sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
1329 !            sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
1330 !            sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
1331 !            do j=1,nQ+2
1332 !              sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
1333 !              sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
1334 !              sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
1335 !               +etot*q(j,t)*weight
1336 !            enddo
1337 !#else
1338 !            sumW(k,iparm)=sumW(k,iparm)+weight
1339 !            sumE(k,iparm)=sumE(k,iparm)+etot*weight
1340 !            sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
1341 !            sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
1342 !            do j=1,nQ+2
1343 !              sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
1344 !              sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
1345 !              sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
1346 !               +etot*q(j,t)*weight
1347 !            enddo
1348 !#endif
1349           enddo
1350           indE = aint(potE(t,iparm)-aint(potEmin))
1351           if (indE.ge.0 .and. indE.le.maxinde) then
1352             if (indE.gt.upindE_p) upindE_p=indE
1353             histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
1354           endif
1355 #ifdef MPI
1356           do ib=1,nT_h(iparm)
1357             potEmin=potEmin_all(ib,iparm)
1358             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1359             hfin_p(ind,ib)=hfin_p(ind,ib)+ &
1360              dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1361              if (rmsrgymap) then
1362                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1363                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1364                hrmsrgy_p(indrgy,indrms,ib)= &
1365                  hrmsrgy_p(indrgy,indrms,ib)+expfac
1366              endif
1367           enddo
1368 #else
1369           do ib=1,nT_h(iparm)
1370             potEmin=potEmin_all(ib,iparm)
1371             expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1372             hfin(ind,ib)=hfin(ind,ib)+ &
1373              dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
1374              if (rmsrgymap) then
1375                indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
1376                indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
1377                hrmsrgy(indrgy,indrms,ib)= &
1378                  hrmsrgy(indrgy,indrms,ib)+expfac
1379              endif
1380           enddo
1381 #endif
1382         enddo ! t
1383 !
1384 ! Thermo and ensemble averages
1385 !
1386         do k=0,nGridT
1387           betaT=startGridT+k*delta_T
1388 !          call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
1389       if (rescale_mode.eq.1) then
1390         quot=betaT/T0
1391         quotl=1.0d0
1392         kfacl=1.0d0
1393         do l=1,5
1394           quotl1=quotl
1395           quotl=quotl*quot
1396           kfacl=kfacl*kfac
1397           denom=kfacl-1.0d0+quotl
1398           fT(l)=kfacl/denom
1399           ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1400           ftbis(l)=l*kfacl*quotl1* &
1401           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1402         enddo
1403 #if defined(FUNCTH)
1404         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1405                  320.0d0
1406         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1407         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1408                /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1409 #elif defined(FUNCT)
1410         fT(6)=betaT/T0
1411         ftprim(6)=1.0d0/T0
1412         ftbis(6)=0.0d0
1413 #else
1414         fT(6)=1.0d0
1415         ftprim(6)=0.0d0
1416         ftbis(6)=0.0d0
1417 #endif
1418       else if (rescale_mode.eq.2) then
1419         quot=betaT/T0
1420         quotl=1.0d0
1421         do l=1,5
1422           quotl1=quotl
1423           quotl=quotl*quot
1424           eplus=dexp(quotl)
1425           eminus=dexp(-quotl)
1426           logfac=1.0d0/dlog(eplus+eminus)
1427           tanhT=(eplus-eminus)/(eplus+eminus)
1428           fT(l)=1.12692801104297249644d0*logfac
1429           ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1430           ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1431          2*l*quotl1/T0*logfac* &
1432          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1433          +ftprim(l)*tanhT)
1434         enddo
1435 #if defined(FUNCTH)
1436         ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/ &
1437                 320.0d0
1438         ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
1439         ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0) &
1440               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
1441 #elif defined(FUNCT)
1442         fT(6)=betaT/T0
1443         ftprim(6)=1.0d0/T0
1444         ftbis(6)=0.0d0
1445 #else
1446         fT(6)=1.0d0
1447         ftprim(6)=0.0d0
1448         ftbis(6)=0.0d0
1449 #endif
1450       else if (rescale_mode.eq.0) then
1451         do l=1,5
1452           fT(l)=1.0d0
1453           ftprim(l)=0.0d0
1454         enddo
1455       else
1456         write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
1457          rescale_mode
1458         call flush(iout)
1459 !        return1
1460       endif
1461
1462 !            write (iout,*) "ftprim",ftprim
1463 !            write (iout,*) "ftbis",ftbis
1464           betaT=1.0d0/(1.987D-3*betaT)
1465 ! 7/10/18 AL Determine the max Botzmann weights for each temerature
1466 !          call sum_ene(1,iparm,ft,etot)
1467       evdw=enetb(1,1,iparm)
1468       evdw_t=enetb(20,1,iparm)
1469 #ifdef SCP14
1470       evdw2_14=enetb(18,1,iparm)
1471       evdw2=enetb(2,1,iparm)+evdw2_14
1472 #else
1473       evdw2=enetb(2,1,iparm)
1474       evdw2_14=0.0d0
1475 #endif
1476 #ifdef SPLITELE
1477       ees=enetb(3,1,iparm)
1478       evdw1=enetb(16,1,iparm)
1479 #else
1480       ees=enetb(3,1,iparm)
1481       evdw1=0.0d0
1482 #endif
1483       ecorr=enetb(4,1,iparm)
1484       ecorr5=enetb(5,1,iparm)
1485       ecorr6=enetb(6,1,iparm)
1486       eel_loc=enetb(7,1,iparm)
1487       eello_turn3=enetb(8,1,iparm)
1488       eello_turn4=enetb(9,1,iparm)
1489       eello_turn6=enetb(10,1,iparm)
1490       ebe=enetb(11,1,iparm)
1491       escloc=enetb(12,1,iparm)
1492       etors=enetb(13,1,iparm)
1493       etors_d=enetb(14,1,iparm)
1494       ehpb=enetb(15,1,iparm)
1495       estr=enetb(17,1,iparm)
1496       esccor=enetb(21,1,iparm)
1497       edihcnstr=enetb(19,1,iparm)
1498 !      eliptran=enetb(22,1,iparm)
1499 !      esaxs=enetb(26,1,iparm)
1500       ecationcation=enetb(41,1,iparm)
1501       ecation_prot=enetb(42,1,iparm)
1502       evdwpp =      enetb(26,1,iparm)
1503       eespp =      enetb(27,1,iparm)
1504       evdwpsb =      enetb(28,1,iparm)
1505       eelpsb =      enetb(29,1,iparm)
1506       evdwsb =      enetb(30,1,iparm)
1507       eelsb =      enetb(31,1,iparm)
1508       estr_nucl =      enetb(32,1,iparm)
1509       ebe_nucl =      enetb(33,1,iparm)
1510       esbloc =      enetb(34,1,iparm)
1511       etors_nucl =      enetb(35,1,iparm)
1512       etors_d_nucl =      enetb(36,1,iparm)
1513       ecorr_nucl =      enetb(37,1,iparm)
1514       ecorr3_nucl =      enetb(38,1,iparm)
1515       epeppho=   enetb(49,1,iparm)
1516       escpho=    enetb(48,1,iparm)
1517       epepbase=  enetb(47,1,iparm)
1518       escbase=   enetb(46,1,iparm)
1519       ecation_nucl= enetb(50,1,iparm)
1520       ehomology_constr=enetb(51,1,iparm)
1521
1522 #ifdef SPLITELE
1523       if (shield_mode.gt.0) then
1524         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1525           +ft(1)*welec*ees &
1526           +ft(1)*wvdwpp*evdw1 &
1527           +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1528 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1529           +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1530           +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1531           +ft(2)*wturn3*eello_turn3 &
1532           +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1533           +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1534             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1535             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1536             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1537             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1538             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1539             +wcorr3_nucl*ecorr3_nucl&
1540             +wscbase*escbase&
1541             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1542             + wcatnucl*ecation_nucl
1543
1544
1545       else
1546         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1547           +wvdwpp*evdw1 &
1548           +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1549 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1550           +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1551           +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1552           +ft(2)*wturn3*eello_turn3 &
1553           +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1554           +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1555             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1556             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1557             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1558             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1559             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1560             +wcorr3_nucl*ecorr3_nucl&
1561             +wscbase*escbase&
1562             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1563             + wcatnucl*ecation_nucl
1564
1565       endif
1566 #else
1567       if (shield_mode.gt.0) then
1568         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1569           +ft(1)*welec*(ees+evdw1) &
1570           +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1571 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1572           +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1573           +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1574           +ft(2)*wturn3*eello_turn3 &
1575           +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1576           +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1577             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1578             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1579             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1580             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1581             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1582             +wcorr3_nucl*ecorr3_nucl&
1583             +wscbase*escbase&
1584             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1585             + wcatnucl*ecation_nucl
1586
1587       else
1588         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1589           +ft(1)*welec*(ees+evdw1) &
1590           +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1591 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1592           +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1593           +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1594           +ft(2)*wturn3*eello_turn3 &
1595           +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1596           +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1597             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1598             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1599             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1600             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1601             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1602             +wcorr3_nucl*ecorr3_nucl&
1603             +wscbase*escbase&
1604             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1605             + wcatnucl*ecation_nucl
1606
1607       endif
1608 #endif
1609
1610           weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
1611 !          write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
1612 #ifdef MPI
1613           do t=2,scount(me1)
1614 #else
1615           do t=2,ntot(islice)
1616 #endif
1617 !           call sum_ene(t,iparm,ft,etot)
1618       evdw=enetb(1,t,iparm)
1619       evdw_t=enetb(20,t,iparm)
1620 #ifdef SCP14
1621       evdw2_14=enetb(18,t,iparm)
1622       evdw2=enetb(2,t,iparm)+evdw2_14
1623 #else
1624       evdw2=enetb(2,t,iparm)
1625       evdw2_14=0.0d0
1626 #endif
1627 #ifdef SPLITELE
1628       ees=enetb(3,t,iparm)
1629       evdw1=enetb(16,t,iparm)
1630 #else
1631       ees=enetb(3,t,iparm)
1632       evdw1=0.0d0
1633 #endif
1634       ecorr=enetb(4,t,iparm)
1635       ecorr5=enetb(5,t,iparm)
1636       ecorr6=enetb(6,t,iparm)
1637       eel_loc=enetb(7,t,iparm)
1638       eello_turn3=enetb(8,t,iparm)
1639       eello_turn4=enetb(9,t,iparm)
1640       eello_turn6=enetb(10,t,iparm)
1641       ebe=enetb(11,t,iparm)
1642       escloc=enetb(12,t,iparm)
1643       etors=enetb(13,t,iparm)
1644       etors_d=enetb(14,t,iparm)
1645       ehpb=enetb(15,t,iparm)
1646       estr=enetb(17,t,iparm)
1647       esccor=enetb(21,t,iparm)
1648       edihcnstr=enetb(19,t,iparm)
1649 !      eliptran=enetb(22,t,iparm)
1650 !      esaxs=enetb(26,t,iparm)
1651       ecationcation=enetb(41,t,iparm)
1652       ecation_prot=enetb(42,t,iparm)
1653       evdwpp =      enetb(26,t,iparm)
1654       eespp =      enetb(27,t,iparm)
1655       evdwpsb =      enetb(28,t,iparm)
1656       eelpsb =      enetb(29,t,iparm)
1657       evdwsb =      enetb(30,t,iparm)
1658       eelsb =      enetb(31,t,iparm)
1659       estr_nucl =      enetb(32,t,iparm)
1660       ebe_nucl =      enetb(33,t,iparm)
1661       esbloc =      enetb(34,t,iparm)
1662       etors_nucl =      enetb(35,t,iparm)
1663       etors_d_nucl =      enetb(36,t,iparm)
1664       ecorr_nucl =      enetb(37,t,iparm)
1665       ecorr3_nucl =      enetb(38,t,iparm)
1666       epeppho=   enetb(49,t,iparm)
1667       escpho=    enetb(48,t,iparm)
1668       epepbase=  enetb(47,t,iparm)
1669       escbase=   enetb(46,t,iparm)
1670       ecation_nucl=enetb(50,t,iparm)
1671       ehomology_constr=enetb(51,t,iparm)
1672
1673 #ifdef SPLITELE
1674       if (shield_mode.gt.0) then
1675         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1676           +ft(1)*welec*ees &
1677           +ft(1)*wvdwpp*evdw1 &
1678           +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1679 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1680           +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1681           +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1682           +ft(2)*wturn3*eello_turn3 &
1683           +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1684           +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1685             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1686             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1687             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1688             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1689             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1690             +wcorr3_nucl*ecorr3_nucl&
1691             +wscbase*escbase&
1692             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1693             + wcatnucl*ecation_nucl
1694
1695       else
1696         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
1697           +wvdwpp*evdw1 &
1698           +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1699 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1700           +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1701           +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1702           +ft(2)*wturn3*eello_turn3 &
1703           +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1704           +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1705             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1706             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1707             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1708             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1709             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1710             +wcorr3_nucl*ecorr3_nucl&
1711             +wscbase*escbase&
1712             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1713             + wcatnucl*ecation_nucl
1714       endif
1715 #else
1716       if (shield_mode.gt.0) then
1717         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1718           +ft(1)*welec*(ees+evdw1) &
1719           +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1720 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1721           +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1722           +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1723           +ft(2)*wturn3*eello_turn3 &
1724           +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1725           +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1726             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1727             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1728             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1729             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1730             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1731             +wcorr3_nucl*ecorr3_nucl&
1732             +wscbase*escbase&
1733             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1734             + wcatnucl*ecation_nucl
1735       else
1736         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
1737           +ft(1)*welec*(ees+evdw1) &
1738           +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1739 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1740           +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1741           +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1742           +ft(2)*wturn3*eello_turn3 &
1743           +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
1744           +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1745             +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1746             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1747             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1748             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1749             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1750             +wcorr3_nucl*ecorr3_nucl&
1751             +wscbase*escbase&
1752             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1753             + wcatnucl*ecation_nucl
1754       endif
1755 #endif
1756
1757             weight=-betaT*(etot-potEmin)+entfac(t)
1758 !            write (iout,*) "k",k," t",t," weight",weight
1759             if (weight.gt.weimax(k,iparm)) weimax(k,iparm)=weight
1760           enddo
1761 #ifdef MPI
1762         enddo
1763 #ifdef DEBUG
1764         write (iout,*) "weimax before REDUCE"
1765         write (iout,*) (weimax(k,iparm),k=0,ngridt)
1766 #endif
1767         do k=0,nGridT
1768           weimax_(k)=weimax(k,iparm)
1769         enddo
1770         call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1, &
1771         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
1772 #ifdef DEBUG
1773         write (iout,*) "weimax"
1774         write (iout,*) (weimax(k,iparm),k=0,ngridt)
1775 #endif
1776         do k=0,nGridT
1777           weimax_(k)=weimax(k,iparm)
1778         enddo
1779         call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1, &
1780         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
1781 #ifdef DEBUG
1782         write (iout,*) "weimax"
1783         write (iout,*) (weimax(k,iparm),k=0,ngridt)
1784 #endif
1785         do k=0,nGridT
1786           temper=startGridT+k*delta_T
1787           betaT=1.0d0/(1.987D-3*temper)
1788 !          call temp_scalfac(temper,ft,ftprim,ftbis,*10)
1789       if (rescale_mode.eq.1) then
1790         quot=temper/T0
1791         quotl=1.0d0
1792         kfacl=1.0d0
1793         do l=1,5
1794           quotl1=quotl
1795           quotl=quotl*quot
1796           kfacl=kfacl*kfac
1797           denom=kfacl-1.0d0+quotl
1798           fT(l)=kfacl/denom
1799           ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
1800           ftbis(l)=l*kfacl*quotl1*&
1801           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
1802         enddo
1803 #if defined(FUNCTH)
1804         ft(6)=(320.0d0+80.0d0*dtanh((temper-320.0d0)/80.0d0))/ &
1805                  320.0d0
1806         ftprim(6)=1.0d0/(320.0d0*dcosh((temper-320.0d0)/80.0d0)**2)
1807         ftbis(6)=-2.0d0*dtanh((temper-320.0d0)/80.0d0) &
1808                /(320.0d0*80.0d0*dcosh((temper-320.0d0)/80.0d0)**3)
1809 #elif defined(FUNCT)
1810         fT(6)=temper/T0
1811         ftprim(6)=1.0d0/T0
1812         ftbis(6)=0.0d0
1813 #else
1814         fT(6)=1.0d0
1815         ftprim(6)=0.0d0
1816         ftbis(6)=0.0d0
1817 #endif
1818       else if (rescale_mode.eq.2) then
1819         quot=temper/T0
1820         quotl=1.0d0
1821         do l=1,5
1822           quotl1=quotl
1823           quotl=quotl*quot
1824           eplus=dexp(quotl)
1825           eminus=dexp(-quotl)
1826           logfac=1.0d0/dlog(eplus+eminus)
1827           tanhT=(eplus-eminus)/(eplus+eminus)
1828           fT(l)=1.12692801104297249644d0*logfac
1829           ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
1830           ftbis(l)=(l-1)*ftprim(l)/(quot*T0)- &
1831          2*l*quotl1/T0*logfac* &
1832          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2) &
1833          +ftprim(l)*tanhT)
1834         enddo
1835 #if defined(FUNCTH)
1836         ft(6)=(320.0d0+80.0d0*dtanh((temper-320.0d0)/80.0d0))/ &
1837                 320.0d0
1838         ftprim(6)=1.0d0/(320.0d0*dcosh((temper-320.0d0)/80.0d0)**2)
1839         ftbis(6)=-2.0d0*dtanh((temper-320.0d0)/80.0d0) &
1840               /(320.0d0*80.0d0*dcosh((temper-320.0d0)/80.0d0)**3)
1841 #elif defined(FUNCT)
1842         fT(6)=temper/T0
1843         ftprim(6)=1.0d0/T0
1844         ftbis(6)=0.0d0
1845 #else
1846         fT(6)=1.0d0
1847         ftprim(6)=0.0d0
1848         ftbis(6)=0.0d0
1849 #endif
1850       else if (rescale_mode.eq.0) then
1851         do l=1,5
1852           fT(l)=1.0d0
1853           ftprim(l)=0.0d0
1854         enddo
1855       else
1856         write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
1857          rescale_mode
1858         call flush(iout)
1859 !        return1
1860       endif
1861
1862
1863
1864           do t=1,scount(me1)
1865 #else
1866           do t=1,ntot(islice)
1867 #endif
1868             ind = ind_point(t)
1869 #ifdef MPI
1870             hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
1871 #else
1872             hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
1873 #endif
1874 !          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
1875 !          call restore_parm(iparm)
1876 !            call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
1877       evdw=enetb(1,t,iparm)
1878       evdw_t=enetb(20,t,iparm)
1879 #ifdef SCP14
1880       evdw2_14=enetb(18,t,iparm)
1881       evdw2=enetb(2,t,iparm)+evdw2_14
1882 #else
1883       evdw2=enetb(2,t,iparm)
1884       evdw2_14=0.0d0
1885 #endif
1886 #ifdef SPLITELE
1887       ees=enetb(3,t,iparm)
1888       evdw1=enetb(16,t,iparm)
1889 #else
1890       ees=enetb(3,t,iparm)
1891       evdw1=0.0d0
1892 #endif
1893       ecorr=enetb(4,t,iparm)
1894       ecorr5=enetb(5,t,iparm)
1895       ecorr6=enetb(6,t,iparm)
1896       eel_loc=enetb(7,t,iparm)
1897       eello_turn3=enetb(8,t,iparm)
1898       eello_turn4=enetb(9,t,iparm)
1899       eello_turn6=enetb(10,t,iparm)
1900       ebe=enetb(11,t,iparm)
1901       escloc=enetb(12,t,iparm)
1902       etors=enetb(13,t,iparm)
1903       etors_d=enetb(14,t,iparm)
1904       ehpb=enetb(15,t,iparm)
1905       estr=enetb(17,t,iparm)
1906       esccor=enetb(21,t,iparm)
1907       edihcnstr=enetb(19,t,iparm)
1908 !      eliptran=enetb(22,t,iparm)
1909 !      esaxs=enetb(26,t,iparm)
1910 !      ehomology_constr=enetb(27,t,iparm)
1911 !      edfadis=enetb(28,t,iparm)
1912 !      edfator=enetb(29,t,iparm)
1913 !      edfanei=enetb(30,t,iparm)
1914 !      edfabet=enetb(31,t,iparm)
1915           ecationcation=enetb(41,i,iparm)
1916           ecation_prot=enetb(42,i,iparm)
1917             evdwpp =      enetb(26,i,iparm)
1918             eespp =      enetb(27,i,iparm)
1919             evdwpsb =      enetb(28,i,iparm)
1920             eelpsb =      enetb(29,i,iparm)
1921             evdwsb =      enetb(30,i,iparm)
1922             eelsb =      enetb(31,i,iparm)
1923             estr_nucl =      enetb(32,i,iparm)
1924             ebe_nucl =      enetb(33,i,iparm)
1925             esbloc =      enetb(34,i,iparm)
1926             etors_nucl =      enetb(35,i,iparm)
1927             etors_d_nucl =      enetb(36,i,iparm)
1928             ecorr_nucl =      enetb(37,i,iparm)
1929             ecorr3_nucl =      enetb(38,i,iparm)
1930             epeppho=   enetb(49,i,iparm)
1931             escpho=    enetb(48,i,iparm)
1932             epepbase=  enetb(47,i,iparm)
1933             escbase=   enetb(46,i,iparm)
1934             ecation_nucl= enetb(50,i,iparm)
1935             ehomology_constr=enetb(51,i,iparm)
1936
1937 #ifdef SPLITELE
1938       if (shield_mode.gt.0) then
1939         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
1940          +ft(1)*welec*ees &
1941          +ft(1)*wvdwpp*evdw1 &
1942          +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
1943 !     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
1944          +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
1945          +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
1946          +ft(2)*wturn3*eello_turn3 &
1947          +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
1948          +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
1949          +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
1950          +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
1951          +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
1952          +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
1953          *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
1954          +wcorr3_nucl*ecorr3_nucl&
1955          +wscbase*escbase&
1956          +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
1957          + wcatnucl*ecation_nucl
1958
1959 !        eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
1960 !!     &            +ftprim(6)*evdw_t
1961 !     &    +ftprim(1)*wscp*evdw2
1962 !     &    +ftprim(1)*welec*ees
1963 !     &    +ftprim(1)*wvdwpp*evdw1
1964 !     &    +ftprim(1)*wtor*etors+
1965 !     &    ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
1966 !     &    ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
1967 !     &    ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+
1968 !     &    ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
1969 !     &    ftprim(1)*wsccor*esccor
1970
1971 !        ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
1972 !     &    +ftbis(1)*wscp*evdw2+
1973 !     &    ftbis(1)*welec*ees
1974 !     &    +ftbis(1)*wvdwpp*evdw
1975 !     &    +ftbis(1)*wtor*etors+
1976 !     &    ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
1977 !     &    ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
1978 !     &    ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+
1979 !     &    ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
1980 !     &    ftbis(1)*wsccor*esccor
1981
1982             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
1983                 +ftprim(1)*wtor*etors+ &
1984                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
1985                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
1986                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
1987                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
1988                  ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
1989                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
1990             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
1991                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
1992                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
1993                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
1994                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
1995                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
1996                  +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
1997
1998
1999
2000       else
2001         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
2002          +wvdwpp*evdw1 &
2003          +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
2004 !c    &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
2005          +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
2006          +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
2007          +ft(2)*wturn3*eello_turn3 &
2008          +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
2009          +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
2010          +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
2011             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
2012             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
2013             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
2014             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
2015             +wcorr3_nucl*ecorr3_nucl&
2016             +wscbase*escbase&
2017             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
2018             + wcatnucl*ecation_nucl
2019
2020
2021             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
2022                 +ftprim(1)*wtor*etors+ &
2023                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
2024                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
2025                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
2026                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
2027                  ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
2028                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
2029             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
2030                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
2031                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
2032                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
2033                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
2034                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
2035                  +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
2036
2037
2038
2039 !       eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
2040 !    &    +ftprim(1)*wtor*etors+
2041 !    &    ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
2042 !    &    ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
2043 !    &    ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+
2044 !    &    ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
2045 !    &    ftprim(1)*wsccor*esccor
2046 !       ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
2047 !    &    ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
2048 !    &    ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
2049 !    &    ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+
2050 !    &    ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
2051 !    &    ftbis(1)*wsccor*esccor
2052       endif
2053 #else
2054       if (shield_mode.gt.0) then
2055         etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 &
2056          +ft(1)*welec*(ees+evdw1) &
2057          +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
2058 !c    &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
2059          +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
2060          +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
2061          +ft(2)*wturn3*eello_turn3 &
2062          +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
2063          +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
2064          +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
2065             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
2066             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
2067             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
2068             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
2069             +wcorr3_nucl*ecorr3_nucl&
2070             +wscbase*escbase&
2071             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
2072             + wcatnucl*ecation_nucl
2073
2074
2075             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
2076                 +ftprim(1)*wtor*etors+ &
2077                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
2078                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
2079                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
2080                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
2081                  ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
2082                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
2083             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
2084                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
2085                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
2086                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
2087                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
2088                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
2089                  +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
2090
2091
2092
2093
2094 !       eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
2095 !    &    +ftprim(1)*welec*(ees+evdw1)
2096 !    &    +ftprim(1)*wtor*etors+
2097 !    &     ftprim(1)*wscp*evdw2+
2098 !    &     ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
2099 !    &     ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
2100 !    &     ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+
2101 !    &     ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
2102 !    &     ftprim(1)*wsccor*esccor
2103 !      ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
2104 !    &    +ftbis(1)*wscp*evdw2
2105 !    &    +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
2106 !    &     ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
2107 !    &     ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
2108 !    &     ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+
2109 !    &     ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
2110 !    &     ftprim(1)*wsccor*esccor
2111       else
2112         etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 &
2113          +ft(1)*welec*(ees+evdw1) &
2114          +wang*ebe+ft(1)*wtor*etors+wscloc*escloc &
2115 !c    &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
2116          +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 &
2117          +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 &
2118          +ft(2)*wturn3*eello_turn3 &
2119          +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
2120          +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
2121          +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
2122             +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
2123             +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
2124             +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
2125             *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
2126             +wcorr3_nucl*ecorr3_nucl&
2127             +wscbase*escbase&
2128             +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
2129             + wcatnucl*ecation_nucl
2130
2131
2132             eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) &
2133                 +ftprim(1)*wtor*etors+ &
2134                  ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ &
2135                  ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+ &
2136                  ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+ &
2137                  ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
2138                  ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
2139                  +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
2140             ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+ &
2141                  ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+ &
2142                  ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+ &
2143                  ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+ &
2144                  ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
2145                  ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
2146                  +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
2147
2148
2149
2150
2151
2152 !       eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
2153 !    &    +ftprim(1)*wtor*etors+
2154 !    &     ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
2155 !    &     ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
2156 !    &     ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eello_turn6+
2157 !    &     ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
2158 !    &     ftprim(1)*wsccor*esccor
2159 !       ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
2160 !    &     ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
2161 !    &     ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
2162 !    &     ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eello_turn6+
2163 !    &     ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
2164 !    &     ftprim(1)*wsccor*esccor
2165       endif
2166 #endif
2167
2168             weight=dexp(-betaT*(etot-potEmin)+entfac(t)-weimax(k,iparm))
2169 #ifdef DEBUG
2170             write (iout,*) "iparm",iparm," t",t," betaT",betaT, &
2171              " etot",etot," entfac",entfac(t)," boltz", &
2172              -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm), &
2173              " weight",weight," ebis",ebis
2174 #endif
2175             etot=etot-temper*eprim
2176 #ifdef MPI
2177             sumW_p(k,iparm)=sumW_p(k,iparm)+weight
2178             sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
2179             sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
2180             sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
2181             do j=1,nQ+2
2182               sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
2183               sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
2184               sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
2185               +etot*q(j,t)*weight
2186             enddo
2187 #else
2188             sumW(k,iparm)=sumW(k,iparm)+weight
2189             sumE(k,iparm)=sumE(k,iparm)+etot*weight
2190             sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
2191             sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
2192             do j=1,nQ+2
2193               sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
2194               sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
2195               sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
2196               +etot*q(j,t)*weight
2197             enddo
2198 #endif
2199           enddo ! t
2200         enddo ! k
2201
2202         do ib=1,nT_h(iparm)
2203           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,&
2204           MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2205           if (rmsrgymap) then
2206           call MPI_Reduce(hrmsrgy_p(0,0,ib),hrmsrgy(0,0,ib),&
2207          (MaxBinRgy+1)*(nbin_rms+1),MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
2208              WHAM_COMM,IERROR)
2209           endif
2210         enddo
2211         call MPI_Reduce(upindE_p,upindE,1,&
2212           MPI_INTEGER,MPI_MAX,Master,WHAM_COMM,IERROR)
2213         call MPI_Reduce(histE_p(0),histE(0),maxindE,&
2214           MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2215
2216         if (me1.eq.master) then
2217
2218         if (histout) then
2219
2220         write (iout,'(6x,$)')
2221         write (iout,'(f20.2,$)') (1.0d0/(1.987D-3*beta_h(ib,iparm)),&
2222          ib=1,nT_h(iparm))
2223         write (iout,*)
2224
2225         write (iout,'(/a)') 'Final histograms'
2226         if (histfile) then
2227           if (nslice.eq.1) then
2228             if (separate_parset) then
2229               write(licz3,"(bz,i3.3)") myparm
2230               histname=prefix(:ilen(prefix))//'_par'//licz3//'.hist'
2231             else
2232               histname=prefix(:ilen(prefix))//'.hist'
2233             endif
2234           else
2235             if (separate_parset) then
2236               write(licz3,"(bz,i3.3)") myparm
2237               histname=prefix(:ilen(prefix))//'_par'//licz3// &
2238                '_slice_'//licz2//'.hist'
2239             else
2240               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.hist'
2241             endif
2242           endif
2243 #if defined(AIX) || defined(PGI)
2244           open (ihist,file=histname,position='append')
2245 #else
2246           open (ihist,file=histname,access='append')
2247 #endif
2248         endif
2249
2250         do t=0,tmax
2251           liczbaW=t
2252           sumH=0.0d0
2253           do ib=1,nT_h(iparm)
2254             sumH=sumH+hfin(t,ib)
2255           enddo
2256           if (sumH.gt.0.0d0) then
2257             do j=1,nQ
2258               jj = mod(liczbaW,nbin1)
2259               liczbaW=liczbaW/nbin1
2260               write (iout,'(f6.3,$)') dmin+(jj+0.5d0)*delta
2261               if (histfile) &
2262                  write (ihist,'(f6.3,$)') dmin+(jj+0.5d0)*delta
2263             enddo
2264             do ib=1,nT_h(iparm)
2265               write (iout,'(e20.10,$)') hfin(t,ib)
2266               if (histfile) write (ihist,'(e20.10,$)') hfin(t,ib)
2267             enddo
2268             write (iout,'(i5)') iparm
2269             if (histfile) write (ihist,'(i5)') iparm
2270           endif
2271         enddo
2272
2273         endif
2274
2275         if (entfile) then
2276           if (nslice.eq.1) then
2277             if (separate_parset) then
2278               write(licz3,"(bz,i3.3)") myparm
2279               histname=prefix(:ilen(prefix))//"_par"//licz3//'.ent'
2280             else
2281               histname=prefix(:ilen(prefix))//'.ent'
2282             endif
2283           else
2284             if (separate_parset) then
2285               write(licz3,"(bz,i3.3)") myparm
2286               histname=prefix(:ilen(prefix))//'par_'//licz3// &
2287                  '_slice_'//licz2//'.ent'
2288             else
2289               histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.ent'
2290             endif
2291           endif
2292 #if defined(AIX) || defined(PGI)
2293           open (ihist,file=histname,position='append')
2294 #else
2295           open (ihist,file=histname,access='append')
2296 #endif
2297           write (ihist,'(a)') "# Microcanonical entropy"
2298           do i=0,upindE
2299             write (ihist,'(f8.0,$)') dint(potEmin)+i
2300             if (histE(i).gt.0.0e0) then
2301               write (ihist,'(f15.5,$)') dlog(histE(i))
2302             else
2303               write (ihist,'(f15.5,$)') 0.0d0
2304             endif
2305           enddo
2306           write (ihist,*)
2307           close(ihist)
2308         endif
2309         write (iout,*) "Microcanonical entropy"
2310         do i=0,upindE
2311           write (iout,'(f8.0,$)') dint(potEmin)+i
2312           if (histE(i).gt.0.0e0) then
2313             write (iout,'(f15.5,$)') dlog(histE(i))
2314           else
2315             write (iout,'(f15.5,$)') 0.0d0
2316           endif
2317           write (iout,*)
2318         enddo
2319         if (rmsrgymap) then
2320           if (nslice.eq.1) then
2321             if (separate_parset) then
2322               write(licz3,"(bz,i3.3)") myparm
2323               histname=prefix(:ilen(prefix))//'_par'//licz3//'.rmsrgy'
2324             else
2325               histname=prefix(:ilen(prefix))//'.rmsrgy'
2326             endif
2327           else
2328             if (separate_parset) then
2329               write(licz3,"(bz,i3.3)") myparm
2330               histname=prefix(:ilen(prefix))//'_par'//licz3// &
2331                '_slice_'//licz2//'.rmsrgy'
2332             else
2333              histname=prefix(:ilen(prefix))//'_slice_'//licz2//'.rmsrgy'
2334             endif
2335           endif
2336 #if defined(AIX) || defined(PGI)
2337           open (ihist,file=histname,position='append')
2338 #else
2339           open (ihist,file=histname,access='append')
2340 #endif
2341           do i=0,nbin_rms
2342             do j=0,nbin_rgy
2343               write(ihist,'(2f8.2,$)') &
2344                 rgymin+deltrgy*j,rmsmin+deltrms*i
2345               do ib=1,nT_h(iparm)
2346                 if (hrmsrgy(j,i,ib).gt.0.0d0) then
2347                   write(ihist,'(e14.5,$)') &
2348                     -dlog(hrmsrgy(j,i,ib))/beta_h(ib,iparm) &
2349                     +potEmin
2350                 else
2351                   write(ihist,'(e14.5,$)') 1.0d6
2352                 endif
2353               enddo
2354               write (ihist,'(i2)') iparm
2355             enddo
2356           enddo
2357           close(ihist)
2358         endif
2359         endif
2360       enddo ! iparm
2361 #ifdef MPI
2362       call MPI_Reduce(hfin_ent_p(0),hfin_ent(0),nbin,&
2363          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2364       call MPI_Reduce(sumW_p(0,1),sumW(0,1),(nGridT+1)*nParmSet,&
2365          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2366       call MPI_Reduce(sumE_p(0,1),sumE(0,1),(nGridT+1)*nParmSet,&
2367          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2368       call MPI_Reduce(sumEbis_p(0,1),sumEbis(0,1),(nGridT+1)*nParmSet,&
2369          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2370       call MPI_Reduce(sumEsq_p(0,1),sumEsq(0,1),(nGridT+1)*nParmSet,&
2371          MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
2372       call MPI_Reduce(sumQ_p(1,0,1),sumQ(1,0,1),&
2373          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
2374          WHAM_COMM,IERROR)
2375       call MPI_Reduce(sumQsq_p(1,0,1),sumQsq(1,0,1),&
2376          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
2377          WHAM_COMM,IERROR)
2378       call MPI_Reduce(sumEQ_p(1,0,1),sumEQ(1,0,1),&
2379          MaxQ1*(nGridT+1)*nParmSet,MPI_DOUBLE_PRECISION,MPI_SUM,Master,&
2380          WHAM_COMM,IERROR)
2381       if (me.eq.master) then
2382 #endif
2383       write (iout,'(/a)') 'Thermal characteristics of folding'
2384       if (nslice.eq.1) then
2385         nazwa=prefix
2386       else
2387         nazwa=prefix(:ilen(prefix))//"_slice_"//licz2
2388       endif
2389       iln=ilen(nazwa)
2390       if (nparmset.eq.1 .and. .not.separate_parset) then
2391         nazwa=nazwa(:iln)//".thermal"
2392       else if (nparmset.eq.1 .and. separate_parset) then
2393         write(licz3,"(bz,i3.3)") myparm
2394         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
2395       endif
2396       do iparm=1,nParmSet
2397       if (nparmset.gt.1) then
2398         write(licz3,"(bz,i3.3)") iparm
2399         nazwa=nazwa(:iln)//"_par_"//licz3//".thermal"
2400       endif
2401       open(34,file=nazwa)
2402       if (separate_parset) then
2403         write (iout,'(a,i3)') "Parameter set",myparm
2404       else
2405         write (iout,'(a,i3)') "Parameter set",iparm
2406       endif
2407       do i=0,NGridT
2408         betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
2409         if (betaT.ge.beta_h(1,iparm)) then
2410           potEmin=potEmin_all(1,iparm)
2411         else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
2412           potEmin=potEmin_all(nT_h(iparm),iparm)
2413         else
2414           do l=1,nT_h(iparm)-1
2415             if (betaT.le.beta_h(l,iparm) .and. &
2416                betaT.gt.beta_h(l+1,iparm)) then
2417               potEmin=potEmin_all(l,iparm)
2418               exit
2419             endif
2420           enddo
2421         endif
2422
2423         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
2424         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ &
2425           sumW(i,iparm)
2426         sumEsq(i,iparm)=(sumEsq(i,iparm)/sumW(i,iparm) &
2427           -sumE(i,iparm)**2)/(1.987D-3*(startGridT+i*delta_T)**2)
2428         do j=1,nQ+2
2429           sumQ(j,i,iparm)=sumQ(j,i,iparm)/sumW(i,iparm)
2430           sumQsq(j,i,iparm)=sumQsq(j,i,iparm)/sumW(i,iparm) &
2431            -sumQ(j,i,iparm)**2
2432           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm) &
2433            -sumQ(j,i,iparm)*sumE(i,iparm)
2434         enddo
2435         sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3* &
2436          (startGridT+i*delta_T))+potEmin
2437         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
2438          sumW(i,iparm),sumE(i,iparm)
2439         write (iout,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
2440         write (iout,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
2441          (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
2442         write (iout,*)
2443         write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T,&
2444          sumW(i,iparm),sumE(i,iparm)
2445         write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2)
2446         write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm),&
2447          (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2)
2448         write (34,*)
2449         call flush(34)
2450       enddo
2451       close(34)
2452       enddo
2453       if (histout) then
2454       do t=0,tmax
2455         if (hfin_ent(t).gt.0.0d0) then
2456           liczbaW=t
2457           jj = mod(liczbaW,nbin1)
2458           write (iout,'(f6.3,e20.10," ent")') dmin+(jj+0.5d0)*delta,&
2459            hfin_ent(t)
2460           if (histfile) write (ihist,'(f6.3,e20.10," ent")') &
2461             dmin+(jj+0.5d0)*delta,&
2462            hfin_ent(t)
2463         endif
2464       enddo
2465       if (histfile) close(ihist)
2466       endif
2467
2468 #ifdef ZSCORE
2469 ! Write data for zscore
2470       if (nslice.eq.1) then
2471         zscname=prefix(:ilen(prefix))//".zsc"
2472       else
2473         zscname=prefix(:ilen(prefix))//"_slice_"//licz2//".zsc"
2474       endif
2475 #if defined(AIX) || defined(PGI)
2476       open (izsc,file=prefix(:ilen(prefix))//'.zsc',position='append')
2477 #else
2478       open (izsc,file=prefix(:ilen(prefix))//'.zsc',access='append')
2479 #endif
2480       write (izsc,'("NQ=",i1," NPARM=",i1)') nQ,nParmSet
2481       do iparm=1,nParmSet
2482         write (izsc,'("NT=",i1)') nT_h(iparm)
2483       do ib=1,nT_h(iparm)
2484         write (izsc,'("TEMP=",f6.1," NR=",i2," SNK=",$)') &
2485           1.0d0/(beta_h(ib,iparm)*1.987D-3),nR(ib,iparm)
2486         jj = min0(nR(ib,iparm),7)
2487         write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=1,jj)
2488         write (izsc,'(a1,$)') (" ",i=22+8*jj+1,79)
2489         write (izsc,'("&")')
2490         if (nR(ib,iparm).gt.7) then
2491           do ii=8,nR(ib,iparm),9
2492             jj = min0(nR(ib,iparm),ii+8)
2493             write (izsc,'(i8,$)') (snk(i,ib,iparm,islice),i=ii,jj)
2494             write (izsc,'(a1,$') (" ",i=(jj-ii+1)*8+1,79)
2495             write (izsc,'("&")')
2496           enddo
2497         endif
2498         write (izsc,'("FI=",$)')
2499         jj=min0(nR(ib,iparm),7)
2500         write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=1,jj)
2501         write (izsc,'(a1,$)') (" ",i=3+10*jj+1,79)
2502         write (izsc,'("&")')
2503         if (nR(ib,iparm).gt.7) then
2504           do ii=8,nR(ib,iparm),9
2505             jj = min0(nR(ib,iparm),ii+8)
2506             write (izsc,'(f10.5,$)') (fi(i,ib,iparm),i=ii,jj)
2507             if (jj.eq.nR(ib,iparm)) then
2508               write (izsc,*)
2509             else
2510               write (izsc,'(a1,$)') (" ",i=10*(jj-ii+1)+1,79)
2511               write (izsc,'(t80,"&")')
2512             endif
2513           enddo
2514         endif
2515         do i=1,nR(ib,iparm)
2516           write (izsc,'("KH=",$)')
2517           write (izsc,'(f7.2,$)') (Kh(j,i,ib,iparm),j=1,nQ)
2518           write (izsc,'(" Q0=",$)')
2519           write (izsc,'(f7.5,$)') (q0(j,i,ib,iparm),j=1,nQ)
2520           write (izsc,*)
2521         enddo
2522       enddo
2523       enddo
2524       close(izsc)
2525 #endif
2526 #ifdef MPI
2527       endif
2528 #endif
2529       return
2530       end subroutine WHAMCALC
2531 !-----------------------------------------------------------------------------
2532       end module wham_calc
2533