bugfix in shield FGPROC>1

diff --git a/source/unres/src_MD-M/COMMON.SHIELD b/source/unres/src_MD-M/COMMON.SHIELD
index 6cbf5a8..2ead23e 100644 (file)
--- a/source/unres/src_MD-M/COMMON.SHIELD
+++ b/source/unres/src_MD-M/COMMON.SHIELD
@@ -1,13 +1,15 @@
        double precision VSolvSphere,VSolvSphere_div,long_r_sidechain,
      & short_r_sidechain,fac_shield,grad_shield_side,grad_shield,
-     & buff_shield,wshield            
+     & buff_shield,wshield,lipscale,sslipi,sslipj,
+     & ssgradlipi,ssgradlipj
        integer  ishield_list,shield_list,ees0plist
        common /shield/ VSolvSphere,VSolvSphere_div,buff_shield,
-     & long_r_sidechain(ntyp),
+     & long_r_sidechain(ntyp),sslipi,sslipj,
+     & ssgradlipi,ssgradlipj,
      & short_r_sidechain(ntyp),fac_shield(maxres),wshield,
      & grad_shield_side(3,maxcontsshi,-1:maxres),
      &  grad_shield(3,-1:maxres),
-     &  grad_shield_loc(3,maxcontsshi,-1:maxres),
+     &  grad_shield_loc(3,maxcontsshi,-1:maxres),lipscale,
      & ishield_list(maxres),shield_list(maxcontsshi,maxres),
      & ees0plist(maxconts,maxres)
 

diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
index 0a638ef..44384a6 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F
@@ -6,7 +6,7 @@ C      if(r.lt.r_cut-rlamb) then
 C        sscale=1.0d0
 C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
 C        gamm=(r-(r_cut-rlamb))/rlamb
-        sscalelip=1.0d0+r*r*(2*r-3.0d0)
+        sscalelip=1.0d0+r*r*(2.0d0*r-3.0d0)
 C      else
 C        sscale=0d0
 C      endif
@@ -20,7 +20,7 @@ C     if(r.lt.r_cut-rlamb) then
 C        sscagrad=0.0d0
 C      else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
 C        gamm=(r-(r_cut-rlamb))/rlamb
-        sscagradlip=r*(6*r-6.0d0)
+        sscagradlip=r*(6.0d0*r-6.0d0)
 C      else
 C        sscagrad=0.0d0
 C      endif
@@ -2683,7 +2683,6 @@ C Uncomment following three lines for Ca-p interactions
        endif
 
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
-
           sss=sscale(1.0d0/(dsqrt(rrij)))
           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
           if (sss.lt.1.0d0) then
@@ -2704,7 +2703,7 @@ C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
              
             fac=-(evdwij+e1)*rrij*(1.0d0-sss)
-            fac=fac-(evdwij)*sssgrad*dsqrt(rrij)
+            fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/expon
             ggg(1)=xj*fac
             ggg(2)=yj*fac
             ggg(3)=zj*fac
@@ -2795,7 +2794,7 @@ C Uncomment following three lines for Ca-p interactions
           if (xj.lt.0) xj=xj+boxxsize
           yj=mod(yj,boxysize)
           if (yj.lt.0) yj=yj+boxysize
-          zj=mod(zi,boxzsize)
+          zj=mod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
       dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
       xj_safe=xj
@@ -2831,7 +2830,7 @@ C Uncomment following three lines for Ca-p interactions
           rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
           sss=sscale(1.0d0/(dsqrt(rrij)))
           sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
-          if (sss.gt.0.0d0) then
+            if (sss.gt.0.0d0) then
 
             fac=rrij**expon2
             e1=fac*fac*aad(itypj,iteli)
@@ -2845,11 +2844,12 @@ C Uncomment following three lines for Ca-p interactions
             evdw2=evdw2+evdwij*sss
             if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
      &          'evdw2',i,j,sss,evdwij
+
 C
 C Calculate contributions to the gradient in the virtual-bond and SC vectors.
 C
             fac=-(evdwij+e1)*rrij*sss
-            fac=fac+(evdwij)*sssgrad*dsqrt(rrij)
+            fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/expon
             ggg(1)=xj*fac
             ggg(2)=yj*fac
             ggg(3)=zj*fac

diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F
index 03bc765..c0670a7 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -3421,6 +3421,7 @@ C
       include 'COMMON.FFIELD'
       include 'COMMON.TIME1'
       include 'COMMON.SPLITELE'
