cluster wham src-M corrections to match unres energy

diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index 4e4a386..2302af9 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -2197,12 +2197,12 @@ cd      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
-           if (i.le.1) cycle
+cAna           if (i.le.1) cycle
            if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
-     &  .or. ((i+2).gt.nres)
-     &  .or. ((i-1).le.0)
-     &  .or. itype(i+2).eq.ntyp1
-     &  .or. itype(i-1).eq.ntyp1
+cAna     &  .or. ((i+2).gt.nres)
+cAna     &  .or. ((i-1).le.0)
+cAna     &  .or. itype(i+2).eq.ntyp1
+cAna     &  .or. itype(i-1).eq.ntyp1
      &) cycle
 C         endif
         if (itel(i).eq.0) goto 1215
@@ -2224,12 +2224,12 @@ C         endif
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
-          if (j.le.1) cycle
+cAna          if (j.le.1) cycle
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
-     & .or.((j+2).gt.nres)
-     & .or.((j-1).le.0)
-     & .or.itype(j+2).eq.ntyp1
-     & .or.itype(j-1).eq.ntyp1
+cAna     & .or.((j+2).gt.nres)
+cAna     & .or.((j-1).le.0)
+cAna     & .or.itype(j+2).eq.ntyp1
+cAna     & .or.itype(j-1).eq.ntyp1
      &) cycle
 C         endif
           if (itel(j).eq.0) goto 1216
@@ -4671,7 +4671,8 @@ c      write (iout,*) "nres",nres
 c     write (*,'(a,i2)') 'EBEND ICG=',icg
 c      write (iout,*) ithet_start,ithet_end
       do i=ithet_start,ithet_end
-        if (i.le.2) cycle
+C        if (itype(i-1).eq.ntyp1) cycle
+c        if (i.le.2) cycle
         if ((itype(i-1).eq.ntyp1).or.itype(i-2).eq.ntyp1
      &  .or.itype(i).eq.ntyp1) cycle
 C Zero the energy function and its derivative at 0 or pi.
@@ -4693,7 +4694,8 @@ C Zero the energy function and its derivative at 0 or pi.
           y(1)=0.0D0
           y(2)=0.0D0
           else
-        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+
+        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
 c          icrc=0
@@ -4953,7 +4955,7 @@ C
       etheta=0.0D0
 c      write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1)
       do i=ithet_start,ithet_end
-        if (i.eq.2) cycle
+c        if (i.eq.2) cycle
 c        print *,i,itype(i-1),itype(i),itype(i-2)
         if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1)
      &  .or.(itype(i).eq.ntyp1)) cycle
@@ -4970,7 +4972,7 @@ C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
           coskt(k)=dcos(k*theti2)
           sinkt(k)=dsin(k*theti2)
         enddo
-        if (i.gt.3 .and. itype(i-3).ne.ntyp1) then
+        if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then
 #ifdef OSF
           phii=phi(i)
           if (phii.ne.phii) phii=150.0
@@ -4984,7 +4986,7 @@ C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
           enddo
         else
           phii=0.0d0
-          ityp1=nthetyp+1
+          ityp1=ithetyp(itype(i-2))
           do k=1,nsingle
             cosph1(k)=0.0d0
             sinph1(k)=0.0d0
@@ -5005,7 +5007,7 @@ C In current verion the ALL DUMMY ATOM POTENTIALS ARE OFF
           enddo
         else
           phii1=0.0d0
-          ityp3=nthetyp+1
+          ityp3=ithetyp(itype(i))
           do k=1,nsingle
             cosph2(k)=0.0d0
             sinph2(k)=0.0d0