projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
poprawka w make (usuniecie nss)
[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 30fc811..3872257 100644 (file)
--- a/source/unres/src_MD-M/ssMD.F
+++ b/source/unres/src_MD-M/ssMD.F
@@ -187,9 +187,10 @@ c      om12=dxi*dxj+dyi*dyj+dzi*dzj
 
       ljXs=sig-sig0ij
       ljA=eps1*eps2rt**2*eps3rt**2
-      ljB=ljA*bb(itypi,itypj)
-      ljA=ljA*aa(itypi,itypj)
-      ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+      ljB=ljA*bb_aq(itypi,itypj)
+      ljA=ljA*aa_aq(itypi,itypj)
+      ljxm=ljXs+(-2.0D0*aa_aq(itypi,itypj)/
+     & bb_aq(itypi,itypj))**(1.0D0/6.0D0)
 
       ssXs=d0cm
       deltat1=1.0d0-om1
@@ -223,7 +224,7 @@ c-------TESTING CODE
 c     Stop and plot energy and derivative as a function of distance
       if (checkstop) then
         ssm=ssC-0.25D0*ssB*ssB/ssA
-        ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+        ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)
         if (ssm.lt.ljm .and.
      &       dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
           nicheck=1000
@@ -248,8 +249,8 @@ c-------END TESTING CODE
         havebond=.false.
         ljd=rij-ljXs
         fac=(1.0D0/ljd)**expon
-        e1=fac*fac*aa(itypi,itypj)
-        e2=fac*bb(itypi,itypj)
+        e1=fac*fac*aa_aq(itypi,itypj)
+        e2=fac*bb_aq(itypi,itypj)
         eij=eps1*eps2rt*eps3rt*(e1+e2)
         eps2der=eij*eps3rt
         eps3der=eij*eps2rt
@@ -314,8 +315,8 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
           eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
         else
           havebond=.false.
-          ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
-          d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
+          ljm=-0.25D0*ljB*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)
+          d_ljm(1)=-0.5D0*bb_aq(itypi,itypj)/aa_aq(itypi,itypj)*ljB
           d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
           d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
      +         alf12/eps3rt)
@@ -1959,6 +1960,7 @@ c      integer itypi,itypj,k,l
       double precision ssA,ssB,ssC,ssXs
       double precision ssxm,ljxm,ssm,ljm
       double precision d_ssxm(1:3),d_ljxm(1:3),d_ssm(1:3),d_ljm(1:3)
+      eij=0.0
       if (dtriss.eq.0) return
       i=resi
       j=resj
@@ -2011,11 +2013,23 @@ 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.
-      eij1=dtriss/(atriss*(rij-rik)**2+btriss*(rij+rik)**2+ctriss)
+      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)**6+ctriss)
+      endif
 C second case jth atom is center
-      eij2=dtriss/(atriss*(rij-rjk)**2+btriss*(rij+rjk)**2+ctriss)
+      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)**6+ctriss)
+      endif
 C the third case kth atom is the center
-      eij3=dtriss/(atriss*(rik-rjk)**2+btriss*(rik+rjk)**2+ctriss)
+      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)**6+ctriss)
+      endif
 C      eij2=0.0
 C      eij3=0.0
 C      eij1=0.0
@@ -2023,8 +2037,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
@@ -2037,8 +2051,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
@@ -2051,8 +2066,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
UNRESPACK 3.x