+      include 'COMMON.SHIELD'
       dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
      &          erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
       double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
@@ -3537,6 +3538,30 @@ c        end if
           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
         call eelecij(i,i+2,ees,evdw1,eel_loc)
         if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
@@ -3589,13 +3614,35 @@ c        if ((zmedi.gt.((0.5d0)*boxzsize)).or.
 c     &       (zmedi.lt.((-0.5d0)*boxzsize))) then
 c        go to 196
 c        endif
-          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
-
+          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.0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0
+       endif
         num_conti=num_cont_hb(i)
 c        write(iout,*) "JESTEM W PETLI"
         call eelecij(i,i+3,ees,evdw1,eel_loc)
@@ -3637,6 +3684,29 @@ c     &  .or. itype(i-1).eq.ntyp1
           if (ymedi.lt.0) ymedi=ymedi+boxysize
           zmedi=mod(zmedi,boxzsize)
           if (zmedi.lt.0) zmedi=zmedi+boxzsize
+       if ((zmedi.gt.bordlipbot)
+     &.and.(zmedi.lt.bordliptop)) then
+C the energy transfer exist
+        if (zmedi.lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((zmedi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+        elseif (zmedi.gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-zmedi)/lipbufthick)
+         sslipi=sscalelip(fracinbuf)
+         ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+        else
+         sslipi=1.0d0
+         ssgradlipi=0.0
+        endif
+       else
+         sslipi=0.0d0
+         ssgradlipi=0.0
+       endif
+C         print *,sslipi,"TU?!"
 C          xmedi=xmedi+xshift*boxxsize
 C          ymedi=ymedi+yshift*boxysize
 C          zmedi=zmedi+zshift*boxzsize
@@ -3771,6 +3841,28 @@ C          zj=c(3,j)+0.5D0*dzj-zmedi
           zj=mod(zj,boxzsize)
           if (zj.lt.0) zj=zj+boxzsize
           if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
+       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
@@ -3865,13 +3957,20 @@ C          fac_shield(j)=0.6
           el2=el2*fac_shield(i)**2*fac_shield(j)**2
           eesij=(el1+el2)
           ees=ees+eesij
+C FOR NOW SHIELD IS NOT USED WITH LIPSCALE
+C     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           else
           fac_shield(i)=1.0
           fac_shield(j)=1.0
           eesij=(el1+el2)
           ees=ees+eesij
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+C          print *,"TUCC",(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
           endif
           evdw1=evdw1+evdwij*sss
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+C          print *,sslipi,sslipj,lipscale**2,
+C     &     (sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
 cd          write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
 cd     &      iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
 cd     &      1.0D0/dsqrt(rrmij),evdwij,eesij,
@@ -3891,7 +3990,9 @@ C Calculate contributions to the Cartesian gradient.
 C
 #ifdef SPLITELE
           facvdw=-6*rrmij*(ev1+evdwij)*sss
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           facel=-3*rrmij*(el1+eesij)
+     &*((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           fac1=fac
           erij(1)=xj*rmij
           erij(2)=yj*rmij
@@ -3990,6 +4091,12 @@ C            gelc_long(k,j-1)=gelc_long(k,j-1)
 C     &                    +grad_shield(k,j)*eesij/fac_shield(j)
           enddo
 C           print *,"bafter", gelc_long(1,i), gelc_long(1,j)
+C Lipidic part for lipscale
+            gelc_long(3,j)=gelc_long(3,j)+
+     &     ssgradlipj*eesij/2.0d0*lipscale**2
+
+            gelc_long(3,i)=gelc_long(3,i)+
+     &     ssgradlipi*eesij/2.0d0*lipscale**2
 
 *
 * Loop over residues i+1 thru j-1.
@@ -4001,8 +4108,13 @@ cgrad            enddo
 cgrad          enddo
           if (sss.gt.0.0) then
           ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           else
           ggg(1)=0.0
           ggg(2)=0.0
@@ -4018,6 +4130,12 @@ c 9/28/08 AL Gradient compotents will be summed only at the end
             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
           enddo
+C Lipidic part for scaling weight
+           gvdwpp(3,j)=gvdwpp(3,j)+
+     &     sss*ssgradlipj*evdwij/2.0d0*lipscale**2
+           gvdwpp(3,i)=gvdwpp(3,i)+
+     &     sss*ssgradlipi*evdwij/2.0d0*lipscale**2
+
 *
 * Loop over residues i+1 thru j-1.
 *
@@ -4029,6 +4147,7 @@ cgrad          enddo
 #else
 C MARYSIA
           facvdw=(ev1+evdwij)*sss
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           facel=(el1+eesij)
           fac1=fac
           fac=-3*rrmij*(facvdw+facvdw+facel)
@@ -4064,12 +4183,22 @@ cgrad            enddo
 cgrad          enddo
 c 9/28/08 AL Gradient compotents will be summed only at the end
           ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
           ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
+     &    *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           do k=1,3
             gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
             gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
           enddo
+           gvdwpp(3,j)=gvdwpp(3,j)+
+     &     sss*ssgradlipj*evdwij/2.0d0*lipscale**2
+           gvdwpp(3,i)=gvdwpp(3,i)+
+     &     sss*ssgradlipi*evdwij/2.0d0*lipscale**2
+
 #endif
 *
 * Angular part
@@ -4088,6 +4217,7 @@ cd   &          (dcosg(k),k=1,3)
           do k=1,3
             ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
      &      fac_shield(i)**2*fac_shield(j)**2
+     &      *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
           enddo
 c          do k=1,3
 c            ghalf=0.5D0*ggg(k)
@@ -4109,10 +4239,12 @@ C                     print *,"before22", gelc_long(1,i), gelc_long(1,j)
      &           +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
      &           + ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))
      &           *fac_shield(i)**2*fac_shield(j)**2   
+     &      *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
             gelc(k,j)=gelc(k,j)
      &           +((ecosa*(dc_norm(k,i)-cosa*dc_norm(k,j))
      &           + ecosg*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1))
      &           *fac_shield(i)**2*fac_shield(j)**2
+     &      *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
             gelc_long(k,j)=gelc_long(k,j)+ggg(k)
             gelc_long(k,i)=gelc_long(k,i)-ggg(k)
           enddo
@@ -4327,6 +4459,8 @@ C           fac_shield(j)=0.6
           endif
           eel_loc_ij=eel_loc_ij
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 C Now derivative over eel_loc
           if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
      &  (shield_mode.gt.0)) then
@@ -4383,6 +4517,8 @@ C Calculate patrial derivative for theta angle
      &     +a32*gmuij1(3)
      &     +a33*gmuij1(4))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 c         write(iout,*) "derivative over thatai"
 c         write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
 c     &   a33*gmuij1(4) 
@@ -4399,6 +4535,8 @@ c     &   a33*gmuij2(4)
          gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
      &      geel_loc_ij*wel_loc
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 
 c  Derivative over j residue
          geel_loc_ji=a22*gmuji1(1)
@@ -4412,6 +4550,7 @@ c     &   a33*gmuji1(4)
         gloc(nphi+j,icg)=gloc(nphi+j,icg)+
      &      geel_loc_ji*wel_loc
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
          geel_loc_ji=
      &     +a22*gmuji2(1)
@@ -4424,6 +4563,8 @@ c     &   a33*gmuji2(4)
          gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
      &      geel_loc_ji*wel_loc
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 #endif
 cd          write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
 
@@ -4440,22 +4581,35 @@ C Partial derivatives in virtual-bond dihedral angles gamma
      &            (a22*muder(1,i)*mu(1,j)+a23*muder(1,i)*mu(2,j)
      &           +a32*muder(2,i)*mu(1,j)+a33*muder(2,i)*mu(2,j))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
           gel_loc_loc(j-1)=gel_loc_loc(j-1)+ 
      &           (a22*mu(1,i)*muder(1,j)+a23*mu(1,i)*muder(2,j)
      &           +a32*mu(2,i)*muder(1,j)+a33*mu(2,i)*muder(2,j))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
           do l=1,3
             ggg(l)=(agg(l,1)*muij(1)+
      &          agg(l,2)*muij(2)+agg(l,3)*muij(3)+agg(l,4)*muij(4))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
             gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
             gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
 cgrad            ghalf=0.5d0*ggg(l)
 cgrad            gel_loc(l,i)=gel_loc(l,i)+ghalf
 cgrad            gel_loc(l,j)=gel_loc(l,j)+ghalf
           enddo
+            gel_loc_long(3,j)=gel_loc_long(3,j)+
+     &     ssgradlipj*eel_loc_ij/2.0d0*lipscale/
+     & ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+            gel_loc_long(3,i)=gel_loc_long(3,i)+
+     &     ssgradlipi*eel_loc_ij/2.0d0*lipscale/
+     & ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 cgrad          do k=i+1,j2
 cgrad            do l=1,3
 cgrad              gel_loc(l,k)=gel_loc(l,k)+ggg(l)
@@ -4466,18 +4620,22 @@ C Remaining derivatives of eello
             gel_loc(l,i)=gel_loc(l,i)+(aggi(l,1)*muij(1)+
      &        aggi(l,2)*muij(2)+aggi(l,3)*muij(3)+aggi(l,4)*muij(4))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
             gel_loc(l,i+1)=gel_loc(l,i+1)+(aggi1(l,1)*muij(1)+
      &     aggi1(l,2)*muij(2)+aggi1(l,3)*muij(3)+aggi1(l,4)*muij(4))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
             gel_loc(l,j)=gel_loc(l,j)+(aggj(l,1)*muij(1)+
      &       aggj(l,2)*muij(2)+aggj(l,3)*muij(3)+aggj(l,4)*muij(4))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
             gel_loc(l,j1)=gel_loc(l,j1)+(aggj1(l,1)*muij(1)+
      &     aggj1(l,2)*muij(2)+aggj1(l,3)*muij(3)+aggj1(l,4)*muij(4))
      &    *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
           enddo
           ENDIF
@@ -4723,7 +4881,42 @@ C Third- and fourth-order contributions from turns
      &    dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
      &    num_conti,j1,j2
       j=i+2
-c      write (iout,*) "eturn3",i,j,j1,j2
+C          xj=(c(1,j)+c(1,j+1))/2.0d0
+C          yj=(c(2,j)+c(2,j+1))/2.0d0
+          zj=(c(3,j)+c(3,j+1))/2.0d0
+C          xj=mod(xj,boxxsize)
+C          if (xj.lt.0) xj=xj+boxxsize
+C          yj=mod(yj,boxysize)
+C          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+          if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
+       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
+C      sslipj=0.0
+C      ssgradlipj=0.0d0
+      
+C      write (iout,*) "eturn3",i,j,j1,j2
       a_temp(1,1)=a22
       a_temp(1,2)=a23
       a_temp(2,1)=a32
@@ -4752,26 +4945,34 @@ c auxalary matrix for i+2 and constant i+1
         call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1))
         call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1))
         if (shield_mode.eq.0) then
