From ffd0e6dcc64b1ef30773e02b4ef72c1b2455fea5 Mon Sep 17 00:00:00 2001 From: Adam Sieradzan Date: Thu, 29 Jan 2015 08:46:52 +0100 Subject: [PATCH] Introduction of lipid energy transfer --- source/unres/src_MD-M/energy_p_new-sep_barrier.F | 29 +++++++ source/unres/src_MD-M/energy_p_new_barrier.F | 101 +++++++++++++++++++++- 2 files changed, 128 insertions(+), 2 deletions(-) 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 --- 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 --- 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 -- 1.7.9.5