Water micro and bere and lang with gly working with D lang not
[unres4.git] / source / wham / wham_calc.F90
index 497ba0c..8a25261 100644 (file)
@@ -79,7 +79,7 @@
               sumQsq_p,sumEQ_p,sumEprim_p !(MaxQ1,0:nGridT,Max_Parm) 
       real(kind=8) :: hfin_p(0:MaxHdim,maxT_h),&
               hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,&
-              hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
+              hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h),weimax_(0:ngridT)
       real(kind=8) :: rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
       real(kind=8) :: potEmin_t,potEmin_t_all(maxT_h,Max_Parm)!,entmin_p,entmax_p
 !      integer :: histent_p(0:2000)
@@ -98,7 +98,8 @@
               hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),&
               potEmin,ent,&
               hfin_ent(0:MaxHdim),vmax,aux,fi_min(MaxR,maxT_h,nParmSet), &
-              potEmin_all(maxT_h,Max_Parm),potEmin_min,entfac_min
+              potEmin_all(maxT_h,Max_Parm),potEmin_min,entfac_min,&
+              weimax(0:nGridT,Max_Parm)
       real(kind=8) :: 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
         ecationcation,ecation_prot, evdwpp,eespp ,evdwpsb,eelpsb, &
         evdwsb, eelsb, estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
         ecorr_nucl, ecorr3_nucl,epeppho, escpho, epepbase,escbase,ecation_nucl,&
-        elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang,eliptran
-
+        elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang,eliptran,&
+        ehomology_constr
 
 
       integer :: ind_point(maxpoint),upindE,indE
             ecat_prottran=enetb(56,i,iparm)
             ecation_protang=enetb(57,i,iparm)
             eliptran=enetb(22,i,iparm)
+            ehomology_constr=enetb(51,i,iparm)
 #ifdef DEBUG
             write (iout,'(3i5,6f5.2,15f12.3)') i,ib,iparm,(ft(l),l=1,6),&
              evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
-             etors,etors_d,eello_turn3,eello_turn4,esccor,ecationcation
+             etors,etors_d,eello_turn3,eello_turn4,esccor,ecationcation,&
+             ehomology_constr
 #endif
 
 !#ifdef SPLITELE
 !            +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 &
