Fixed a MAJOR error in WHAM; energies were computed inconsistently

diff --git a/source/wham/src/enecalc1.F b/source/wham/src/enecalc1.F
index 3d2b6c6..5ef4d77 100644 (file)
--- a/source/wham/src/enecalc1.F
+++ b/source/wham/src/enecalc1.F
@@ -204,10 +204,16 @@ c        call pdbout(ii+1,beta_h(ib,ipar),efree,energia(0),0.0d0,rmsdev)
      &         iii+1,indstart(me1)+iii," T",
      &         1.0d0/(1.987D-3*beta_h(ib,ipar))
              call enerprint(energia(0),fT)
-             itmp=ipdb
-             ipdb=iout
-             call pdbout(iii+1,beta_h(ib,ipar),
-     &                   eini,energia(0),0.0d0,rmsdev)
+            write (iout,'(4f10.5,2i5)') 0.0,energia(0),0.0,
+     &       1.0d0/(beta_h(ib,ipar)*1.987D-3),
+     &       0,0
+            write(iout,'(8f10.5)')
+     &       ((c(l,k),l=1,3),k=1,nres),
+     &       ((c(l,k+nres),l=1,3),k=nnt,nct)
+c             itmp=ipdb
+c             ipdb=iout
+c             call pdbout(iii+1,beta_h(ib,ipar),
+c     &                   eini,energia(0),0.0d0,rmsdev)
              write (iout,*)
              ipdb=itmp
 

diff --git a/source/wham/src/wham_calc1.F b/source/wham/src/wham_calc1.F
index ec46232..d9c3de2 100644 (file)
--- a/source/wham/src/wham_calc1.F
+++ b/source/wham/src/wham_calc1.F
@@ -313,7 +313,7 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
             edfanei=enetb(25,i,iparm)
             edfabet=enetb(26,i,iparm)
 #ifdef DEBUG
-            write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
+            write (iout,'(3i5,6f5.2,16f12.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,
      &       ehomology_constr,edfadis,edfator,edfanei,edfabet
@@ -674,12 +674,13 @@ 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)
+            ehomology_constr=enetb(22,i,iparm)
             edfadis=enetb(23,i,iparm)
             edfator=enetb(24,i,iparm)
             edfanei=enetb(25,i,iparm)
             edfabet=enetb(26,i,iparm)
 #ifdef DEBUG
-            write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
+            write (iout,'(3i5,6f5.2,17f12.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,edihcnstr,
      &       ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
@@ -709,12 +710,16 @@ c              write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
      &      +wbond*estr+ehomology_constr+wdfa_dist*edfadis
      &      +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet
 #endif
-c            write (iout,*) "i",i," ib",ib,
-c     &      " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
-c     &      " entfac",entfac(i)
+#ifdef DEBUG
+            write (iout,*) "i",i," ib",ib,
+     &      " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot,
+     &      " entfac",entfac(i)," ecorr",etot-entfac(i)/beta_h(ib,iparm)
             etot=etot-entfac(i)/beta_h(ib,iparm)
+#endif
             if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
-c            write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
+#ifdef DEBUG
+            write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm)
+#endif
           enddo   ! ib
         enddo     ! iparm
       enddo       ! i
@@ -928,6 +933,10 @@ c      write (iout,*) "me1",me1," scount",scount(me1)
           esccor=enetb(19,t,iparm)
           edihcnstr=enetb(20,t,iparm)
           ehomology_constr=enetb(22,t,iparm)
+          edfadis=enetb(23,t,iparm)
+          edfator=enetb(24,t,iparm)
+          edfanei=enetb(25,t,iparm)
+          edfabet=enetb(26,t,iparm)
           do k=0,nGridT
             betaT=startGridT+k*delta_T
             temper=betaT
@@ -1076,6 +1085,11 @@ c            write (iout,*) ib," PotEmin",potEmin
 #endif
             weight=dexp(-betaT*(etot-potEmin)+entfac(t))
 #ifdef DEBUG
+            write (iout,'(3i5,6f5.2,17f12.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,edihcnstr,
+     &       ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+
+     &       wdfa_nei*edfanei+wdfa_beta*edfabet
             write (iout,*) "iparm",iparm," t",t," temper",temper,
      &        " etot",etot," entfac",entfac(t),
      &        " efree",etot-entfac(t)/betaT," potEmin",potEmin,