update
[unres.git] / source / wham / src-M / wham_calc1.F
index 15d6716..35a9278 100644 (file)
@@ -18,9 +18,8 @@
       parameter (NGridT=400)
       integer MaxBinRms,MaxBinRgy
       parameter (MaxBinRms=100,MaxBinRgy=100)
-      integer MaxHdim
-c      parameter (MaxHdim=200000)
-      parameter (MaxHdim=200)
+c      integer MaxHdim
+c      parameter (MaxHdim=200)
       integer maxinde
       parameter (maxinde=200)
 #ifdef MPI
@@ -36,6 +35,7 @@ c      parameter (MaxHdim=200000)
       include "COMMON.SBRIDGE"
       include "COMMON.PROT"
       include "COMMON.ENEPS"
+      include "COMMON.SHIELD"
       integer MaxPoint,MaxPointProc
       parameter (MaxPoint=MaxStr,
      & MaxPointProc=MaxStr_Proc)
@@ -50,7 +50,8 @@ c      parameter (MaxHdim=200000)
       double precision energia(0:max_ene)
 #ifdef MPI
       integer tmax_t,upindE_p
-      double precision fi_p(MaxR,MaxT_h,Max_Parm)
+      double precision fi_p(MaxR,MaxT_h,Max_Parm),
+     & fimax_p(MaxR,MaxT_h,Max_Parm)
       double precision sumW_p(0:nGridT,Max_Parm),
      & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
      & sumQ_p(MaxQ1,0:nGridT,Max_Parm),
@@ -63,6 +64,8 @@ c      parameter (MaxHdim=200000)
      & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
       double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
       double precision potEmin_t,entmin_p,entmax_p
+      double precision ePMF,ePMF_q
+      double precision weimax_(0:ngridT)
       integer histent_p(0:2000)
       logical lprint /.true./
 #endif
@@ -74,17 +77,19 @@ c      parameter (MaxHdim=200000)
      & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
      & weight,econstr
       double precision fi(MaxR,maxT_h,Max_Parm),
+     & fimax(MaxR,maxT_h,Max_Parm),
      & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
      & hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
      & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
      & potEmin,ent,
-     & hfin_ent(0:MaxHdim),vmax,aux
+     & hfin_ent(0:MaxHdim),vmax,aux,weimax(0:nGridT,Max_Parm)
       double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
      &  eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
      &  eplus,eminus,logfac,tanhT,tt
       double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
      &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
-     &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor
+     &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
+     &  eliptran
 
       integer ind_point(maxpoint),upindE,indE
       character*16 plik
@@ -112,11 +117,13 @@ c      parameter (MaxHdim=200000)
       do t=0,MaxN
         htot(t)=0
       enddo
+C#define DEBUG
 #ifdef MPI
       do i=1,scount(me1)
 #else
       do i=1,ntot(islice)
 #endif
+c        write (iout,*) "i",i," potE",(potE(i,j),j=1,nParmset)
         do j=1,nParmSet
           if (potE(i,j).le.potEmin) potEmin=potE(i,j)
         enddo
@@ -134,7 +141,7 @@ c      parameter (MaxHdim=200000)
           endif
 c          write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
 c     &      ind_point(i)
-          call flush(iout)
+c          call flush(iout)
           if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
             write (iout,*) "Error - index exceeds range for point",i,
      &      " q=",q(j,i)," ind",ind_point(i)
@@ -154,8 +161,8 @@ c     &      ind_point(i)
         call flush(iout)
 #endif
       enddo ! i
-      call flush(iout)
 
+      write (iout,*) "potEmin before reduce",potEmin
       nbin=nbin1**nQ-1
       write (iout,'(a)') "Numbers of counts in Q bins"
       do t=0,tmax
@@ -194,7 +201,8 @@ c     &      ind_point(i)
      &  MPI_MIN,WHAM_COMM,IERROR)
       call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
      &  MPI_MAX,WHAM_COMM,IERROR)
-      potEmin=potEmin_t/2
+c      potEmin=potEmin_t/2
+      potEmin=potEmin_t
       rgymin=rgymin_t
       rgymax=rgymax_t
       rmsmin=rmsmin_t
@@ -227,7 +235,7 @@ c        write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
         do iparm=1,nParmSet
 #ifdef DEBUG
           write (iout,'(2i5,21f8.2)') i,iparm,