-        fac_shield(i)=1.0
-        fac_shield(j)=1.0
+        fac_shield(i)=1.0d0
+        fac_shield(j)=1.0d0
 C        else
 C        fac_shield(i)=0.4
 C        fac_shield(j)=0.6
         endif
 C         if (j.eq.78)
 C     &   write(iout,*) i,j,fac_shield(i),fac_shield(j)
-        eello_turn3=eello_turn3+0.5d0*(pizda(1,1)+pizda(2,2))
+        eello_turn3=eello_turn3+
+C     &  1.0*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+     &0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
-        eello_t3=0.5d0*(pizda(1,1)+pizda(2,2))
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+        eello_t3=
+     &0.5d0*(pizda(1,1)+pizda(2,2))
      &  *fac_shield(i)*fac_shield(j)
 #ifdef NEWCORR
 C Derivatives in theta
         gloc(nphi+i,icg)=gloc(nphi+i,icg)
      &  +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)
      &  +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 #endif
 
 C        if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
@@ -4828,6 +5029,8 @@ C Derivatives in gamma(i)
         call matmat2(a_temp(1,1),auxmat3(1,1),pizda(1,1))
         gel_loc_turn3(i)=gel_loc_turn3(i)+0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 C Derivatives in gamma(i+1)
         call matmat2(EUg(1,1,i+1),EUgder(1,1,i+2),auxmat2(1,1))
         call transpose2(auxmat2(1,1),auxmat3(1,1))
