X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fsc_move.F;fp=source%2Funres%2Fsrc_MD-M%2Fsc_move.F;h=c552ee06068b36564aa2c6568b1040f3fad88897;hb=8df996bdc4bc56846b3799e9a2d47c24a9273125;hp=274767b11386a7e30b9159e119d545d6762522fa;hpb=ad2687da4ea0c0ddde6c400a60c4817a689b3dff;p=unres.git

diff --git a/source/unres/src_MD-M/sc_move.F b/source/unres/src_MD-M/sc_move.F
index 274767b..c552ee0 100644
--- a/source/unres/src_MD-M/sc_move.F
+++ b/source/unres/src_MD-M/sc_move.F
@@ -724,6 +724,35 @@ c     if (icall.eq.0) lprn=.true.
         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
+       if ((zi.gt.bordlipbot)
+     &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+        if (zi.lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zi.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipi=1.0d0
+         ssgradlipi=0.0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0
+       endif
+
         dxi=dc_norm(1,nres+i)
         dyi=dc_norm(2,nres+i)
         dzi=dc_norm(3,nres+i)
@@ -757,14 +786,85 @@ c           chip12=0.0D0
 c           alf1=0.0D0
 c           alf2=0.0D0
 c           alf12=0.0D0
-            xj=c(1,nres+j)-xi
-            yj=c(2,nres+j)-yi
-            zj=c(3,nres+j)-zi
+C            xj=c(1,nres+j)-xi
+C            yj=c(2,nres+j)-yi
+C            zj=c(3,nres+j)-zi
+            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-
+     &        ((positi-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-positi)/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)
+            sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
+            sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+
 C Calculate angle-dependent terms of energy and contributions to their
 C derivatives.
             call sc_angular
@@ -783,16 +883,16 @@ cd     &        rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq)
 c---------------------------------------------------------------
             rij_shift=1.0D0/rij_shift 
             fac=rij_shift**expon
-            e1=fac*fac*aa(itypi,itypj)
-            e2=fac*bb(itypi,itypj)
+            e1=fac*fac*aa
+            e2=fac*bb
             evdwij=eps1*eps2rt*eps3rt*(e1+e2)
             eps2der=evdwij*eps3rt
             eps3der=evdwij*eps2rt
             evdwij=evdwij*eps2rt*eps3rt
             evdw=evdw+evdwij
             if (lprn) then
-            sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
-            epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+            sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+            epsi=bb**2/aa
 cd            write (iout,'(2(a3,i3,2x),17(0pf7.3))')
 cd     &        restyp(itypi),i,restyp(itypj),j,
 cd     &        epsi,sigm,chi1,chi2,chip1,chip2,
@@ -809,10 +909,13 @@ C Calculate gradient components.
             fac=-expon*(e1+evdwij)*rij_shift
             sigder=fac*sigder
             fac=rij*fac
+            fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
 C Calculate the radial part of the gradient
             gg(1)=xj*fac
             gg(2)=yj*fac
             gg(3)=zj*fac
+            gg_lipi(3)=ssgradlipi*evdwij
+            gg_lipj(3)=ssgradlipj*evdwij
 C Calculate angular part of the gradient.
             call sc_grad
           ENDIF