modyfication for tube

diff --git a/source/cluster/wham/src-M/probabl.F b/source/cluster/wham/src-M/probabl.F
index 7bc3f1a..67341dd 100644 (file)
--- a/source/cluster/wham/src-M/probabl.F
+++ b/source/cluster/wham/src-M/probabl.F
@@ -296,12 +296,13 @@ c        write (iout,*) "qfree",qfree
         write (iout,*) "ncon", ncon,maxstr_proc
         do i=1,min0(ncon,maxstr_proc)-1 
           sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree 
-#ifdef DEBUG
+C#ifdef DEBUG
+          write (iout,*) "tu szukaj ponizej 7"
           write (iout,*) i,ib,beta_h(ib),
      &     1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),
      &     totfree(list_conf(i)),
      &     -entfac(list_conf(i)),fdimless(i),sumprob
-#endif
+C#endif
           if (sumprob.gt.prob_limit) goto 122
 c          if (sumprob.gt.1.00d0) goto 122
           nlist=nlist+1

diff --git a/source/cluster/wham/src-M/read_coords.F b/source/cluster/wham/src-M/read_coords.F
index 2911e60..55f52e3 100644 (file)
--- a/source/cluster/wham/src-M/read_coords.F
+++ b/source/cluster/wham/src-M/read_coords.F
@@ -393,7 +393,7 @@ c------------------------------------------------------------------------------
       chalen=int((nct-nnt+2)/symetr)
       call int_from_cart1(.false.)
       do j=nnt+1,nct
-        if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)
+        if ((vbld(j).lt.2.0d0 .or. vbld(j).gt.6.0d0)
      &      .and.(itype(j).ne.ntyp1)) then
          if (j.gt.2) then
           if (itel(j).ne.0 .and. itel(j-1).ne.0) then
@@ -419,7 +419,7 @@ c------------------------------------------------------------------------------
       enddo
       do j=nnt,nct
         itj=itype(j)
-        if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0
+        if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(iabs(itj))).gt.5.0d0
      &  .and. itype(j).ne.ntyp1) then
           write (iout,*) "Conformation",jjj,jj+1
           write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j)

diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F
index 2635a2f..3a45231 100644 (file)
--- a/source/unres/src_MD-M/MD_A-MTS.F
+++ b/source/unres/src_MD-M/MD_A-MTS.F
@@ -1653,6 +1653,7 @@ c Removing the velocity of the center of mass
 #endif
       call zerograd
       call etotal(potEcomp)
+      call enerprint(potEcomp)
       if (large) call enerprint(potEcomp)
 #ifdef TIMING_ENE
 #ifdef MPI

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 21c0197..f557ebc 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -11536,7 +11536,8 @@ C the vector between center of side_chain and peptide group
      &pept_group(3),costhet_grad(3),cosphi_grad_long(3),
      &cosphi_grad_loc(3),pep_side_norm(3),side_calf_norm(3)
 C the line belowe needs to be changed for FGPROC>1
-      do i=1,nres-1
+      do i=iatscp_s,iatscp_e
+C      do i=1,nres-1
       if ((itype(i).eq.ntyp1).and.itype(i+1).eq.ntyp1) cycle
       ishield_list(i)=0
 Cif there two consequtive dummy atoms there is no peptide group between them
@@ -11729,8 +11730,8 @@ C lets ommit dummy atoms for now
 C now calculate distance from center of tube and direction vectors
       vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
           if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
-      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxxsize)
-          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize
+      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
       vectube(1)=vectube(1)-tubecenter(1)
       vectube(2)=vectube(2)-tubecenter(2)
 
@@ -11779,8 +11780,8 @@ C      .or.(iti.eq.10)
           vectube(1)=mod(vectube(1),boxxsize)
           if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
           vectube(2)=c(2,i+nres)