@@ -4835,6 +5038,8 @@ C Derivatives in gamma(i+1)
         gel_loc_turn3(i+1)=gel_loc_turn3(i+1)
      &    +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 C Cartesian derivatives
         do l=1,3
 c            ghalf1=0.5d0*agg(l,1)
@@ -4849,6 +5054,7 @@ c            ghalf4=0.5d0*agg(l,4)
           gcorr3_turn(l,i)=gcorr3_turn(l,i)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
           a_temp(1,1)=aggi1(l,1)!+agg(l,1)
           a_temp(1,2)=aggi1(l,2)!+agg(l,2)
@@ -4858,6 +5064,7 @@ c            ghalf4=0.5d0*agg(l,4)
           gcorr3_turn(l,i+1)=gcorr3_turn(l,i+1)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
           a_temp(1,1)=aggj(l,1)!+ghalf1
           a_temp(1,2)=aggj(l,2)!+ghalf2
           a_temp(2,1)=aggj(l,3)!+ghalf3
@@ -4866,6 +5073,8 @@ c            ghalf4=0.5d0*agg(l,4)
           gcorr3_turn(l,j)=gcorr3_turn(l,j)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)
           a_temp(2,1)=aggj1(l,3)
@@ -4874,7 +5083,18 @@ c            ghalf4=0.5d0*agg(l,4)
           gcorr3_turn(l,j1)=gcorr3_turn(l,j1)
      &      +0.5d0*(pizda(1,1)+pizda(2,2))
      &   *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
         enddo
