projects / unres.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
pdbread & boxshift
[unres.git] / source / unres / src-HCD-5D / energy_p_new-sep_barrier.F
diff --git a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
index c6c6832..96f7777 100644 (file)
--- a/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
+++ b/source/unres/src-HCD-5D/energy_p_new-sep_barrier.F
@@ -85,6 +85,7 @@ c      include 'COMMON.CONTACTS'
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & sigij,r0ij,rcut,sss1,sssgrad1,sqrij
       double precision sscale,sscagrad
+      double precision boxshift
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -97,6 +98,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -106,9 +108,13 @@ cd   &                  'iend=',iend(i,iint)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             rij=xj*xj+yj*yj+zj*zj
             sqrij=dsqrt(rrij)
             eps0ij=eps(itypi,itypj)
@@ -187,6 +193,7 @@ c      include 'COMMON.CONTACTS'
       double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
      & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
       double precision sscale,sscagrad
+      double precision boxshift
 c      write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -199,6 +206,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C Change 12/1/95
         num_conti=0
 C
@@ -210,9 +218,13 @@ cd   &                  'iend=',iend(i,iint)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
 C Change 12/1/95 to calculate four-body interactions
             rij=xj*xj+yj*yj+zj*zj
             sqrij=dsqrt(rij)
@@ -285,6 +297,7 @@ C
      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
       double precision sscale,sscagrad
+      double precision boxshift
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -297,6 +310,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -304,9 +318,13 @@ c        do iint=1,nint_gr(i)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             fac_augm=rrij**expon
             e_augm=augm(itypi,itypj)*fac_augm
@@ -384,6 +402,7 @@ C
      & fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
       logical scheck
       double precision sscale,sscagrad
+      double precision boxshift
 c     print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
 c      do i=iatsc_s,iatsc_e
@@ -396,6 +415,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
 C
 C Calculate SC interaction energy.
 C
@@ -403,9 +423,13 @@ c        do iint=1,nint_gr(i)
 c          do j=istart(i,iint),iend(i,iint)
             itypj=iabs(itype(j))
             if (itypj.eq.ntyp1) cycle
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
             fac_augm=rrij**expon
             e_augm=augm(itypi,itypj)*fac_augm
@@ -478,6 +502,7 @@ C
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision sss1,sssgrad1
       double precision sscale,sscagrad
+      double precision boxshift
 c     double precision rrsave(maxdim)
       logical lprn
       evdw=0.0D0
@@ -499,6 +524,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -523,9 +549,13 @@ c            dscj_inv=dsc_inv(itypj)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -604,6 +634,7 @@ C
       integer itypi,itypj,itypi1,iint,ind,ikont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision sscale,sscagrad
+      double precision boxshift
 c     double precision rrsave(maxdim)
       logical lprn
       evdw=0.0D0
@@ -625,6 +656,7 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -649,9 +681,13 @@ c            dscj_inv=dsc_inv(itypj)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -731,6 +767,7 @@ C
      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
       double precision subchap,sss1,sssgrad1
+      double precision boxshift
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -748,12 +785,8 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -787,79 +820,27 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
             xj=c(1,nres+j)
             yj=c(2,nres+j)
             zj=c(3,nres+j)
-            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
-            if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
-C the energy transfer exist
-              if (zj.lt.buflipbot) then
-C what fraction I am in
-                fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-                sslipj=sscalelip(fracinbuf)
-                ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
-              else if (zi.gt.bufliptop) then
-                fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
-                sslipj=sscalelip(fracinbuf)
-                ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
-             else
-               sslipj=1.0d0
-               ssgradlipj=0.0
-             endif
-           else
-             sslipj=0.0d0
-             ssgradlipj=0.0
-           endif
-           aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &       +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
-           bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
-     &       +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
-           dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-           xj_safe=xj
-           yj_safe=yj
-           zj_safe=zj
-           subchap=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
-                   subchap=1
-                 endif
-               enddo
-             enddo
-           enddo
-           if (subchap.eq.1) then
-             xj=xj_temp-xi
-             yj=yj_temp-yi
-             zj=zj_temp-zi
-           else
-             xj=xj_safe-xi
-             yj=yj_safe-yi
-             zj=zj_safe-zi
-           endif
-           dxj=dc_norm(1,nres+j)
-           dyj=dc_norm(2,nres+j)
-           dzj=dc_norm(3,nres+j)
-           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-           rij=dsqrt(rrij)
-           sss1=sscale(1.0d0/rij,r_cut_int)
-           if (sss1.eq.0.0d0) cycle
-           sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
-           if (sss.lt.1.0d0) then
+            call to_box(xj,yj,zj)
+            call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+            aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &        +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &        +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
+            dxj=dc_norm(1,nres+j)
+            dyj=dc_norm(2,nres+j)
+            dzj=dc_norm(3,nres+j)
+            rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
+            rij=dsqrt(rrij)
+            sss1=sscale(1.0d0/rij,r_cut_int)
+            if (sss1.eq.0.0d0) cycle
+            sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
+            if (sss.lt.1.0d0) then
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
-             sssgrad=
+              sssgrad=
      &         sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
               sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
               call sc_angular