-     &     (enetb(k,i,iparm),k=1,21)
+     &     (enetb(k,i,iparm),k=1,22)
 #endif
           call restore_parm(iparm)
 #ifdef DEBUG
@@ -312,6 +320,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             estr=enetb(18,i,iparm)
             esccor=enetb(19,i,iparm)
             edihcnstr=enetb(20,i,iparm)
+            eliptran=enetb(22,i,iparm)
+
 #ifdef DEBUG
             write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
      &       evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
@@ -319,30 +329,61 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 #endif
 
 #ifdef SPLITELE
+            if (shield_mode.gt.0) then
+            etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+     &      +ft(1)*welec*ees
+     &      +ft(1)*wvdwpp*evdw1
+     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+!     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &      +ft(2)*wturn3*eello_turn3
+     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+     &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &      +wbond*estr+wliptran*eliptran
+             else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
      &      +wvdwpp*evdw1
      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+c     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
      &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr
+     &      +wbond*estr+wliptran*eliptran
+             endif
 #else
+      if (shield_mode.gt.0) then
+            etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+     &      +ft(1)*welec*(ees+evdw1)
+     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &      +ft(2)*wturn3*eello_turn3
+     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+     &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &      +wbond*estr+wliptran*eliptran
+            else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
      &      +ft(1)*welec*(ees+evdw1)
      &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+c     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &      +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
      &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
      &      +ft(2)*wturn3*eello_turn3
      &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
      &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr
+     &      +wbond*estr+wliptran*eliptran
+           endif
+
 #endif
 #ifdef DEBUG
             write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
      &        etot,potEmin
 #endif
+#define DEBUG
 #ifdef DEBUG
             if (iparm.eq.1 .and. ib.eq.1) then
             write (iout,*)"Conformation",i
@@ -353,6 +394,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             call enerprint(energia(0),fT)
             endif
 #endif
+#undef DEBUG
             do kk=1,nR(ib,iparm)
               Econstr=0.0d0
               do j=1,nQ
@@ -360,6 +402,11 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
                 Econstr=Econstr+Kh(j,kk,ib,iparm)
      &           *(dd-q0(j,kk,ib,iparm))**2
               enddo
+c Adaptive potential contribution
+              if (adaptive) then
+                call PMF_energy(q(1,i),ib,kk,iparm,ePMF,ePMF_q)
+                Econstr=Econstr+ePMF
+              endif
               v(i,kk,ib,iparm)=
      &          -beta_h(ib,iparm)*(etot-potEmin+Econstr)
 #ifdef DEBUG
@@ -406,6 +453,31 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
           write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
 #endif
         enddo
