corrected cluster for lipid

diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F
index e7759b2..206d327 100644 (file)
--- a/source/cluster/wham/src-M/energy_p_new.F
+++ b/source/cluster/wham/src-M/energy_p_new.F
@@ -2237,6 +2237,31 @@ C         endif
           if (ymedi.lt.0) ymedi=ymedi+boxysize
           zmedi=mod(zmedi,boxzsize)
           if (zmedi.lt.0) zmedi=zmedi+boxzsize
+          zmedi2=mod(zmedi,boxzsize)
+          if (zmedi2.lt.0) zmedi2=zmedi2+boxzsize
+       if ((zmedi2.gt.bordlipbot)
+     &.and.(zmedi2.lt.bordliptop)) then
+C the energy transfer exist
+        if (zmedi2.lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((zmedi2-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zmedi2.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-zmedi2)/lipbufthick)
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipi=1.0d0
+         ssgradlipi=0.0d0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0d0
+       endif
+
         num_conti=0
 c        write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
         do j=ielstart(i),ielend(i)
@@ -2280,6 +2305,29 @@ C End diagnostics
           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
+
       dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
       xj_safe=xj
       yj_safe=yj
@@ -2349,12 +2397,13 @@ C#undef DEBUG
           el1=el1*fac_shield(i)**2*fac_shield(j)**2
           el2=el2*fac_shield(i)**2*fac_shield(j)**2
           eesij=(el1+el2)
-     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+C     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           ees=ees+eesij
           else
           fac_shield(i)=1.0
           fac_shield(j)=1.0
           eesij=(el1+el2)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           ees=ees+eesij
           endif
 C          ees=ees+eesij