@@ -946,15 +927,13 @@ C
       include 'COMMON.CONTROL'
       include "COMMON.SPLITELE"
       logical lprn
-      integer xshift,yshift,zshift
       double precision evdw
       integer itypi,itypj,itypi1,iint,ind,ikont
       double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
       double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
-     & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
-     & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+     & sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
-      double precision subchap
+      double precision boxshift
       evdw=0.0D0
 ccccc      energy_dec=.false.
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
@@ -972,12 +951,8 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
-          xi=mod(xi,boxxsize)
-          if (xi.lt.0) xi=xi+boxxsize
-          yi=mod(yi,boxysize)
-          if (yi.lt.0) yi=yi+boxysize
-          zi=mod(zi,boxzsize)
-          if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1011,67 +986,15 @@ c            write (iout,*) "i",i," j", j," itype",itype(i),itype(j)
             xj=c(1,nres+j)
             yj=c(2,nres+j)
             zj=c(3,nres+j)
-            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
-            if ((zj.gt.bordlipbot).and.(zj.lt.bordliptop)) then
-C the energy transfer exist
-              if (zj.lt.buflipbot) then
-C what fraction I am in
-                fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
-C lipbufthick is thickenes of lipid buffore
-                sslipj=sscalelip(fracinbuf)
-                ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
-              elseif (zi.gt.bufliptop) then
-                fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
-                sslipj=sscalelip(fracinbuf)
-                ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
-              else
-                sslipj=1.0d0
-                ssgradlipj=0.0
-              endif
-            else
-              sslipj=0.0d0
-              ssgradlipj=0.0
-            endif
+            call to_box(xj,yj,zj)
+            call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
             aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &        +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
             bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
      &        +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
-            dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-            xj_safe=xj
-            yj_safe=yj
-            zj_safe=zj
-            subchap=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
-                    subchap=1
-                  endif
-                enddo
-              enddo
-            enddo
-            if (subchap.eq.1) then
-              xj=xj_temp-xi
-              yj=yj_temp-yi
-              zj=zj_temp-zi
-            else
-              xj=xj_safe-xi
-              yj=yj_safe-yi
-              zj=zj_safe-zi
-            endif
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -1191,6 +1114,8 @@ c      do i=iatsc_s,iatsc_e
         xi=c(1,nres+i)
         yi=c(2,nres+i)
         zi=c(3,nres+i)
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -1217,9 +1142,18 @@ c            dscj_inv=dsc_inv(itypj)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+            aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -1317,6 +1251,7 @@ C
      & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
      & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
       double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+      double precision boxshift
       evdw=0.0D0
 c     print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
       evdw=0.0D0
@@ -1336,6 +1271,8 @@ c      do i=iatsc_s,iatsc_e
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
+        call to_box(xi,yi,zi)
+        call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
 c        dsci_inv=dsc_inv(itypi)
         dsci_inv=vbld_inv(i+nres)
 C
@@ -1359,9 +1296,18 @@ c            dscj_inv=dsc_inv(itypj)
             alf1=alp(itypi)
             alf2=alp(itypj)
             alf12=0.5D0*(alf1+alf2)
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+            xj=c(1,nres+j)
+            yj=c(2,nres+j)
+            zj=c(3,nres+j)
+            call to_box(xj,yj,zj)
+            call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+            aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+     &       +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+            xj=boxshift(xj-xi,boxxsize)
+            yj=boxshift(yj-yi,boxysize)
+            zj=boxshift(zj-zi,boxzsize)
             dxj=dc_norm(1,nres+j)
             dyj=dc_norm(2,nres+j)
             dzj=dc_norm(3,nres+j)
@@ -1612,12 +1558,7 @@ C     &  .or. itype(i+4).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*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)
         num_conti=0
         call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -1641,12 +1582,7 @@ C     &    .or. itype(i-1).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*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)
 #ifdef FOURBODY
         num_conti=num_cont_hb(i)
 #endif
@@ -1677,12 +1613,7 @@ C     &  .or. itype(i-1).eq.ntyp1
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*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',i,' ielstart',ielstart(i),' ielend',ielend
 #ifdef FOURBODY
         num_conti=num_cont_hb(i)
@@ -1764,6 +1695,7 @@ C-------------------------------------------------------------------------------
       double precision sss1,sssgrad1
       double precision sscale,sscagrad
       double precision scalar
+      double precision boxshift
 
 c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
 #ifdef MOMENT
@@ -1796,44 +1728,10 @@ C      print *,"WCHODZE2"
       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)
-      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
-
+      call to_box(xj,yj,zj)
+      xj=boxshift(xj-xmedi,boxxsize)
+      yj=boxshift(yj-ymedi,boxysize)
+      zj=boxshift(zj-zmedi,boxzsize)
       rij=xj*xj+yj*yj+zj*zj
       rrmij=1.0D0/rij
       rij=dsqrt(rij)
@@ -2734,6 +2632,7 @@ c      write (iout,*) "evdwpp_short"
       double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
      & dist_temp, dist_init,sss_grad
       double precision sscale,sscagrad
