readpdb correction

[unres.git] / source / unres / src_MD-M / ssMD.F

diff --git a/source/unres/src_MD-M/ssMD.F b/source/unres/src_MD-M/ssMD.F
index de0abb0..51ad513 100644 (file)
--- a/source/unres/src_MD-M/ssMD.F
+++ b/source/unres/src_MD-M/ssMD.F
@@ -2012,22 +2012,22 @@ C The first case the ith atom is the center
 C Energy function is E=d/(a*(x-y)**2+b*(x+y)**2+c) where x is first
 C distance y is second distance the a,b,c,d are parameters derived for
 C this problem d parameter was set as a penalty currenlty set to 1.
-      if ((iabs(j-i).eq.1).or.(iabs(i-k).eq.1)) then
+      if ((iabs(j-i).le.2).or.(iabs(i-k).le.2)) then
       eij1=0.0d0
       else
-      eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
+      eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**6+ctriss)
       endif
 C second case jth atom is center
-      if ((iabs(j-i).eq.1).or.(iabs(j-k).eq.1)) then
+      if ((iabs(j-i).le.2).or.(iabs(j-k).le.2)) then
       eij2=0.0d0
       else
-      eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
+      eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**6+ctriss)
       endif
 C the third case kth atom is the center
-      if ((iabs(i-k).eq.1).or.(iabs(j-k).eq.1)) then
+      if ((iabs(i-k).le.2).or.(iabs(j-k).le.2)) then
       eij3=0.0d0
       else
-      eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
+      eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**6+ctriss)
       endif
 C      eij2=0.0
 C      eij3=0.0
@@ -2036,8 +2036,8 @@ C      eij1=0.0
 C      write(iout,*)i,j,k,eij
 C The energy penalty calculated now time for the gradient part 
 C derivative over rij
-      fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
-     &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))  
+      fac=-eij1**2/dtriss*(2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5)
+     &-eij2**2/dtriss*(2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)  
             gg(1)=xij*fac/rij
             gg(2)=yij*fac/rij
             gg(3)=zij*fac/rij
@@ -2050,8 +2050,9 @@ C derivative over rij
         gvdwc(l,j)=gvdwc(l,j)+gg(l)
       enddo
 C now derivative over rik
-      fac=-eij1**2/dtriss*(-2.0*atriss*(rij-rik)+2.0*btriss*(rij+rik))
-     &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
+      fac=-eij1**2/dtriss*
+     &(-2.0*atriss*(rij-rik)+6.0*btriss*(rij+rik)**5)
+     &-eij3**2/dtriss*(2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
             gg(1)=xik*fac/rik
             gg(2)=yik*fac/rik
             gg(3)=zik*fac/rik
@@ -2064,8 +2065,9 @@ C now derivative over rik
         gvdwc(l,k)=gvdwc(l,k)+gg(l)
       enddo
 C now derivative over rjk
-      fac=-eij2**2/dtriss*(-2.0*atriss*(rij-rjk)+2.0*btriss*(rij+rjk))-
-     &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+2.0*btriss*(rik+rjk))
+      fac=-eij2**2/dtriss*
+     &(-2.0*atriss*(rij-rjk)+6.0*btriss*(rij+rjk)**5)-
+     &eij3**2/dtriss*(-2.0*atriss*(rik-rjk)+6.0*btriss*(rik+rjk)**5)
             gg(1)=xjk*fac/rjk
             gg(2)=yjk*fac/rjk
             gg(3)=zjk*fac/rjk