source/wham/src-HCD-5D/wham_calc1.F

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