Adam's wham update

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

diff --git a/source/wham/src-HCD/elecont.f b/source/wham/src-HCD/elecont.f
index fb105a4..86db2df 100644 (file)
--- a/source/wham/src-HCD/elecont.f
+++ b/source/wham/src-HCD/elecont.f
@@ -18,10 +18,11 @@
      &  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)
+     &  appc(2,2),bppc(2,2),epp_(2,2),rpp_(2,2)
       double precision elcutoff,elecutoff_14
       integer ncont,icont(2,maxcont),xshift,yshift,zshift,isubchap
       double precision econt(maxcont)
+      double precision boxshift
 *
 * Load the constants of peptide bond - peptide bond interactions.
 * Type 1 - ordinary peptide bond, type 2 - alkylated peptide bond (e.g.
@@ -29,8 +30,8 @@
 *
 * as of 7/06/91.
 *
-c      data epp    / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
-c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
+c      data epp_   / 0.3045d0, 0.3649d0, 0.3649d0, 0.5743d0/
+      data rpp_   / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
       data elpp6c  /-0.2379d0,-0.2056d0,-0.2056d0,-0.0610d0/
       data elpp3c  / 0.0503d0, 0.0000d0, 0.0000d0, 0.0692d0/
       data elcutoff /-0.3d0/,elecutoff_14 /-0.5d0/
@@ -40,7 +41,7 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
      &  "Constants of electrostatic interaction energy expression."
       do i=1,2
         do j=1,2
-        rri=rpp(i,j)**6
+        rri=rpp_(i,j)**6
         appc(i,j)=epp(i,j)*rri*rri 
         bppc(i,j)=-2.0*epp(i,j)*rri
         ael6c(i,j)=elpp6c(i,j)*4.2**6
@@ -62,12 +63,8 @@ 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
+        call to_box(xmedi,ymedi,zmedi)
+c        write (iout,*) "i",xmedi,ymedi,zmedi
         do 4 j=i+2,ien-1
           jj=iperm(j,ipermmin)
           ind=ind+1
@@ -86,46 +83,13 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
           xj=c(1,jj)+0.5*dxj
           yj=c(2,jj)+0.5*dyj
           zj=c(3,jj)+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
+c          write (iout,*) "j",xj,yj,zj
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xmedi,boxxsize)
+          yj=boxshift(yj-ymedi,boxysize)
+          zj=boxshift(zj-zmedi,boxzsize)
+c          write (iout,*) "j",xj,yj,zj
           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
@@ -152,6 +116,7 @@ c      data rpp    / 4.5088d0, 4.5395d0, 4.5395d0, 4.4846d0/
           endif
           ees=ees+eesij
           evdw=evdw+evdwij*sss
+c          write (iout,*) "i",i," j",j," rij",dsqrt(rij)," eesij",eesij
     4   continue
     1 continue
       if (lprint) then