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