X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2FssMD.F;h=f298c98617dc5e14ef304513f8dcfbcab916d532;hb=14b55af23c7da34bfdc10f315551c8e01a3a906a;hp=283adf3218abc79d91803e7c4d3d5897d550b8f6;hpb=13279fc00959f1e7a4bf024637996dcf4c729e99;p=unres.git

diff --git a/source/wham/src-M/ssMD.F b/source/wham/src-M/ssMD.F
index 283adf3..f298c98 100644
--- a/source/wham/src-M/ssMD.F
+++ b/source/wham/src-M/ssMD.F
@@ -150,11 +150,76 @@ c-------END TESTING CODE
       dyi=dc_norm(2,nres+i)
       dzi=dc_norm(3,nres+i)
       dsci_inv=vbld_inv(i+nres)
-
+      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-
+     &        ((zi-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-zi)/lipbufthick)
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipi=1.0d0
+         ssgradlipi=0.0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0
+       endif
       itypj=itype(j)
-      xj=c(1,nres+j)-c(1,nres+i)
-      yj=c(2,nres+j)-c(2,nres+i)
-      zj=c(3,nres+j)-c(3,nres+i)
+      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-
+     &        ((zj-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipj=sscalelip(fracinbuf)
+         ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zj.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-zj)/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
+      xj=xj-xi
+      yj=yj-yi
+      zj=zj-zi
       dxj=dc_norm(1,nres+j)
       dyj=dc_norm(2,nres+j)
       dzj=dc_norm(3,nres+j)
@@ -187,9 +252,9 @@ c      om12=dxi*dxj+dyi*dyj+dzi*dzj
 
       ljXs=sig-sig0ij
       ljA=eps1*eps2rt**2*eps3rt**2
-      ljB=ljA*bb(itypi,itypj)
-      ljA=ljA*aa(itypi,itypj)
-      ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+      ljB=ljA*bb
+      ljA=ljA*aa
+      ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
 
       ssXs=d0cm
       deltat1=1.0d0-om1
@@ -223,7 +288,7 @@ c-------TESTING CODE
 c     Stop and plot energy and derivative as a function of distance
       if (checkstop) then
         ssm=ssC-0.25D0*ssB*ssB/ssA
-        ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+        ljm=-0.25D0*ljB*bb/aa
         if (ssm.lt.ljm .and.
      &       dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
           nicheck=1000
@@ -248,9 +313,10 @@ c-------END TESTING CODE
         havebond=.false.
         ljd=rij-ljXs
         fac=(1.0D0/ljd)**expon
-        e1=fac*fac*aa(itypi,itypj)
-        e2=fac*bb(itypi,itypj)
+        e1=fac*fac*aa
+        e2=fac*bb
         eij=eps1*eps2rt*eps3rt*(e1+e2)
+C        write(iout,*) eij,'TU?1'
         eps2der=eij*eps3rt
         eps3der=eij*eps2rt
         eij=eij*eps2rt*eps3rt
@@ -267,7 +333,7 @@ c-------END TESTING CODE
         havebond=.true.
         ssd=rij-ssXs
         eij=ssA*ssd*ssd+ssB*ssd+ssC
-
+C        write(iout,*) 'TU?2',ssc,ssd
         ed=2*akcm*ssd+akct*deltat12
         pom1=akct*ssd
         pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
@@ -303,6 +369,7 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
           h1=h_base(f1,hd1)
           h2=h_base(f2,hd2)
           eij=ssm*h1+Ht*h2
+C         write(iout,*) eij,'TU?3'
           delta_inv=1.0d0/(xm-ssxm)
           deltasq_inv=delta_inv*delta_inv
           fac=ssm*hd1-Ht*hd2
@@ -314,8 +381,8 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
           eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
         else
           havebond=.false.
-          ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
-          d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
+          ljm=-0.25D0*ljB*bb/aa
+          d_ljm(1)=-0.5D0*bb/aa*ljB
           d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
           d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
      +         alf12/eps3rt)
@@ -325,6 +392,7 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE
           h1=h_base(f1,hd1)
           h2=h_base(f2,hd2)
           eij=Ht*h1+ljm*h2
+C        write(iout,*) 'TU?4',ssA
           delta_inv=1.0d0/(ljxm-xm)
           deltasq_inv=delta_inv*delta_inv
           fac=Ht*hd1-ljm*hd2
@@ -396,7 +464,7 @@ c$$$        if (ed.gt.0.0d0) havebond=.true.
 c-------END SECOND METHOD, CONTINUOUS SECOND DERIVATIVE
 
       endif
-
+      write(iout,*) 'havebond',havebond
       if (havebond) then
 #ifndef CLUST
 #ifndef WHAM