+         gshieldc_t3(3,i)=gshieldc_t3(3,i)+
+     &     ssgradlipi*eello_t3/4.0d0*lipscale
+         gshieldc_t3(3,j)=gshieldc_t3(3,j)+
+     &     ssgradlipj*eello_t3/4.0d0*lipscale
+         gshieldc_t3(3,i-1)=gshieldc_t3(3,i-1)+
+     &     ssgradlipi*eello_t3/4.0d0*lipscale
+         gshieldc_t3(3,j-1)=gshieldc_t3(3,j-1)+
+     &     ssgradlipj*eello_t3/4.0d0*lipscale
+
+C         print *,ssgradlipi,ssgradlipj,eello_t3,sslipi,sslipj
       return
       end
 C-------------------------------------------------------------------------------
@@ -4923,6 +5143,37 @@ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 cd        call checkint_turn4(i,a_temp,eello_turn4_num)
 c        write (iout,*) "eturn4 i",i," j",j," j1",j1," j2",j2
 c        write(iout,*)"WCHODZE W PROGRAM"
+          zj=(c(3,j)+c(3,j+1))/2.0d0
+C          xj=mod(xj,boxxsize)
+C          if (xj.lt.0) xj=xj+boxxsize
+C          yj=mod(yj,boxysize)
+C          if (yj.lt.0) yj=yj+boxysize
+          zj=mod(zj,boxzsize)
+          if (zj.lt.0) zj=zj+boxzsize
+C          if ((zj.lt.0).or.(xj.lt.0).or.(yj.lt.0)) write (*,*) "CHUJ"
+       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
+
         a_temp(1,1)=a22
         a_temp(1,2)=a23
         a_temp(2,1)=a32
@@ -5004,6 +5255,8 @@ C        fac_shield(j)=0.4
         endif
         eello_turn4=eello_turn4-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
         eello_t4=-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
 c             write(iout,*)'chujOWO', auxvec(1),b1(1,iti2)
@@ -5063,13 +5316,17 @@ cd     &    ' eello_turn4_num',8*eello_turn4_num
         gloc(nphi+i,icg)=gloc(nphi+i,icg)
      &                  -(gs13+gsE13+gsEE1)*wturn4
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
         gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)
      &                    -(gs23+gs21+gsEE2)*wturn4
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
         gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)
      &                    -(gs32+gsE31+gsEE3)*wturn4
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
 
 c         gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)-
 c     &   gs2
@@ -5087,6 +5344,8 @@ C Derivatives in gamma(i)
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i)=gel_loc_turn4(i)-(s1+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 C Derivatives in gamma(i+1)
         call transpose2(EUgder(1,1,i+2),e2tder(1,1))
         call matvec2(ae3(1,1),Ub2der(1,i+2),auxvec(1)) 
@@ -5096,6 +5355,8 @@ C Derivatives in gamma(i+1)
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i+1)=gel_loc_turn4(i+1)-(s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 C Derivatives in gamma(i+2)
         call transpose2(EUgder(1,1,i+3),e3tder(1,1))
         call matvec2(e1a(1,1),Ub2der(1,i+3),auxvec(1))
@@ -5108,6 +5369,8 @@ C Derivatives in gamma(i+2)
         s3=0.5d0*(pizda(1,1)+pizda(2,2))
         gel_loc_turn4(i+2)=gel_loc_turn4(i+2)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
 C Cartesian derivatives
 C Derivatives of this turn contributions in DC(i+2)
         if (j.lt.nres-1) then
@@ -5128,6 +5391,8 @@ C Derivatives of this turn contributions in DC(i+2)
             ggg(l)=-(s1+s2+s3)
             gcorr4_turn(l,i+2)=gcorr4_turn(l,i+2)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
           enddo
         endif
 C Remaining derivatives of this turn contribution
@@ -5147,6 +5412,8 @@ C Remaining derivatives of this turn contribution
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,i)=gcorr4_turn(l,i)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
           a_temp(1,1)=aggi1(l,1)
           a_temp(1,2)=aggi1(l,2)
           a_temp(2,1)=aggi1(l,3)
