correction in wham and UNRES for lipid and correlation

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

diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F
index eb6146a..2b99394 100644 (file)
--- a/source/wham/src-M/energy_p_new.F
+++ b/source/wham/src-M/energy_p_new.F
@@ -895,7 +895,7 @@ C
       include 'COMMON.SBRIDGE'
       logical lprn
       common /srutu/icall
-      integer icant
+      integer icant,xshift,yshift,zshift
       external icant
       do i=1,210
         do j=1,2
@@ -1034,6 +1034,8 @@ C lipbufthick is thickenes of lipid buffore
      &  +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
       bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &  +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+C      write(iout,*) "tu,", i,j,aa_lip(itypi,itypj),bb_lip(itypi,itypj)
+
 C        write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
 C checking the distance
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
@@ -2098,15 +2100,15 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         gcorr_loc(i)=0.0d0
       enddo
       do i=iatel_s,iatel_e
-          if (i.eq.1) then 
+C          if (i.eq.1) then 
            if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1) cycle
-          else
-        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
-     &  .or. itype(i+2).eq.ntyp1
-     &  .or. itype(i-1).eq.ntyp1
+C     &  .or. itype(i+2).eq.ntyp1) cycle
+C          else
+C        if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C     &  .or. itype(i+2).eq.ntyp1
+C     &  .or. itype(i-1).eq.ntyp1
      &) cycle
-         endif
+C         endif
         if (itel(i).eq.0) goto 1215
         dxi=dc(1,i)
         dyi=dc(2,i)
@@ -2126,16 +2128,16 @@ C      write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
         num_conti=0
 C        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
-          if (j.eq.1) then
-           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
-     & .or.itype(j+2).eq.ntyp1
-     &) cycle  
-          else     
+          if (j.le.1) cycle
+C           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+C     & .or.itype(j+2).eq.ntyp1
+C     &) cycle  
+C          else     
           if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
-     & .or.itype(j+2).eq.ntyp1
-     & .or.itype(j-1).eq.ntyp1
+C     & .or.itype(j+2).eq.ntyp1
+C     & .or.itype(j-1).eq.ntyp1
      &) cycle
-         endif
+C         endif
 C
 C) cycle
           if (itel(j).eq.0) goto 1216
@@ -3028,6 +3030,18 @@ C Third- and fourth-order contributions from turns
      &    aggj(3,4),aggj1(3,4),a_temp(2,2)
       common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,j1,j2
       if (j.eq.i+2) then
+      if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+C     & .or.((i+5).gt.nres)
+C     & .or.((i-1).le.0)
+C end of changes suggested by Ana
+     &    .or. itype(i+2).eq.ntyp1
+     &    .or. itype(i+3).eq.ntyp1
+C     &    .or. itype(i+5).eq.ntyp1
+C     &    .or. itype(i).eq.ntyp1
+C     &    .or. itype(i-1).eq.ntyp1
+     &    ) goto 178
+
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Third-order contributions
@@ -3158,14 +3172,15 @@ C Cartesian derivatives
       else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
       if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
 C changes suggested by Ana to avoid out of bounds
-     & .or.((i+5).gt.nres)
-     & .or.((i-1).le.0)
+C     & .or.((i+5).gt.nres)
+C     & .or.((i-1).le.0)
 C end of changes suggested by Ana
      &    .or. itype(i+3).eq.ntyp1
      &    .or. itype(i+4).eq.ntyp1
-     &    .or. itype(i+5).eq.ntyp1
+C     &    .or. itype(i+5).eq.ntyp1
      &    .or. itype(i).eq.ntyp1
-     &    .or. itype(i-1).eq.ntyp1) goto 178
+C     &    .or. itype(i-1).eq.ntyp1
+     &    ) goto 178
 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 C
 C               Fourth-order contributions
@@ -3851,6 +3866,7 @@ C     &       gnmr1(vbld(i),-1.0d0,distchainmax)
 C        else
          if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
         diff = vbld(i)-vbldpDUM
+C         write(iout,*) i,diff
          else
           diff = vbld(i)-vbldp0
 c          write (iout,*) i,vbld(i),vbldp0,diff,AKP*diff*diff