poprawka w triss

[unres.git] / source / wham / src-M / elecont.f

diff --git a/source/wham/src-M/elecont.f b/source/wham/src-M/elecont.f
index 5de56cb..720f860 100644 (file)
--- a/source/wham/src-M/elecont.f
+++ b/source/wham/src-M/elecont.f
@@ -14,11 +14,12 @@
       double precision rri,xi,yi,zi,dxi,dyi,dzi,xmedi,ymedi,zmedi,
      &  xj,yj,zj,dxj,dyj,dzj,aaa,bbb,ael6i,ael3i,rrmij,rmij,r3ij,r6ij,
      &  vrmij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,evdwij,el1,el2,
-     &  eesij,ees,evdw,ene
+     &  eesij,ees,evdw,ene, rij,zj_temp,xj_temp,yj_temp,
+     & sscale,sscagrad,dist_temp,xj_safe,yj_safe,zj_safe,dist_init
       double precision elpp6c(2,2),elpp3c(2,2),ael6c(2,2),ael3c(2,2),
      &  appc(2,2),bppc(2,2)
       double precision elcutoff,elecutoff_14
-      integer ncont,icont(2,maxcont)
+      integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap
       double precision econt(maxcont)
 *
 * Load the constants of peptide bond - peptide bond interactions.
@@ -59,6 +60,12 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
         xmedi=xi+0.5*dxi
         ymedi=yi+0.5*dyi
         zmedi=zi+0.5*dzi
+          xmedi=mod(xmedi,boxxsize)
+          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
         do 4 j=i+2,ien-1
           ind=ind+1
           iteli=itel(i)
@@ -73,9 +80,49 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
           dxj=c(1,j+1)-c(1,j)
           dyj=c(2,j+1)-c(2,j)
           dzj=c(3,j+1)-c(3,j)
-          xj=c(1,j)+0.5*dxj-xmedi
-          yj=c(2,j)+0.5*dyj-ymedi
-          zj=c(3,j)+0.5*dzj-zmedi
+          xj=c(1,j)+0.5*dxj
+          yj=c(2,j)+0.5*dyj
+          zj=c(3,j)+0.5*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)
+          if (zj.lt.0) zj=zj+boxzsize
+      dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**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-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**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
+            sss=sscale(sqrt(rij))
+            sssgrad=sscagrad(sqrt(rij))
           rrmij=1.0/(xj*xj+yj*yj+zj*zj)
           rmij=sqrt(rrmij)
           r3ij=rrmij*rmij
@@ -101,7 +148,7 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
             econt(ncont)=eesij
           endif
           ees=ees+eesij
-          evdw=evdw+evdwij
+          evdw=evdw+evdwij*sss
     4   continue
     1 continue
       if (lprint) then