-          vectube(2)=mod(vectube(2),boxxsize)
-          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize
+          vectube(2)=mod(vectube(2),boxysize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
 
       vectube(1)=vectube(1)-tubecenter(1)
       vectube(2)=vectube(2)-tubecenter(2)
@@ -11864,8 +11865,8 @@ C lets ommit dummy atoms for now
 C now calculate distance from center of tube and direction vectors
       vectube(1)=mod((c(1,i)+c(1,i+1))/2.0d0,boxxsize)
           if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
-      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxxsize)
-          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize
+      vectube(2)=mod((c(2,i)+c(2,i+1))/2.0d0,boxysize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
       vectube(1)=vectube(1)-tubecenter(1)
       vectube(2)=vectube(2)-tubecenter(2)
 
@@ -11914,8 +11915,8 @@ C in UNRES uncomment the line below as GLY has no side-chain...
           vectube(1)=mod(vectube(1),boxxsize)
           if (vectube(1).lt.0) vectube(1)=vectube(1)+boxxsize
           vectube(2)=c(2,i+nres)
-          vectube(2)=mod(vectube(2),boxxsize)
-          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxxsize
+          vectube(2)=mod(vectube(2),boxysize)
+          if (vectube(2).lt.0) vectube(2)=vectube(2)+boxysize
 
       vectube(1)=vectube(1)-tubecenter(1)
       vectube(2)=vectube(2)-tubecenter(2)
@@ -11937,10 +11938,10 @@ C lipbufthick is thickenes of lipid buffore
          ssgradtube=-sscagradlip(fracinbuf)/tubebufthick
          print *,ssgradtube, sstube,tubetranene(itype(i))
          enetube(i+nres)=enetube(i+nres)+sstube*tubetranene(itype(i))
-         gg_tube_SC(3,i)=gg_tube_SC(3,i)
-     &+ssgradtube*tubetranene(itype(i))
-         gg_tube(3,i-1)= gg_tube(3,i-1)
-     &+ssgradtube*tubetranene(itype(i))
+C         gg_tube_SC(3,i)=gg_tube_SC(3,i)
+C     &+ssgradtube*tubetranene(itype(i))
+C         gg_tube(3,i-1)= gg_tube(3,i-1)
+C     &+ssgradtube*tubetranene(itype(i))
 C         print *,"doing sccale for lower part"
         elseif (positi.gt.buftubetop) then
          fracinbuf=1.0d0-

diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F
index 5001b6b..8fa2bad 100644 (file)
--- a/source/unres/src_MD-M/parmread.F
+++ b/source/unres/src_MD-M/parmread.F
@@ -1361,7 +1361,7 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
         enddo
       enddo
       write(iout,*) "tube param"
-      read(itube,*) epspeptube,sigmapeptube
+      read(itube,*) epspeptube,sigmapeptube,tubetranenepep
       sigmapeptube=sigmapeptube**6
       sigeps=dsign(1.0D0,epspeptube)
       epspeptube=dabs(epspeptube)
@@ -1369,7 +1369,7 @@ c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
       pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
       write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
       do i=1,ntyp
-       read(itube,*) epssctube,sigmasctube
+       read(itube,*) epssctube,sigmasctube,tubetranene(i)
        sigmasctube=sigmasctube**6
        sigeps=dsign(1.0D0,epssctube)
        epssctube=dabs(epssctube)

diff --git a/source/wham/src-M/enecalc1.F b/source/wham/src-M/enecalc1.F
index a5d25b3..3a5dbf3 100644 (file)
--- a/source/wham/src-M/enecalc1.F
+++ b/source/wham/src-M/enecalc1.F
@@ -744,7 +744,7 @@ c------------------------------------------------------------------------------
       call int_from_cart1(.false.)
       do j=nnt+1,nct
         if (itype(j-1).ne.ntyp1 .and. itype(j).ne.ntyp1 .and. 
-     &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0)) then
+     &    (vbld(j).lt.2.0d0 .or. vbld(j).gt.6.0d0)) then
           if (iprint.gt.0) 
      &    write (iout,*) "Bad CA-CA bond length",j," ",vbld(j),
      &      " for conformation",ii
@@ -769,7 +769,7 @@ c------------------------------------------------------------------------------
       do j=nnt,nct
         itj=itype(j)
         if (itype(j).ne.10 .and.itype(j).ne.ntyp1 .and. 
-     &     (vbld(nres+j)-dsc(iabs(itj))).gt.2.0d0) then
+     &     (vbld(nres+j)-dsc(iabs(itj))).gt.5.0d0) then
           if (iprint.gt.0) 
      &    write (iout,*) "Bad CA-SC bond length",j," ",vbld(nres+j),
      &     " for conformation",ii

diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F
index f77a245..f0180ba 100644 (file)
--- a/source/wham/src-M/wham_calc1.F
+++ b/source/wham/src-M/wham_calc1.F
@@ -51,7 +51,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),
+     &  fi_p_min(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,7 +64,8 @@ c      parameter (MaxHdim=200000)
      & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,
      & 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 potEmin_t,entmin_p,entmax_p,
+     & potEmin_t_all(maxT_h,Max_Parm)
       integer histent_p(0:2000)
       logical lprint /.true./
 #endif
@@ -75,11 +77,12 @@ 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),
+     & fi_min(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
+     & potEmin,ent,potEmin_all(maxT_h,Max_Parm),potEmin_min,
+     & hfin_ent(0:MaxHdim),vmax,aux,entfac_min
       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
@@ -107,6 +110,11 @@ c      parameter (MaxHdim=200000)
       dmin=0.0d0
       tmax=0
       potEmin=1.0d10
+       do i=1,nParmset
+         do j=1,nT_h(i)
+           potEmin_all(j,i)=1.0d10
+         enddo
+       enddo
       rgymin=1.0d10
       rmsmin=1.0d10
       rgymax=0.0d0
@@ -120,9 +128,9 @@ C#define DEBUG
 #else
       do i=1,ntot(islice)
 #endif
-        do j=1,nParmSet
-          if (potE(i,j).le.potEmin) potEmin=potE(i,j)
-        enddo
+C        do j=1,nParmSet
+C          if (potE(i,j).le.potEmin) potEmin=potE(i,j)
+C        enddo
         if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i)
         if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i)
         if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i)
@@ -187,8 +195,8 @@ c     &      ind_point(i)
       call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,
      &  WHAM_COMM,IERROR)
       tmax=tmax_t
-      call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
-     &  MPI_MIN,WHAM_COMM,IERROR)
+C      call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,
+C     &  MPI_MIN,WHAM_COMM,IERROR)
       call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,
      &  MPI_MIN,WHAM_COMM,IERROR)
       call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,
@@ -197,7 +205,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
+
       rgymin=rgymin_t
       rgymax=rgymax_t
       rmsmin=rmsmin_t
@@ -392,7 +401,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &           *(dd-q0(j,kk,ib,iparm))**2
               enddo
               v(i,kk,ib,iparm)=
-     &          -beta_h(ib,iparm)*(etot-potEmin+Econstr)
+     &          -beta_h(ib,iparm)*(etot+Econstr)
 #ifdef DEBUG
               write (iout,'(4i5,4e15.5)') i,kk,ib,iparm,
      &         etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
@@ -408,6 +417,8 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
 ! Compute new free-energy values corresponding to the righ-hand side of the 
 ! equation and their derivatives.
         write (iout,*) "------------------------fi"
+        entfac_min=1.0d10
+
 #ifdef MPI
         do t=1,scount(me1)
 #else
@@ -433,10 +444,52 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             enddo
           enddo
           entfac(t)=-dlog(denom)-vmax
+          if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
 #ifdef DEBUG
           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
+              fi_p_min(ii,iib,iparm)=-1.0d10
+              do t=1,scount(me)
+                aux=v(t,ii,iib,iparm)+entfac(t)
+                if (aux.gt.fi_p_min(ii,iib,iparm))
+     &                                   fi_p_min(ii,iib,iparm)=aux
+              enddo
+#else
+              do t=1,ntot(islice)
+                aux=v(t,ii,iib,iparm)+entfac(t)
+                if (aux.gt.fi_min(ii,iib,iparm)) 
+     &                                     fi_min(ii,iib,iparm)=aux
+              enddo
+#endif
+            enddo ! ii
+          enddo ! iib
+        enddo ! iparm
+#ifdef MPI
+#ifdef DEBUG
+        write (iout,*) "fi_min before AllReduce"
+        do i=1,nParmSet
+          do j=1,nT_h(i)
+            write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i))
+          enddo
+        enddo
+#endif
+        call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet,
+     &         MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
+#ifdef DEBUG
+        write (iout,*) "fi_min after AllReduce"
+        do i=1,nParmSet
+          do j=1,nT_h(i)
+            write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
+          enddo
+        enddo
+#endif
+#endif
         do iparm=1,nParmSet
           do iib=1,nT_h(iparm)
             do ii=1,nR(iib,iparm)
