X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new-sep.F;h=a7f0bab4843b7cca557d03b4df0c09bcd587fd4e;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=0b8f27b07646e94854c1853316a020b691decf2c;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git

diff --git a/source/unres/src_MD-M/energy_p_new-sep.F b/source/unres/src_MD-M/energy_p_new-sep.F
index 0b8f27b..a7f0bab 100644
--- a/source/unres/src_MD-M/energy_p_new-sep.F
+++ b/source/unres/src_MD-M/energy_p_new-sep.F
@@ -2,6 +2,7 @@ C-----------------------------------------------------------------------
       double precision function sscale(r)
       double precision r,gamm
       include "COMMON.SPLITELE"
+      include "COMMON.CHAIN"
       if(r.lt.r_cut-rlamb) then
         sscale=1.0d0
       else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
@@ -13,6 +14,23 @@ C-----------------------------------------------------------------------
       return
       end
 C-----------------------------------------------------------------------
+C-----------------------------------------------------------------------
+      double precision function sscagrad(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+      include "COMMON.CHAIN"
+      if(r.lt.r_cut-rlamb) then
+        sscagrad=0.0d0
+      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+        gamm=(r-(r_cut-rlamb))/rlamb
+        sscagrad=gamm*(6*gamm-6.0d0)/rlamb
+      else
+        sscagrad=0.0d0
+      endif
+      return
+      end
+C-----------------------------------------------------------------------
+
       subroutine elj_long(evdw)
 C
 C This subroutine calculates the interaction energy of nonbonded side chains
@@ -2110,6 +2128,12 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=mod(xmedi,boxysize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=mod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=mod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
@@ -2125,13 +2149,50 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
           dx_normj=dc_norm(1,j)
           dy_normj=dc_norm(2,j)
           dz_normj=dc_norm(3,j)
-          xj=c(1,j)+0.5D0*dxj-xmedi
-          yj=c(2,j)+0.5D0*dyj-ymedi
-          zj=c(3,j)+0.5D0*dzj-zmedi
-          rij=xj*xj+yj*yj+zj*zj
+          xj=c(1,j)+0.5D0*dxj
+          yj=c(2,j)+0.5D0*dyj
+          zj=c(3,j)+0.5D0*dzj
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+      dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      isubchap=0
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            isubchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (isubchap.eq.1) then
+          xj=xj_temp-xmedi
+          yj=yj_temp-ymedi
+          zj=zj_temp-zmedi
+       else
+          xj=xj_safe-xmedi
+          yj=yj_safe-ymedi
+          zj=zj_safe-zmedi
+       endif
+         rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
           sss=sscale(rij/rpp(iteli,itelj))
+            sssgrad=sscagrad(sqrt(rij))
           if (sss.lt.1.0d0) then
             rmij=1.0D0/rij
             r3ij=rrmij*rmij
@@ -2150,9 +2211,9 @@ C
 C Calculate contributions to the Cartesian gradient.
 C
             facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
-            ggg(1)=facvdw*xj
-            ggg(2)=facvdw*yj
-            ggg(3)=facvdw*zj
+            ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj
+            ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj
+            ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj
 
             do k=1,3
               ghalf=0.5D0*ggg(k)
@@ -2198,6 +2259,7 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
 #else
       double precision scal_el /0.5d0/
 #endif
+c      write (iout,*) "evdwpp_short"
       evdw1=0.0D0
       do i=iatel_s,iatel_e
         dxi=dc(1,i)
@@ -2209,6 +2271,12 @@ c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
+          xmedi=mod(xmedi,boxysize)
+          if (xmedi.lt.0) xmedi=xmedi+boxxsize
+          ymedi=mod(ymedi,boxysize)
+          if (ymedi.lt.0) ymedi=ymedi+boxysize
+          zmedi=mod(zmedi,boxzsize)
+          if (zmedi.lt.0) zmedi=zmedi+boxzsize
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
@@ -2224,13 +2292,50 @@ c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
           dx_normj=dc_norm(1,j)
           dy_normj=dc_norm(2,j)
           dz_normj=dc_norm(3,j)
-          xj=c(1,j)+0.5D0*dxj-xmedi
-          yj=c(2,j)+0.5D0*dyj-ymedi
-          zj=c(3,j)+0.5D0*dzj-zmedi
+          xj=c(1,j)+0.5D0*dxj
+          yj=c(2,j)+0.5D0*dyj
+          zj=c(3,j)+0.5D0*dzj
+          xj=mod(xj,boxxsize)
+          if (xj.lt.0) xj=xj+boxxsize
+          yj=mod(yj,boxysize)
+          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+     dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+      xj_safe=xj
+      yj_safe=yj
+      zj_safe=zj
+      isubchap=0
+      do xshift=-1,1
+      do yshift=-1,1
+      do zshift=-1,1
+          xj=xj_safe+xshift*boxxsize
+          yj=yj_safe+yshift*boxysize
+          zj=zj_safe+zshift*boxzsize
+          dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+          if(dist_temp.lt.dist_init) then
+            dist_init=dist_temp
+            xj_temp=xj
+            yj_temp=yj
+            zj_temp=zj
+            isubchap=1
+          endif
+       enddo
+       enddo
+       enddo
+       if (isubchap.eq.1) then
+          xj=xj_temp-xmedi
+          yj=yj_temp-ymedi
+          zj=zj_temp-zmedi
+       else
+          xj=xj_safe-xmedi
+          yj=yj_safe-ymedi
+          zj=zj_safe-zmedi
+       endif
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
           sss=sscale(rij/rpp(iteli,itelj))
+            sssgrad=sscagrad(sqrt(rij))
           if (sss.gt.0.0d0) then
             rmij=1.0D0/rij
             r3ij=rrmij*rmij
@@ -2249,9 +2354,9 @@ C
 C Calculate contributions to the Cartesian gradient.
 C
             facvdw=-6*rrmij*(ev1+evdwij)*sss
-            ggg(1)=facvdw*xj
-            ggg(2)=facvdw*yj
-            ggg(3)=facvdw*zj
+            ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+            ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+            ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
 
             do k=1,3
               ghalf=0.5D0*ggg(k)
@@ -2289,9 +2394,11 @@ C
       include 'COMMON.FFIELD'
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
       dimension ggg(3)
       evdw2=0.0D0
       evdw2_14=0.0d0
+      if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
 cd    print '(a)','Enter ESCP'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
       do i=iatscp_s,iatscp_e
@@ -2405,11 +2512,13 @@ C
       include 'COMMON.FFIELD'
       include 'COMMON.IOUNITS'
       include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
       dimension ggg(3)
       evdw2=0.0D0
       evdw2_14=0.0d0
 cd    print '(a)','Enter ESCP'
 cd    write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
+      if (energy_dec) write (iout,*) "escp_short:",r_cut,rlamb
       do i=iatscp_s,iatscp_e
         iteli=itel(i)
         xi=0.5D0*(c(1,i)+c(1,i+1))