Introduction of lipid energy transfer

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 7d2d27f..292bbee 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
@@ -1,4 +1,33 @@
 C-----------------------------------------------------------------------
+      double precision function sscalelip(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+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
+        sscale=1.0d0+r*r*(2*r-3.0d0)
+C      else
+C        sscale=0d0
+C      endif
+      return
+      end
+C-----------------------------------------------------------------------
+      double precision function sscagradlip(r)
+      double precision r,gamm
+      include "COMMON.SPLITELE"
+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
+        sscagrad=r*(6*r-6.0d0)
+C      else
+C        sscagrad=0.0d0
+C      endif
+      return
+      end
+
+C-----------------------------------------------------------------------
       double precision function sscale(r)
       double precision r,gamm
       include "COMMON.SPLITELE"

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 180f1bf..a8c0323 100644 (file)
--- a/source/unres/src_MD-M/energy_p_new_barrier.F
+++ b/source/unres/src_MD-M/energy_p_new_barrier.F
@@ -261,6 +261,12 @@ C  after the equilibration time
          Uconst=0.0d0
          Uconst_back=0.0d0
       endif
+C 01/27/2015 added by adasko
+C the energy component below is energy transfer into lipid environment 
+C based on partition function
+      if (wliptran.gt.0) then
+        call Eliptransfer
+      endif
 #ifdef TIMING
       time_enecalc=time_enecalc+MPI_Wtime()-time00
 #endif
@@ -302,6 +308,7 @@ C
       energia(17)=estr
       energia(20)=Uconst+Uconst_back
       energia(21)=esccor
+      energia(22)=eliptrans
 c    Here are the energies showed per procesor if the are more processors 
 c    per molecule then we sum it up in sum_energy subroutine 
 c      print *," Processor",myrank," calls SUM_ENERGY"
@@ -393,20 +400,21 @@ cMS$ATTRIBUTES C ::  proc_proc
       estr=energia(17)
       Uconst=energia(20)
       esccor=energia(21)
+      energia(22)=eliptrans
 #ifdef SPLITELE
       etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
      & +wang*ebe+wtor*etors+wscloc*escloc
      & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
-     & +wbond*estr+Uconst+wsccor*esccor
+     & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
 #else
       etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1)
      & +wang*ebe+wtor*etors+wscloc*escloc
      & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5
      & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3
      & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d
-     & +wbond*estr+Uconst+wsccor*esccor
+     & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran
 #endif
       energia(0)=etot
 c detecting NaNQ
@@ -658,6 +666,7 @@ c      enddo
      &                wturn6*gcorr6_turn(j,i)+
      &                wsccor*gsccorc(j,i)
      &               +wscloc*gscloc(j,i)
+     &               +wliptran*gliptranc(j,i)
 #else
           gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
      &                wel_loc*gel_loc(j,i)+
@@ -677,12 +686,14 @@ c      enddo
      &                wturn6*gcorr6_turn(j,i)+
      &                wsccor*gsccorc(j,i)
      &               +wscloc*gscloc(j,i)
+     &               +wliptran*gliptranc(j,i)
 #endif
           gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+
      &                  wbond*gradbx(j,i)+
      &                  wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+
      &                  wsccor*gsccorx(j,i)
      &                 +wscloc*gsclocx(j,i)
+     &                 +wliptran*gliptranx(j,i)
         enddo
       enddo 
 #ifdef DEBUG
@@ -9661,4 +9672,90 @@ crc      print *,((prod(i,j),i=1,2),j=1,2)
 
       return
       end
+CCC----------------------------------------------
+      subroutine Eliptransfer(eliptran)
+      include 'DIMENSIONS'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.LOCAL'
+      include 'COMMON.CHAIN'
+      include 'COMMON.DERIV'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.CALC'
+      include 'COMMON.CONTROL'
+      include 'COMMON.SPLITELE'
+      include 'COMMON.SBRIDGE'
+C structure of box:
+C      water
+C--bordliptop-- buffore starts
+C--bufliptop--- here true lipid starts
+C      lipid
+C--buflipbot--- lipid ends buffore starts
+C--bordlipbot--buffore ends
+      eliptran=0.0
+      do i=1,nres
+C first for peptide groups
+c for each residue check if it is in lipid or lipid water border area
+       if ((mod(c(3,i),boxzsize).gt.bordlipbot)
+     &.and.(mod(c(3,i),boxzsize).lt.bordliptop)) then
+C the energy transfer exist
+        if (mod(c(3,i),boxzsize).lt.buflipbot) then
+C what fraction I am in
+         fracinbuf=1.0d0-
+     &        ((mod(c(3,i),boxzsize)-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         ssslip=sscale(fracinbuf)
+         ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+         eliptran=eliptran+sslip
+         gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran
+C         print *,"doing sccale for lower part"
+        elseif (mod(c(3,i),boxzsize).gt.bufliptop) then
+         fracinbuf=1.0d0-((bordliptop-mod(c(3,i),boxzsize))/lipbufthick)
+         ssslip=sscale(fracinbuf)
+         ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+         eliptran=eliptran+sslip
+         gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran
+          print *, "doing sscalefor top part"
+        else
+         eliptran=eliptran+1.0d0
+         print *,"I am in true lipid"
+        endif
+C       else
+C       eliptran=elpitran+0.0 ! I am in water
+       endif
+       enddo
+C now multiply all by the peptide group transfer factor
+       eliptran=eliptran*pepliptran
+C now the same for side chains
+       do i=1,nres
+c for each residue check if it is in lipid or lipid water border area
+       if ((mod(c(3,i+nres),boxzsize).gt.bordlipbot)
+     & .and.(mod(c(3,i+nres),boxzsize).lt.bordliptop)) then
+C the energy transfer exist
+        if (mod(c(3,i+nres),boxzsize).lt.buflipbot) then
+         fracinbuf=1.0d0-
+     &     ((mod(c(3,i+nres),boxzsize)-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+         ssslip=sscale(fracinbuf)
+         ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+         eliptran=eliptran+sslip*liptranene(itype(i))
+         gliptranx(3,i)=gliptranx(3,i)+ssgradlip*liptranene(itype(i))
+         print *,"doing sccale for lower part"
+        elseif (mod(c(3,i+nres),boxzsize).gt.bufliptop) then
+         fracinbuf=1.0d0-
+     &((bordliptop-mod(c(3,i+nres),boxzsize))/lipbufthick)
+         ssslip=sscale(fracinbuf)
+         ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+         eliptran=eliptran+sslip*liptranene(itype(i))
+         gliptranx(3,i)=gliptranx(3,i)+ssgradlip*liptranene(itype(i))
+          print *, "doing sscalefor top part"
+        else
+         eliptran=eliptran+liptranene(itype(i))
+         print *,"I am in true lipid"
+        endif
+C       else
+C       eliptran=elpitran+0.0 ! I am in water
+       enddo