+
+        do iparm=1,nParmSet
+          do iib=1,nT_h(iparm)
+            do ii=1,nR(iib,iparm)
+#ifdef MPI
+              fimax_p(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
+              do t=2,scount(me)
+                if(v(t,ii,iib,iparm)+entfac(t).gt.fimax_p(ii,iib,iparm))
+     &            fimax_p(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
+              enddo
+#else
+              fimax(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
+              do t=2,ntot(islice)
+                if(v(t,ii,iib,iparm)+entfac(t).gt.fimax(ii,iib,iparm))
+     &            fimax(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
+              enddo
+#endif
+            enddo ! ii
+          enddo ! iib
+        enddo ! iparm
+#ifdef MPI
+        call MPI_AllReduce(fimax_p(1,1,1),fimax(1,1,1),
+     &   maxR*MaxT_h*nParmSet,MPI_DOUBLE_PRECISION,
+     &   MPI_MAX,WHAM_COMM,IERROR)
+#endif
         do iparm=1,nParmSet
           do iib=1,nT_h(iparm)
             do ii=1,nR(iib,iparm)
@@ -413,7 +485,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
               fi_p(ii,iib,iparm)=0.0d0
               do t=1,scount(me)
                 fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
-     &            +dexp(v(t,ii,iib,iparm)+entfac(t))
+     &            +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
 #ifdef DEBUG
               write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
      &         v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
@@ -423,7 +495,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
               fi(ii,iib,iparm)=0.0d0
               do t=1,ntot(islice)
                 fi(ii,iib,iparm)=fi(ii,iib,iparm)
-     &            +dexp(v(t,ii,iib,iparm)+entfac(t))
+     &            +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
               enddo
 #endif
             enddo ! ii
@@ -441,10 +513,10 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
           enddo
         enddo
 #endif
-        write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
-     &   maxR*MaxT_h*nParmSet
-        write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
-     &   " WHAM_COMM",WHAM_COMM
+c        write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
+c     &   maxR*MaxT_h*nParmSet
+c        write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
+c     &   " WHAM_COMM",WHAM_COMM
         call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
      &   MPI_DOUBLE_PRECISION,
      &   MPI_SUM,Master,WHAM_COMM,IERROR)
@@ -464,7 +536,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
         do iparm=1,nParmSet
           do ib=1,nT_h(iparm)
             do i=1,nR(ib,iparm)
-              fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
+              fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fimax(i,ib,iparm)
               avefi=avefi+fi(i,ib,iparm)
             enddo
           enddo
@@ -607,12 +679,12 @@ c        ient=-dlog(entfac(t))-entmin
         write(iout,'(2f15.4)') entmin+i,histent(i)
       enddo
 #endif
-      
-#ifdef MPI
-c      write (iout,*) "me1",me1," scount",scount(me1)
-
       do iparm=1,nParmSet
 
+        call restore_parm(iparm)
+c
+C Histograms
+c
 #ifdef MPI
         do ib=1,nT_h(iparm)
           do t=0,tmax
@@ -642,182 +714,102 @@ c      write (iout,*) "me1",me1," scount",scount(me1)
             enddo
           enddo
         enddo
-
+#ifdef MPI
         do t=1,scount(me1)
 #else
         do t=1,ntot(islice)
 #endif
-          ind = ind_point(t)
+          indE = aint(potE(t,iparm)-aint(potEmin))
+          if (indE.ge.0 .and. indE.le.maxinde) then
+            if (indE.gt.upindE_p) upindE_p=indE
+            histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
+          endif
 #ifdef MPI
-          hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
+          do ib=1,nT_h(iparm)
+            expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+            hfin_p(ind,ib)=hfin_p(ind,ib)+
+     &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+             if (rmsrgymap) then
+               indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy) 
+               indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
+               hrmsrgy_p(indrgy,indrms,ib)=
+     &           hrmsrgy_p(indrgy,indrms,ib)+expfac
+             endif
+          enddo
 #else
-          hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
+          do ib=1,nT_h(iparm)
+            expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+            hfin(ind,ib)=hfin(ind,ib)+
+     &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+             if (rmsrgymap) then
+               indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
+               indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
+               hrmsrgy(indrgy,indrms,ib)=
+     &           hrmsrgy(indrgy,indrms,ib)+expfac
+             endif
+          enddo
 #endif
-c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
-          call restore_parm(iparm)
-          evdw=enetb(21,t,iparm)
-          evdw_t=enetb(1,t,iparm)
-#ifdef SCP14
-          evdw2_14=enetb(17,t,iparm)
-          evdw2=enetb(2,t,iparm)+evdw2_14
+        enddo ! t
+c
+c Thermo and ensemble averages
+c
+        do k=0,nGridT
+          betaT=startGridT+k*delta_T
+          call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
+c            write (iout,*) "ftprim",ftprim
+c            write (iout,*) "ftbis",ftbis
+          betaT=1.0d0/(1.987D-3*betaT)
+c 7/10/18 AL Determine the max Botzmann weights for each temerature
+          call sum_ene(1,iparm,ft,etot)
+          weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
+c          write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
+#ifdef MPI
+          do t=2,scount(me1)
 #else
-          evdw2=enetb(2,t,iparm)
-          evdw2_14=0.0d0
+          do t=2,ntot(islice)
 #endif
-#ifdef SPLITELE
-          ees=enetb(3,t,iparm)
-          evdw1=enetb(16,t,iparm)
-#else
-          ees=enetb(3,t,iparm)
-          evdw1=0.0d0
+            call sum_ene(t,iparm,ft,etot)
+            weight=-betaT*(etot-potEmin)+entfac(t)
+c            write (iout,*) "k",k," t",t," weight",weight
+            if (weight.gt.weimax(k,iparm)) weimax(k,iparm)=weight
+          enddo
+#ifdef MPI
+        enddo
+#ifdef DEBUG
+        write (iout,*) "weimax before REDUCE"
+        write (iout,*) (weimax(k,iparm),k=0,ngridt)
 #endif
-          ecorr=enetb(4,t,iparm)
-          ecorr5=enetb(5,t,iparm)
-          ecorr6=enetb(6,t,iparm)
-          eel_loc=enetb(7,t,iparm)
-          eello_turn3=enetb(8,t,iparm)
-          eello_turn4=enetb(9,t,iparm)
-          eturn6=enetb(10,t,iparm)
-          ebe=enetb(11,t,iparm)
-          escloc=enetb(12,t,iparm)
-          etors=enetb(13,t,iparm)
-          etors_d=enetb(14,t,iparm)
-          ehpb=enetb(15,t,iparm)
-          estr=enetb(18,t,iparm)
-          esccor=enetb(19,t,iparm)
-          edihcnstr=enetb(20,t,iparm)
-          edihcnstr=0.0d0
-          do k=0,nGridT
-            betaT=startGridT+k*delta_T
-            temper=betaT
-c            fT=T0/betaT
-c            ft=2*T0/(T0+betaT)
-            if (rescale_mode.eq.1) then
-              quot=betaT/T0
-              quotl=1.0d0
-              kfacl=1.0d0
-              do l=1,5
-                quotl1=quotl
-                quotl=quotl*quot
-                kfacl=kfacl*kfac
-                denom=kfacl-1.0d0+quotl
-                fT(l)=kfacl/denom
-                ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
-                ftbis(l)=l*kfacl*quotl1*
-     &           (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
-              enddo
-#if defined(FUNCTH)
-              ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
-     &                  320.0d0
-              ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
-             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
-     &              /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
-#elif defined(FUNCT)
-              fT(6)=betaT/T0
-              ftprim(6)=1.0d0/T0
-              ftbis(6)=0.0d0
-#else
-              fT(6)=1.0d0
-              ftprim(6)=0.0d0
-              ftbis(6)=0.0d0
+        do k=0,nGridT
+          weimax_(k)=weimax(k,iparm)
+        enddo
+        call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1,
+     &   MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
+#ifdef DEBUG
+        write (iout,*) "weimax"
+        write (iout,*) (weimax(k,iparm),k=0,ngridt)
 #endif
-            else if (rescale_mode.eq.2) then
-              quot=betaT/T0
-              quotl=1.0d0
-              do l=1,5
-                quotl1=quotl
-                quotl=quotl*quot
-                eplus=dexp(quotl)
-                eminus=dexp(-quotl)
-                logfac=1.0d0/dlog(eplus+eminus)
-                tanhT=(eplus-eminus)/(eplus+eminus)
-                fT(l)=1.12692801104297249644d0*logfac
-                ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
-                ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
-     &          2*l*quotl1/T0*logfac*
-     &          (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
-     &          +ftprim(l)*tanhT)
-              enddo
-#if defined(FUNCTH)
-              ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
-     &                 320.0d0
-              ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
-             ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
-     &               /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
-#elif defined(FUNCT)
-              fT(6)=betaT/T0
-              ftprim(6)=1.0d0/T0
-              ftbis(6)=0.0d0
+        do k=0,nGridT
+          temper=startGridT+k*delta_T
+          betaT=1.0d0/(1.987D-3*temper)
+          call temp_scalfac(temper,ft,ftprim,ftbis,*10)
+          do t=1,scount(me1)
 #else
-              fT(6)=1.0d0
-              ftprim(6)=0.0d0
-              ftbis(6)=0.0d0
+          do t=1,ntot(islice)
 #endif
-            else if (rescale_mode.eq.0) then
-              do l=1,5
-                fT(l)=1.0d0
-                ftprim(l)=0.0d0
-              enddo
-            else
-              write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
-     &          rescale_mode
-              call flush(iout)
-              return1
-            endif
-c            write (iout,*) "ftprim",ftprim
-c            write (iout,*) "ftbis",ftbis
-            betaT=1.0d0/(1.987D-3*betaT)
-#ifdef SPLITELE
-            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
-     &      +wvdwpp*evdw1
-     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
-     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
-     &      +ft(2)*wturn3*eello_turn3
-     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
-     &      +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr
-            eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
-     &            +ftprim(1)*wtor*etors+
-     &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
-     &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
-     &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
-     &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
-     &            ftprim(1)*wsccor*esccor
-            ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
-     &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
-     &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
-     &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
-     &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
-     &            ftbis(1)*wsccor*esccor
+            ind = ind_point(t)
+#ifdef MPI
+            hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
 #else
-            etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
-     &      +ft(1)*welec*(ees+evdw1)
-     &      +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
-     &      +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
-     &      +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
-     &      +ft(2)*wturn3*eello_turn3
-     &      +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
-     &      +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
-     &      +wbond*estr
-            eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
-     &           +ftprim(1)*wtor*etors+
-     &            ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
-     &            ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
-     &            ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
-     &            ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
-     &            ftprim(1)*wsccor*esccor
-            ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
-     &            ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
-     &            ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
-     &            ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
-     &            ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
-     &            ftprim(1)*wsccor*esccor
+            hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
 #endif
-            weight=dexp(-betaT*(etot-potEmin)+entfac(t))
+c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
+c          call restore_parm(iparm)
+            call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
+            weight=dexp(-betaT*(etot-potEmin)+entfac(t)-weimax(k,iparm))
 #ifdef DEBUG
             write (iout,*) "iparm",iparm," t",t," betaT",betaT,
-     &        " etot",etot," entfac",entfac(t),
+     &        " etot",etot," entfac",entfac(t)," boltz",
+     &        -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm),
      &        " weight",weight," ebis",ebis
 #endif
             etot=etot-temper*eprim
@@ -844,38 +836,8 @@ c            write (iout,*) "ftbis",ftbis
      &         +etot*q(j,t)*weight
             enddo
 #endif
-          enddo
-          indE = aint(potE(t,iparm)-aint(potEmin))
-          if (indE.ge.0 .and. indE.le.maxinde) then
-            if (indE.gt.upindE_p) upindE_p=indE
-            histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
-          endif
-#ifdef MPI
-          do ib=1,nT_h(iparm)
-            expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
-            hfin_p(ind,ib)=hfin_p(ind,ib)+
-     &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
-             if (rmsrgymap) then
-               indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy) 
-               indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
-               hrmsrgy_p(indrgy,indrms,ib)=
-     &           hrmsrgy_p(indrgy,indrms,ib)+expfac
-             endif
-          enddo
-#else
-          do ib=1,nT_h(iparm)
-            expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
-            hfin(ind,ib)=hfin(ind,ib)+
-     &       dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
-             if (rmsrgymap) then
-               indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
-               indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
-               hrmsrgy(indrgy,indrms,ib)=
-     &           hrmsrgy(indrgy,indrms,ib)+expfac
-             endif
-          enddo
-#endif
-        enddo ! t
+          enddo ! t
+        enddo ! k
         do ib=1,nT_h(iparm)
           if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
      &    MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
@@ -1094,7 +1056,7 @@ c            write (iout,*) "ftbis",ftbis
           sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
      &     -sumQ(j,i,iparm)*sumE(i,iparm)
         enddo
-        sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
+        sumW(i,iparm)=(-dlog(sumW(i,iparm))-weimax(i,iparm))*(1.987D-3*
      &   (startGridT+i*delta_T))+potEmin
         write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
      &   sumW(i,iparm),sumE(i,iparm)
@@ -1190,5 +1152,359 @@ c            write (iout,*) "ftbis",ftbis
 #endif
 
       return
-
+C#undef DEBUG
+   10 return1
+      end
+c------------------------------------------------------------------------
+      subroutine temp_scalfac(betaT,ft,ftprim,ftbis,*)
+      implicit none
+      include "DIMENSIONS"
+      include "DIMENSIONS.FREE" 
+      include "COMMON.CONTROL"
+      include "COMMON.FREE"
+      include "COMMON.IOUNITS"
+      double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
+     &  kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
+     &  logfac,tanhT,betaT,denom,eplus,eminus
+      integer l
+      if (rescale_mode.eq.1) then
+        quot=betaT/T0
+        quotl=1.0d0
+        kfacl=1.0d0
+        do l=1,5
+          quotl1=quotl
+          quotl=quotl*quot
+          kfacl=kfacl*kfac
+          denom=kfacl-1.0d0+quotl
+          fT(l)=kfacl/denom
+          ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
+          ftbis(l)=l*kfacl*quotl1*
+     &     (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
+        enddo
+#if defined(FUNCTH)
+        ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+     &            320.0d0
+        ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+        ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+     &          /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+        fT(6)=betaT/T0
+        ftprim(6)=1.0d0/T0
+        ftbis(6)=0.0d0
+#else
+        fT(6)=1.0d0
+        ftprim(6)=0.0d0
+        ftbis(6)=0.0d0
+#endif
+      else if (rescale_mode.eq.2) then
+        quot=betaT/T0
+        quotl=1.0d0
+        do l=1,5
+          quotl1=quotl
+          quotl=quotl*quot
+          eplus=dexp(quotl)
+          eminus=dexp(-quotl)
+          logfac=1.0d0/dlog(eplus+eminus)
+          tanhT=(eplus-eminus)/(eplus+eminus)
+          fT(l)=1.12692801104297249644d0*logfac
+          ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
+          ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
+     &    2*l*quotl1/T0*logfac*
+     &    (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
+     &    +ftprim(l)*tanhT)
+        enddo
+#if defined(FUNCTH)
+        ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+     &           320.0d0
+        ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+        ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+     &         /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+        fT(6)=betaT/T0
+        ftprim(6)=1.0d0/T0
+        ftbis(6)=0.0d0
+#else
+        fT(6)=1.0d0
+        ftprim(6)=0.0d0
+        ftbis(6)=0.0d0
+#endif
+      else if (rescale_mode.eq.0) then
+        do l=1,5
+          fT(l)=1.0d0
+          ftprim(l)=0.0d0
+        enddo
+      else
+        write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
+     &    rescale_mode
+        call flush(iout)
+        return1
+      endif
+      return
+      end
+c--------------------------------------------------------------------
+      subroutine sum_ene(t,iparm,ft,etot)
+      implicit none
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
+      include 'COMMON.CONTROL'
+      include 'COMMON.FFIELD'
+      include "COMMON.SBRIDGE"
+      include "COMMON.ENERGIES"
+      include "COMMON.IOUNITS"
+      integer t,iparm
+      double precision fT(6)
+      double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
+     &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
+     &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
+     &  eliptran
+      evdw=enetb(21,t,iparm)
+      evdw_t=enetb(1,t,iparm)
+#ifdef SCP14
+      evdw2_14=enetb(17,t,iparm)
+      evdw2=enetb(2,t,iparm)+evdw2_14
+#else
+      evdw2=enetb(2,t,iparm)
+      evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+      ees=enetb(3,t,iparm)
+      evdw1=enetb(16,t,iparm)
+#else
+      ees=enetb(3,t,iparm)
+      evdw1=0.0d0
+#endif
+      ecorr=enetb(4,t,iparm)
+      ecorr5=enetb(5,t,iparm)
+      ecorr6=enetb(6,t,iparm)
+      eel_loc=enetb(7,t,iparm)
+      eello_turn3=enetb(8,t,iparm)
+      eello_turn4=enetb(9,t,iparm)
+      eturn6=enetb(10,t,iparm)
+      ebe=enetb(11,t,iparm)
+      escloc=enetb(12,t,iparm)
+      etors=enetb(13,t,iparm)
+      etors_d=enetb(14,t,iparm)
+      ehpb=enetb(15,t,iparm)
+      estr=enetb(18,t,iparm)
+      esccor=enetb(19,t,iparm)
+      edihcnstr=enetb(20,t,iparm)
+      edihcnstr=0.0d0
+#ifdef SPLITELE
+      if (shield_mode.gt.0) then
+        etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+     &    +ft(1)*welec*ees
+     &    +ft(1)*wvdwpp*evdw1
+     &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &    +ft(2)*wturn3*eello_turn3
+     &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+     &    +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &    +wbond*estr+wliptran*eliptran
+      else
+        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+     &    +wvdwpp*evdw1
+     &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &    +ft(2)*wturn3*eello_turn3
+     &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+     &    +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &    +wbond*estr+wliptran*eliptran
+      endif
+#else
+      if (shield_mode.gt.0) then
+        etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+     &    +ft(1)*welec*(ees+evdw1)
+     &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &    +ft(2)*wturn3*eello_turn3
+     &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+     &    +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &    +wbond*estr+wliptran*eliptran
+      else
+        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+     &    +ft(1)*welec*(ees+evdw1)
+     &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &    +ft(2)*wturn3*eello_turn3
+     &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+     &    +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &    +wbond*estr+wliptran*eliptran
+      endif
+#endif
+      return
+      end
+c--------------------------------------------------------------------
+      subroutine sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
+      implicit none
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.FREE'
+      include 'COMMON.CONTROL'
+      include 'COMMON.FFIELD'
+      include "COMMON.SBRIDGE"
+      include 'COMMON.ENERGIES'
+      include "COMMON.IOUNITS"
+      integer t,iparm
+      double precision fT(6),fTprim(6),fTbis(6),
+     &  eprim,ebis,temper
+      double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
+     &  escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
+     &  eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
+     &  eliptran
+      evdw=enetb(21,t,iparm)
+      evdw_t=enetb(1,t,iparm)
+#ifdef SCP14
+      evdw2_14=enetb(17,t,iparm)
+      evdw2=enetb(2,t,iparm)+evdw2_14
+#else
+      evdw2=enetb(2,t,iparm)
+      evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+      ees=enetb(3,t,iparm)
+      evdw1=enetb(16,t,iparm)
+#else
+      ees=enetb(3,t,iparm)
+      evdw1=0.0d0
+#endif
+      ecorr=enetb(4,t,iparm)
+      ecorr5=enetb(5,t,iparm)
+      ecorr6=enetb(6,t,iparm)
+      eel_loc=enetb(7,t,iparm)
+      eello_turn3=enetb(8,t,iparm)
+      eello_turn4=enetb(9,t,iparm)
+      eturn6=enetb(10,t,iparm)
+      ebe=enetb(11,t,iparm)
+      escloc=enetb(12,t,iparm)
+      etors=enetb(13,t,iparm)
+      etors_d=enetb(14,t,iparm)
+      ehpb=enetb(15,t,iparm)
+      estr=enetb(18,t,iparm)
+      esccor=enetb(19,t,iparm)
+      edihcnstr=enetb(20,t,iparm)
+      edihcnstr=0.0d0
+#ifdef SPLITELE
+      if (shield_mode.gt.0) then
+        etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+     &    +ft(1)*welec*ees
+     &    +ft(1)*wvdwpp*evdw1
+     &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &    +ft(2)*wturn3*eello_turn3
+     &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+     &    +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &    +wbond*estr+wliptran*eliptran
+        eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
+C     &            +ftprim(6)*evdw_t
+     &    +ftprim(1)*wscp*evdw2
+     &    +ftprim(1)*welec*ees
+     &    +ftprim(1)*wvdwpp*evdw1
+     &    +ftprim(1)*wtor*etors+
+     &    ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+     &    ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+     &    ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+     &    ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+     &    ftprim(1)*wsccor*esccor
+        ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
+     &    +ftbis(1)*wscp*evdw2+
+     &    ftbis(1)*welec*ees
+     &    +ftbis(1)*wvdwpp*evdw
+     &    +ftbis(1)*wtor*etors+
+     &    ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+     &    ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+     &    ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+     &    ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+     &    ftbis(1)*wsccor*esccor
+      else
+        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+     &    +wvdwpp*evdw1
+     &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &    +ft(2)*wturn3*eello_turn3
+     &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+     &    +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &    +wbond*estr+wliptran*eliptran
+        eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
+     &    +ftprim(1)*wtor*etors+
+     &    ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+     &    ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+     &    ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+     &    ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+     &    ftprim(1)*wsccor*esccor
+        ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
+     &    ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+     &    ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+     &    ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+     &    ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+     &    ftbis(1)*wsccor*esccor
+      endif
+#else
+      if (shield_mode.gt.0) then
+        etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+     &    +ft(1)*welec*(ees+evdw1)
+     &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &    +ft(2)*wturn3*eello_turn3
+     &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+     &    +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &    +wbond*estr+wliptran*eliptran
+        eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
+     &    +ftprim(1)*welec*(ees+evdw1)
+     &    +ftprim(1)*wtor*etors+
+     &     ftprim(1)*wscp*evdw2+
+     &     ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+     &     ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+     &     ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+     &     ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+     &     ftprim(1)*wsccor*esccor
+       ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
+     &    +ftbis(1)*wscp*evdw2
+     &    +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
+     &     ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+     &     ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+     &     ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+     &     ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+     &     ftprim(1)*wsccor*esccor
+      else
+        etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+     &    +ft(1)*welec*(ees+evdw1)
+     &    +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c     &    +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+     &    +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+     &    +ft(2)*wturn3*eello_turn3
+     &    +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+     &    +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+     &    +wbond*estr+wliptran*eliptran
+        eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
+     &    +ftprim(1)*wtor*etors+
+     &     ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+     &     ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+     &     ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+     &     ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+     &     ftprim(1)*wsccor*esccor
+        ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
+     &     ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+     &     ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+     &     ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+     &     ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+     &     ftprim(1)*wsccor*esccor
+      endif
+#endif
+      return
       end