From 6f877b0a7b5b4f292dd52b2b447ececb7e469044 Mon Sep 17 00:00:00 2001
From: Pawel Krupa <vetinari@piasek4.chem.univ.gda.pl>
Date: Thu, 16 Apr 2015 22:58:27 +0200
Subject: [PATCH] Fixed energy differs in WHAM-M (EBE).

---
 source/wham/src-M/energy_p_new.F |   31 ++++++++++++++++++-------------
 1 file changed, 18 insertions(+), 13 deletions(-)

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index 8f932cc..cac8518 100644
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -3179,8 +3179,8 @@ C
      & delthe0,sig0inv,sigtc,sigsqtc,delthec,it
       double precision y(2),z(2)
       delta=0.02d0*pi
-      time11=dexp(-2*time)
-      time12=1.0d0
+c      time11=dexp(-2*time)
+c      time12=1.0d0
       etheta=0.0D0
 c      write (iout,*) "nres",nres
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
@@ -3206,8 +3206,8 @@ C Zero the energy function and its derivative at 0 or pi.
         if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
-          icrc=0
-          call proc_proc(phii,icrc)
+c          icrc=0
+c          call proc_proc(phii,icrc)
           if (icrc.eq.1) phii=150.0
 #else
           phii=phi(i)
@@ -3221,8 +3221,8 @@ C Zero the energy function and its derivative at 0 or pi.
         if (i.lt.nres .and. itype(i).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
-          icrc=0
-          call proc_proc(phii1,icrc)
+c          icrc=0
+c          call proc_proc(phii1,icrc)
           if (icrc.eq.1) phii1=150.0
           phii1=pinorm(phii1)
           z(1)=cos(phii1)
@@ -3290,7 +3290,7 @@ c     &    rad2deg*phii,rad2deg*phii1,ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2
         gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett)
- 1215   continue
+c 1215   continue
       enddo
 C Ufff.... We've done all this!!! 
       return
@@ -3433,7 +3433,9 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-        if (itype(i-1).eq.ntyp1) cycle
+c        if (itype(i-1).eq.ntyp1) cycle
+        if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or.
+     &(itype(i).eq.ntyp1)) cycle
         if (iabs(itype(i+1)).eq.20) iblock=2
         if (iabs(itype(i+1)).ne.20) iblock=1
         dethetai=0.0d0
@@ -3445,7 +3447,7 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-2).ne.ntyp1) then
+        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -3459,13 +3461,14 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+c          ityp1=nthetyp+1
           do k=1,nsingle
+            ityp1=ithetyp((itype(i-2)))
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
           enddo 
         endif
-        if (i.lt.nres .and. itype(i).ne.ntyp1) then
+        if (i.lt.nres .and. itype(i+1).ne.ntyp1) then
 #ifdef OSF
           phii1=phi(i+1)
           if (phii1.ne.phii1) phii1=150.0
@@ -3480,7 +3483,8 @@ c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+c          ityp3=nthetyp+1
+          ityp3=ithetyp((itype(i)))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0
@@ -3597,7 +3601,8 @@ c        call flush(iout)
         etheta=etheta+ethetai
         if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii
         if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1
-        gloc(nphi+i-2,icg)=wang*dethetai
+c        gloc(nphi+i-2,icg)=wang*dethetai
+        gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai
       enddo
       return
       end
-- 
1.7.9.5