@@ -5162,6 +5429,8 @@ C Remaining derivatives of this turn contribution
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,i+1)=gcorr4_turn(l,i+1)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
           a_temp(1,1)=aggj(l,1)
           a_temp(1,2)=aggj(l,2)
           a_temp(2,1)=aggj(l,3)
@@ -5177,6 +5446,8 @@ C Remaining derivatives of this turn contribution
           s3=0.5d0*(pizda(1,1)+pizda(2,2))
           gcorr4_turn(l,j)=gcorr4_turn(l,j)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
           a_temp(1,1)=aggj1(l,1)
           a_temp(1,2)=aggj1(l,2)
           a_temp(2,1)=aggj1(l,3)
@@ -5193,7 +5464,16 @@ C Remaining derivatives of this turn contribution
 c          write (iout,*) "s1",s1," s2",s2," s3",s3," s1+s2+s3",s1+s2+s3
           gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
      &  *fac_shield(i)*fac_shield(j)
+     &*((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
         enddo
+         gshieldc_t4(3,i)=gshieldc_t4(3,i)+
+     &     ssgradlipi*eello_t4/4.0d0*lipscale
+         gshieldc_t4(3,j)=gshieldc_t4(3,j)+
+     &     ssgradlipj*eello_t4/4.0d0*lipscale
+         gshieldc_t4(3,i-1)=gshieldc_t4(3,i-1)+
+     &     ssgradlipi*eello_t4/4.0d0*lipscale
+         gshieldc_t4(3,j-1)=gshieldc_t4(3,j-1)+
+     &     ssgradlipj*eello_t4/4.0d0*lipscale
       return
       end
 C-----------------------------------------------------------------------------

diff --git a/source/unres/src_MD-M/energy_split-sep.F b/source/unres/src_MD-M/energy_split-sep.F
index b3f924b..f79deea 100644 (file)
--- a/source/unres/src_MD-M/energy_split-sep.F
+++ b/source/unres/src_MD-M/energy_split-sep.F
@@ -293,7 +293,7 @@ C
       call sum_energy(energia,.true.)
 C      call enerprint
 C      write (iout,*) "Exit ETOTAL_LONG"
-      call flush(iout)
+C      call flush(iout)
       return
       end
 c------------------------------------------------------------------------------
@@ -627,6 +627,6 @@ C      write (iout,*) "ETOTAL_SHORT before SUM_ENERGY"
       call flush(iout)
       call sum_energy(energia,.true.)
 C      write (iout,*) "Exit ETOTAL_SHORT"
-      call flush(iout)
+C      call flush(iout)
       return
       end

diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F
index e40717f..d0c10c7 100644 (file)
--- a/source/unres/src_MD-M/readrtns_CSA.F
+++ b/source/unres/src_MD-M/readrtns_CSA.F
@@ -766,10 +766,16 @@ C 12/1/95 Added weight for the multi-body term WCORR
       call reada(weightcard,"BTRISS",btriss,0.021D0)
       call reada(weightcard,"CTRISS",ctriss,1.001D0)
       call reada(weightcard,"DTRISS",dtriss,1.001D0)
+      call reada(weightcard,"LIPSCALE",lipscale,1.3D0)
       write (iout,*) "ATRISS=", atriss
       write (iout,*) "BTRISS=", btriss
       write (iout,*) "CTRISS=", ctriss
       write (iout,*) "DTRISS=", dtriss
+       if (shield_mode.gt.0) then
+        lipscale=0.0d0
+        write (iout,*) "liscale not used in shield mode"
+       endif
+      write (iout,*) "lipid scaling factor", lipscale
       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
       do i=1,maxres
         dyn_ss_mask(i)=.false.