+      double precision boxshift
       integer ikont
       evdw1=0.0D0
 C      print *,"WCHODZE"
@@ -2754,12 +2653,7 @@ c      do i=iatel_s_vdw,iatel_e_vdw
         xmedi=c(1,i)+0.5d0*dxi
         ymedi=c(2,i)+0.5d0*dyi
         zmedi=c(3,i)+0.5d0*dzi
-        xmedi=mod(xmedi,boxxsize)
-        if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
-        ymedi=mod(ymedi,boxysize)
-        if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
-        zmedi=mod(zmedi,boxzsize)
-        if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
+        call to_box(xmedi,ymedi,zmedi)
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
 c     &   ' ielend',ielend_vdw(i)
@@ -2781,43 +2675,10 @@ c        do j=ielstart_vdw(i),ielend_vdw(i)
           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)
-          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
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xmedi,boxxsize)
+          yj=boxshift(yj-ymedi,boxysize)
+          zj=boxshift(zj-zmedi,boxzsize)
           rij=xj*xj+yj*yj+zj*zj
           rrmij=1.0D0/rij
           rij=dsqrt(rij)
@@ -2881,14 +2742,12 @@ C
       logical lprint_short
       common /shortcheck/ lprint_short
       double precision ggg(3)
-      integer xshift,yshift,zshift
       integer i,iint,j,k,iteli,itypj,subchap
       double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
      & fac,e1,e2,rij
       double precision evdw2,evdw2_14,evdwij
-      double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
-     & dist_temp, dist_init
       double precision sscale,sscagrad
+      double precision boxshift
       integer ikont
       if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
       evdw2=0.0D0
@@ -2907,12 +2766,7 @@ c      do i=iatscp_s,iatscp_e
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-        xi=mod(xi,boxxsize)
-        if (xi.lt.0) xi=xi+boxxsize
-        yi=mod(yi,boxysize)
-        if (yi.lt.0) yi=yi+boxysize
-        zi=mod(zi,boxzsize)
-        if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
 
 c        do iint=1,nscp_gr(i)
 
@@ -2928,44 +2782,10 @@ C Uncomment following three lines for Ca-p interactions
           yj=c(2,j)
           zj=c(3,j)
 c corrected by AL
-          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
-c end correction
-          dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-          xj_safe=xj
-          yj_safe=yj
-          zj_safe=zj
-          subchap=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
-            subchap=1
-          endif
-          enddo
-          enddo
-          enddo
-          if (subchap.eq.1) then
-            xj=xj_temp-xi
-            yj=yj_temp-yi
-            zj=zj_temp-zi
-          else
-            xj=xj_safe-xi
-            yj=yj_safe-yi
-            zj=zj_safe-zi
-          endif
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xi,boxxsize)
+          yj=boxshift(yj-yi,boxysize)
+          zj=boxshift(zj-zi,boxzsize)
 
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 
@@ -3063,6 +2883,7 @@ C
      & dist_temp, dist_init
       double precision ggg(3)
       double precision sscale,sscagrad
+      double precision boxshift
       evdw2=0.0D0
       evdw2_14=0.0d0
 cd    print '(a)','Enter ESCP'
@@ -3076,12 +2897,7 @@ c     & ' iatscp_e=',iatscp_e
         xi=0.5D0*(c(1,i)+c(1,i+1))
         yi=0.5D0*(c(2,i)+c(2,i+1))
         zi=0.5D0*(c(3,i)+c(3,i+1))
-        xi=mod(xi,boxxsize)
-        if (xi.lt.0) xi=xi+boxxsize
-        yi=mod(yi,boxysize)
-        if (yi.lt.0) yi=yi+boxysize
-        zi=mod(zi,boxzsize)
-        if (zi.lt.0) zi=zi+boxzsize
+        call to_box(xi,yi,zi)
 
 c        if (lprint_short) 
 c     &    write (iout,*) "i",i," itype",itype(i),itype(i+1),
@@ -3102,49 +2918,13 @@ C Uncomment following three lines for Ca-p interactions
           yj=c(2,j)
           zj=c(3,j)
 c corrected by AL
-          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
-c end correction
-          dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
-c          if (lprint_short) then
-c            write (iout,*) i,j,xi,yi,zi,xj,yj,zj
-c            write (iout,*) "dist_init",dsqrt(dist_init)
-c          endif
-          xj_safe=xj
-          yj_safe=yj
-          zj_safe=zj
-          subchap=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
-            subchap=1
-          endif
-          enddo
-          enddo
-          enddo
-c          if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
-          if (subchap.eq.1) then
-             xj=xj_temp-xi
-             yj=yj_temp-yi
-             zj=zj_temp-zi
-          else
-             xj=xj_safe-xi
-             yj=yj_safe-yi
-             zj=zj_safe-zi
-          endif
+          xj=c(1,j)
+          yj=c(2,j)
+          zj=c(3,j)
+          call to_box(xj,yj,zj)
+          xj=boxshift(xj-xi,boxxsize)
+          yj=boxshift(yj-yi,boxysize)
+          zj=boxshift(zj-zi,boxzsize)
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
 c          sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
 c          sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
UNRESPACK 3.x