+!            +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
 !            +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 !            +wbond*estr
 !#else
 !            +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(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
 !            +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
 !            +wbond*estr
 !#endif
             call enerprint(energia(0))
             endif
 #endif
+            write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
+            write (iout,*) "waga_homology", waga_homology
+#ifdef DEBUG
+            write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
+#endif
+            if (homol_nset.gt.1) then
+
+              do kk=1,nR(ib,iparm)
+                Econstr=waga_homology(kk)*ehomology_constr
+                v(i,kk,ib,iparm)= &
+                 -beta_h(ib,iparm)*(etot+Econstr)
+#ifdef DEBUG
+                write (iout,'(4i5,4e15.5)') i,kk,ib,iparm, &
+                etot,Econstr,v(i,kk,ib,iparm)
+#endif
+              enddo ! kk
+
+            else
+
+              etot=etot+ehomology_constr
+
+
             do kk=1,nR(ib,iparm)
               Econstr=0.0d0
               do j=1,nQ
                etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
 #endif
             enddo ! kk
+           endif
           enddo   ! ib
         enddo     ! iparm
       enddo       ! i
             ecat_prottran=enetb(56,i,iparm)
             ecation_protang=enetb(57,i,iparm)
             eliptran=enetb(22,i,iparm)
+            ehomology_constr=enetb(51,i,iparm)
 #ifdef SPLITELE
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees &
             +wvdwpp*evdw1 &
             ecat_prottran=enetb(56,t,iparm)
             ecation_protang=enetb(57,t,iparm)
             eliptran=enetb(22,t,iparm)
+            ehomology_constr=enetb(51,t,iparm)
+
           do k=0,nGridT
             betaT=startGridT+k*delta_T
             temper=betaT
                  +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
 
 #endif
-            weight=dexp(-betaT*(etot-potEmin)+entfac(t))
-#ifdef DEBUG
-            write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
-              " etot",etot," entfac",entfac(t),&
-              " weight",weight," ebis",ebis
-#endif
-            etot=etot-temper*eprim
-#ifdef MPI
-            sumW_p(k,iparm)=sumW_p(k,iparm)+weight
-            sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
-            sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
-            sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
-            do j=1,nQ+2
-              sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
-              sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
-              sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
-               +etot*q(j,t)*weight
-            enddo
-#else
-            sumW(k,iparm)=sumW(k,iparm)+weight
-            sumE(k,iparm)=sumE(k,iparm)+etot*weight
-            sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
-            sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
-            do j=1,nQ+2
-              sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
-              sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
-              sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
-               +etot*q(j,t)*weight
-            enddo
-#endif
+!            weight=dexp(-betaT*(etot-potEmin)+entfac(t))
+!#ifdef DEBUG
+!            write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
+!              " etot",etot," entfac",entfac(t),&
+!              " weight",weight," ebis",ebis
+!#endif
+!            etot=etot-temper*eprim
+!#ifdef MPI
+!            sumW_p(k,iparm)=sumW_p(k,iparm)+weight
+!            sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
+!            sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
+!            sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
+!            do j=1,nQ+2
+!              sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
+!              sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
+!              sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
+!               +etot*q(j,t)*weight
+!            enddo
+!#else
+!            sumW(k,iparm)=sumW(k,iparm)+weight
+!            sumE(k,iparm)=sumE(k,iparm)+etot*weight
+!            sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
+!            sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
+!            do j=1,nQ+2
+!              sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
+!              sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
+!              sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
+!               +etot*q(j,t)*weight
+!            enddo
+!#endif
           enddo
           indE = aint(potE(t,iparm)-aint(potEmin))
           if (indE.ge.0 .and. indE.le.maxinde) then
           enddo
 #endif
         enddo ! t
+!
+! Thermo and ensemble averages
+!
+        do k=0,nGridT
+          betaT=startGridT+k*delta_T
+!          call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
+      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
+
+!            write (iout,*) "ftprim",ftprim
+!            write (iout,*) "ftbis",ftbis
+          betaT=1.0d0/(1.987D-3*betaT)
+! 7/10/18 AL Determine the max Botzmann weights for each temerature
+!          call sum_ene(1,iparm,ft,etot)
+      evdw=enetb(1,1,iparm)
+      evdw_t=enetb(20,1,iparm)
+#ifdef SCP14
+      evdw2_14=enetb(18,1,iparm)
+      evdw2=enetb(2,1,iparm)+evdw2_14
+#else
+      evdw2=enetb(2,1,iparm)
+      evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+      ees=enetb(3,1,iparm)
+      evdw1=enetb(16,1,iparm)
+#else
+      ees=enetb(3,1,iparm)
+      evdw1=0.0d0
+#endif
+      ecorr=enetb(4,1,iparm)
+      ecorr5=enetb(5,1,iparm)
+      ecorr6=enetb(6,1,iparm)
+      eel_loc=enetb(7,1,iparm)
+      eello_turn3=enetb(8,1,iparm)
+      eello_turn4=enetb(9,1,iparm)
+      eello_turn6=enetb(10,1,iparm)
+      ebe=enetb(11,1,iparm)
+      escloc=enetb(12,1,iparm)
+      etors=enetb(13,1,iparm)
+      etors_d=enetb(14,1,iparm)
+      ehpb=enetb(15,1,iparm)
+      estr=enetb(17,1,iparm)
+      esccor=enetb(21,1,iparm)
+      edihcnstr=enetb(19,1,iparm)
+!      eliptran=enetb(22,1,iparm)
+!      esaxs=enetb(26,1,iparm)
+      ecationcation=enetb(41,1,iparm)
+      ecation_prot=enetb(42,1,iparm)
+      evdwpp =      enetb(26,1,iparm)
+      eespp =      enetb(27,1,iparm)
+      evdwpsb =      enetb(28,1,iparm)
+      eelpsb =      enetb(29,1,iparm)
+      evdwsb =      enetb(30,1,iparm)
+      eelsb =      enetb(31,1,iparm)
+      estr_nucl =      enetb(32,1,iparm)
+      ebe_nucl =      enetb(33,1,iparm)
+      esbloc =      enetb(34,1,iparm)
+      etors_nucl =      enetb(35,1,iparm)
+      etors_d_nucl =      enetb(36,1,iparm)
+      ecorr_nucl =      enetb(37,1,iparm)
+      ecorr3_nucl =      enetb(38,1,iparm)
+      epeppho=   enetb(49,1,iparm)
+      escpho=    enetb(48,1,iparm)
+      epepbase=  enetb(47,1,iparm)
+      escbase=   enetb(46,1,iparm)
+      ecation_nucl= enetb(50,1,iparm)
+      ehomology_constr=enetb(51,1,iparm)
+
+#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*eello_turn6+ft(2)*wel_loc*eel_loc &
+          +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+            + wcatnucl*ecation_nucl
+
+          
+      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
+          +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*eello_turn6+ft(2)*wel_loc*eel_loc &
+          +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+            + wcatnucl*ecation_nucl
+
+      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 &
+!     &    +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+          +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+            + wcatnucl*ecation_nucl
+
+      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
+          +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+          +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+            + wcatnucl*ecation_nucl
+
+      endif
+#endif
+
+          weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
+!          write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
+#ifdef MPI
+          do t=2,scount(me1)
+#else
+          do t=2,ntot(islice)
+#endif
+!           call sum_ene(t,iparm,ft,etot)
+      evdw=enetb(1,t,iparm)
+      evdw_t=enetb(20,t,iparm)
+#ifdef SCP14
+      evdw2_14=enetb(18,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)
+      eello_turn6=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(17,t,iparm)
+      esccor=enetb(21,t,iparm)
+      edihcnstr=enetb(19,t,iparm)
+!      eliptran=enetb(22,t,iparm)
+!      esaxs=enetb(26,t,iparm)
+      ecationcation=enetb(41,t,iparm)
+      ecation_prot=enetb(42,t,iparm)
+      evdwpp =      enetb(26,t,iparm)
+      eespp =      enetb(27,t,iparm)
+      evdwpsb =      enetb(28,t,iparm)
+      eelpsb =      enetb(29,t,iparm)
+      evdwsb =      enetb(30,t,iparm)
+      eelsb =      enetb(31,t,iparm)
+      estr_nucl =      enetb(32,t,iparm)
+      ebe_nucl =      enetb(33,t,iparm)
+      esbloc =      enetb(34,t,iparm)
+      etors_nucl =      enetb(35,t,iparm)
+      etors_d_nucl =      enetb(36,t,iparm)
+      ecorr_nucl =      enetb(37,t,iparm)
+      ecorr3_nucl =      enetb(38,t,iparm)
+      epeppho=   enetb(49,t,iparm)
+      escpho=    enetb(48,t,iparm)
+      epepbase=  enetb(47,t,iparm)
+      escbase=   enetb(46,t,iparm)
+      ecation_nucl=enetb(50,t,iparm)
+      ehomology_constr=enetb(51,t,iparm)
+
+#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*eello_turn6+ft(2)*wel_loc*eel_loc &
+          +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+            + wcatnucl*ecation_nucl
+
+      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
+          +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*eello_turn6+ft(2)*wel_loc*eel_loc &
+          +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+            + wcatnucl*ecation_nucl
+      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 &
+!     &    +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+          +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+            + wcatnucl*ecation_nucl
+      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
+          +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+          +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+            +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+            + wcatnucl*ecation_nucl
+      endif
+#endif
+
+            weight=-betaT*(etot-potEmin)+entfac(t)
+!            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
+        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
+        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
+        do k=0,nGridT
+          temper=startGridT+k*delta_T
+          betaT=1.0d0/(1.987D-3*temper)
+!          call temp_scalfac(temper,ft,ftprim,ftbis,*10)
+      if (rescale_mode.eq.1) then
+        quot=temper/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((temper-320.0d0)/80.0d0))/ &
+                 320.0d0
+        ftprim(6)=1.0d0/(320.0d0*dcosh((temper-320.0d0)/80.0d0)**2)
+        ftbis(6)=-2.0d0*dtanh((temper-320.0d0)/80.0d0) &
+               /(320.0d0*80.0d0*dcosh((temper-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+        fT(6)=temper/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=temper/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((temper-320.0d0)/80.0d0))/ &
+                320.0d0
+        ftprim(6)=1.0d0/(320.0d0*dcosh((temper-320.0d0)/80.0d0)**2)
+        ftbis(6)=-2.0d0*dtanh((temper-320.0d0)/80.0d0) &
+              /(320.0d0*80.0d0*dcosh((temper-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+        fT(6)=temper/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
+
+
+
+          do t=1,scount(me1)
+#else
+          do t=1,ntot(islice)
+#endif
+            ind = ind_point(t)
+#ifdef MPI
+            hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
+#else
+            hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
+#endif
+!          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
+!          call restore_parm(iparm)
+!            call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
+      evdw=enetb(1,t,iparm)
+      evdw_t=enetb(20,t,iparm)
+#ifdef SCP14
+      evdw2_14=enetb(18,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)
+      eello_turn6=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(17,t,iparm)
+      esccor=enetb(21,t,iparm)
+      edihcnstr=enetb(19,t,iparm)
+!      eliptran=enetb(22,t,iparm)
+!      esaxs=enetb(26,t,iparm)
+!      ehomology_constr=enetb(27,t,iparm)
+!      edfadis=enetb(28,t,iparm)
+!      edfator=enetb(29,t,iparm)
+!      edfanei=enetb(30,t,iparm)
+!      edfabet=enetb(31,t,iparm)
+          ecationcation=enetb(41,i,iparm)
+          ecation_prot=enetb(42,i,iparm)
+            evdwpp =      enetb(26,i,iparm)
+            eespp =      enetb(27,i,iparm)
+            evdwpsb =      enetb(28,i,iparm)
+            eelpsb =      enetb(29,i,iparm)
+            evdwsb =      enetb(30,i,iparm)
+            eelsb =      enetb(31,i,iparm)
+            estr_nucl =      enetb(32,i,iparm)
+            ebe_nucl =      enetb(33,i,iparm)
+            esbloc =      enetb(34,i,iparm)
+            etors_nucl =      enetb(35,i,iparm)
+            etors_d_nucl =      enetb(36,i,iparm)
+            ecorr_nucl =      enetb(37,i,iparm)
+            ecorr3_nucl =      enetb(38,i,iparm)
+            epeppho=   enetb(49,i,iparm)
+            escpho=    enetb(48,i,iparm)
+            epepbase=  enetb(47,i,iparm)
+            escbase=   enetb(46,i,iparm)
+            ecation_nucl= enetb(50,i,iparm)
+            ehomology_constr=enetb(51,i,iparm)
+
+#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*eello_turn6+ft(2)*wel_loc*eel_loc &
+         +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+         +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+         +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+         +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+         +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+         *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+         +wcorr3_nucl*ecorr3_nucl&
+         +wscbase*escbase&
+         +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+         + wcatnucl*ecation_nucl
+
+!        eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
+!!     &            +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*eello_turn6+
+!     &    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*eello_turn6+
+!     &    ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+!     &    ftbis(1)*wsccor*esccor
+
+            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*eello_turn6+ &
+                 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
+                 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
+            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*eello_turn6+ &
+                 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
+                 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
+
+
+      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*eello_turn6+ft(2)*wel_loc*eel_loc &
+         +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+         +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+            + wcatnucl*ecation_nucl
+
+
+            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*eello_turn6+ &
+                 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
+                 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
+            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*eello_turn6+ &
+                 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
+                 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
+
+
+!       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*eello_turn6+
+!    &    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*eello_turn6+
+!    &    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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+         +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+         +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+            + wcatnucl*ecation_nucl
+
+
+            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*eello_turn6+ &
+                 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
+                 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
+            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*eello_turn6+ &
+                 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
+                 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
+
+
+
+!       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*eello_turn6+
+!    &     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*eello_turn6+
+!    &     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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+         +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+         +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation& 
+            +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+            +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+            +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+            *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+            +wcorr3_nucl*ecorr3_nucl&
+            +wscbase*escbase&
+            +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+            + wcatnucl*ecation_nucl
+
+
+            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*eello_turn6+ &
+                 ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
+                 ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
+            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*eello_turn6+ &
+                 ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
+                 ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+                 +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
+
+
+
+
+!       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*eello_turn6+
+!    &     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*eello_turn6+
+!    &     ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+!    &     ftprim(1)*wsccor*esccor
+      endif
+#endif
+
+            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)," boltz", &
+             -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm), &
+             " weight",weight," ebis",ebis
+#endif
+            etot=etot-temper*eprim
+#ifdef MPI
+            sumW_p(k,iparm)=sumW_p(k,iparm)+weight
+            sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
+            sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
+            sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
+            do j=1,nQ+2
+              sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
+              sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
+              sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
+              +etot*q(j,t)*weight
+            enddo
+#else
+            sumW(k,iparm)=sumW(k,iparm)+weight
+            sumE(k,iparm)=sumE(k,iparm)+etot*weight
+            sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
+            sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
+            do j=1,nQ+2
+              sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
+              sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
+              sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
+              +etot*q(j,t)*weight
+            enddo
+#endif
+          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)