@@ -444,23 +497,23 @@ 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)-fi_min(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)
+              write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,
+     &         v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm),
+     &         fi_p(ii,iib,iparm)
 #endif
               enddo
 #else
               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)-fi_min(ii,iib,iparm))
               enddo
 #endif
             enddo ! ii
           enddo ! iib
         enddo ! iparm
-
 #ifdef MPI
 #ifdef DEBUG
         write (iout,*) "fi before MPI_Reduce me",me,' master',master
@@ -495,7 +548,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))-fi_min(i,ib,iparm)
               avefi=avefi+fi(i,ib,iparm)
             enddo
           enddo
@@ -543,6 +596,198 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
       enddo ! iter
 
    20 continue
+
+#ifdef MPI
+      do i=1,scount(me1)
+#else
+      do i=1,ntot(islice)
+#endif
+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)
+#endif
+          call restore_parm(iparm)
+#ifdef DEBUG
+          write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,
+     &      wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,
+     &      wtor_d,wsccor,wbond
+#endif
+          do ib=1,nT_h(iparm)
+            if (rescale_mode.eq.1) then
+              quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+              quotl=1.0d0
+              kfacl=1.0d0
+              do l=1,5
+                quotl1=quotl
+                quotl=quotl*quot
+                kfacl=kfacl*kfac
+                fT(l)=kfacl/(kfacl-1.0d0+quotl)
+              enddo
+#if defined(FUNCTH)
+              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+              ft(6)=1.0d0
+#endif
+            else if (rescale_mode.eq.2) then
+              quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+              quotl=1.0d0
+              do l=1,5
+                quotl=quotl*quot
+                fT(l)=1.12692801104297249644d0/
+     &             dlog(dexp(quotl)+dexp(-quotl))
+              enddo
+#if defined(FUNCTH)
+              tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+              ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+              ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+              ft(6)=1.0d0
+#endif
+c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
+            else if (rescale_mode.eq.0) then
+              do l=1,6
+                fT(l)=1.0d0
+              enddo
+            else
+              write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
+     &          rescale_mode
+              call flush(iout)
+              return1
+            endif
+          evdw=enetb(21,i,iparm)
+          evdw_t=enetb(1,i,iparm)
+#ifdef SCP14
+          evdw2_14=enetb(17,i,iparm)
+          evdw2=enetb(2,i,iparm)+evdw2_14
+#else
+          evdw2=enetb(2,i,iparm)
+          evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+          ees=enetb(3,i,iparm)
+          evdw1=enetb(16,i,iparm)
+#else
+          ees=enetb(3,i,iparm)
+          evdw1=0.0d0
+#endif
+          ecorr=enetb(4,i,iparm)
+          ecorr5=enetb(5,i,iparm)
+          ecorr6=enetb(6,i,iparm)
+          eel_loc=enetb(7,i,iparm)
+          eello_turn3=enetb(8,i,iparm)
+          eello_turn4=enetb(9,i,iparm)
+          eturn6=enetb(10,i,iparm)
+          ebe=enetb(11,i,iparm)
+          escloc=enetb(12,i,iparm)
+          etors=enetb(13,i,iparm)
+          etors_d=enetb(14,i,iparm)
+          ehpb=enetb(15,i,iparm)
+          estr=enetb(18,i,iparm)
+          esccor=enetb(19,i,iparm)
+          edihcnstr=enetb(20,i,iparm)
+C          edihcnstr=0.0d0
+          eliptran=enetb(22,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
+     &      +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
+     &      +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
+     &      +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+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
+     &      +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
+            etot=etot-entfac(i)/beta_h(ib,iparm)
+            if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
+
+          enddo   ! ib
+        enddo     ! iparm
+      enddo       ! i
+
+
+#ifdef DEBUG
+      write (iout,*) "The potEmin array before reduction"
+      do i=1,nParmSet
+        write (iout,*) "Parameter set",i
+        do j=1,nT_h(i)
+          write (iout,*) j,PotEmin_all(j,i)
+        enddo
+      enddo
+      write (iout,*) "potEmin_min",potEmin_min
+#endif
+#ifdef MPI
+C Determine the minimum energes for all parameter sets and temperatures
+      call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1),
+     &  maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
+      do i=1,nParmSet
+        do j=1,nT_h(i)
+          potEmin_all(j,i)=potEmin_t_all(j,i)
+        enddo
+      enddo
+#endif
+      potEmin_min=potEmin_all(1,1)
+      do i=1,nParmSet
+        do j=1,nT_h(i)
+          if (potEmin_all(j,i).lt.potEmin_min)
+     &             potEmin_min=potEmin_all(j,i)
+        enddo
+      enddo
+#ifdef DEBUG
+      write (iout,*) "The potEmin array"
+      do i=1,nParmSet
+        write (iout,*) "Parameter set",i
+        do j=1,nT_h(i)
+          write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)),
+     &      PotEmin_all(j,i)
+        enddo
+      enddo
+      write (iout,*) "potEmin_min",potEmin_min
+#endif
+
+
 ! Now, put together the histograms from all simulations, in order to get the
 ! unbiased total histogram.
 #ifdef MPI
