Adams corrections

[unres.git] / source / cluster / wham / src-HCD / energy_p_new.F

diff --git a/source/cluster/wham/src-HCD/energy_p_new.F b/source/cluster/wham/src-HCD/energy_p_new.F
index 5cc851c..27e944b 100644 (file)
--- a/source/cluster/wham/src-HCD/energy_p_new.F
+++ b/source/cluster/wham/src-HCD/energy_p_new.F
@@ -1775,13 +1775,17 @@ C to calculate the el-loc multibody terms of various order.
 C
 c      write(iout,*) 'SET_MATRICES nphi=',nphi,nres
       do i=3,nres+1
-        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+        ii=ireschain(i-2)
+        if (ii.eq.0) cycle
+        innt=chain_border(1,ii)
+        inct=chain_border(2,ii)
+        if (i.gt. innt+2 .and. i.lt.inct+2) then
           iti = itype2loc(itype(i-2))
         else
           iti=nloctyp
         endif
 c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
-        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+        if (i.gt. innt+1 .and. i.lt.inct+1) then
           iti1 = itype2loc(itype(i-1))
         else
           iti1=nloctyp
@@ -1856,6 +1860,19 @@ c        b2tilde(2,i-2)=-b2(2,i-2)
         write (iout,*) 'theta=', theta(i-1)
 #endif
 #else
+        if (i.gt. innt+2 .and. i.lt.inct+2) then
+c        if (i.gt. nnt+2 .and. i.lt.nct+2) then
+          iti = itype2loc(itype(i-2))
+        else
+          iti=nloctyp
+        endif
+c        write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti
+c        if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then
+        if (i.gt. nnt+1 .and. i.lt.nct+1) then
+          iti1 = itype2loc(itype(i-1))
+        else
+          iti1=nloctyp
+        endif
 c        if (i.gt. nnt+2 .and. i.lt.nct+2) then
 c          iti = itype2loc(itype(i-2))
 c        else
@@ -4520,6 +4537,10 @@ c
       estr1=0.0d0
 c      write (iout,*) "distchainmax",distchainmax
       do i=nnt+1,nct
+#ifdef FIVEDIAG
+        if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) cycle
+        diff = vbld(i)-vbldp0
+#else
         if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
 C          estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
 C          do j=1,3
@@ -4537,6 +4558,9 @@ C         write(iout,*) i,diff
           diff = vbld(i)-vbldp0
 c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff
          endif
+#endif
+        if (energy_dec) write (iout,'(a7,i5,4f7.3)')
+     &     "estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
           estr=estr+diff*diff
           do j=1,3
             gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
@@ -4555,8 +4579,9 @@ c
           nbi=nbondterm(iti)
           if (nbi.eq.1) then
             diff=vbld(i+nres)-vbldsc0(1,iti)
-C            write (iout,*) i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
-C     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
+            if (energy_dec) write (iout,*) 
+     &      i,iti,vbld(i+nres),vbldsc0(1,iti),diff,
+     &      AKSC(1,iti),AKSC(1,iti)*diff*diff
             estr=estr+0.5d0*AKSC(1,iti)*diff*diff
             do j=1,3
               gradbx(j,i)=AKSC(1,iti)*diff*dc(j,i+nres)/vbld(i+nres)