@@ -717,7 +962,9 @@ c          write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
           estr=enetb(18,t,iparm)
           esccor=enetb(19,t,iparm)
           edihcnstr=enetb(20,t,iparm)
-          edihcnstr=0.0d0
+C          edihcnstr=0.0d0
+          eliptran=enetb(22,i,iparm)
+
           do k=0,nGridT
             betaT=startGridT+k*delta_T
             temper=betaT
@@ -798,6 +1045,36 @@ c            ft=2*T0/(T0+betaT)
 c            write (iout,*) "ftprim",ftprim
 c            write (iout,*) "ftbis",ftbis
             betaT=1.0d0/(1.987D-3*betaT)
+            if (betaT.ge.beta_h(1,iparm)) then
+              potEmin=potEmin_all(1,iparm)+
+     &         (potEmin_all(1,iparm)-potEmin_all(2,iparm))/
+     &         (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))*
+     &         (1.0/betaT-1.0/beta_h(1,iparm))
+#ifdef DEBUG
+              write(iout,*) "first",temper,potEmin
+#endif
+            else if (betaT.le.beta_h(nT_h(iparm),iparm)) then
+              potEmin=potEmin_all(nT_h(iparm),iparm)+
+     &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/
+     &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))*
+     &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm))
+#ifdef DEBUG
+              write (iout,*) "last",temper,potEmin
+#endif
+            else
+              do l=1,nT_h(iparm)-1
+                if (betaT.le.beta_h(l,iparm) .and.
+     &              betaT.gt.beta_h(l+1,iparm)) then
+                  potEmin=potEmin_all(l,iparm)
+#ifdef DEBUG
+                  write (iout,*) "l",l,
+     &             betaT,1.0d0/(1.987D-3*beta_h(l,iparm)),
+     &             1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin
+#endif
+                  exit
+                endif
+              enddo
+            endif
 #ifdef SPLITELE
       if (shield_mode.gt.0) then
             etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
@@ -882,7 +1159,7 @@ C     &            +ftprim(6)*evdw_t
      &            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
+     &            ftbis(1)*wsccor*esccor
        else
             etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
      &      +ft(1)*welec*(ees+evdw1)
@@ -905,7 +1182,7 @@ C     &            +ftprim(6)*evdw_t
      &            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
+     &            ftbis(1)*wsccor*esccor
 
        endif
 
@@ -948,6 +1225,7 @@ C     &            +ftprim(6)*evdw_t
           endif
 #ifdef MPI
           do ib=1,nT_h(iparm)
+            potEmin=potEmin_all(ib,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))
@@ -960,6 +1238,7 @@ C     &            +ftprim(6)*evdw_t
           enddo
 #else
           do ib=1,nT_h(iparm)
+            potEmin=potEmin_all(ib,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))
@@ -1178,6 +1457,21 @@ C     &            +ftprim(6)*evdw_t
         write (iout,'(a,i3)') "Parameter set",iparm
       endif
       do i=0,NGridT
+        betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
+        if (betaT.ge.beta_h(1,iparm)) then
+          potEmin=potEmin_all(1,iparm)
+        else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
+          potEmin=potEmin_all(nT_h(iparm),iparm)
+        else
+          do l=1,nT_h(iparm)-1
+            if (betaT.le.beta_h(l,iparm) .and.
+     &          betaT.gt.beta_h(l+1,iparm)) then
+              potEmin=potEmin_all(l,iparm)
+              exit
+            endif
+          enddo
+        endif
+
         sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
         sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/
      &    sumW(i,iparm)