C
C Calculate the SC local energy.
C
+C print *,"TU DOCHODZE?"
call esc(escloc)
c print *,"Processor",myrank," computed USC"
C
else
esccor=0.0d0
endif
+C print *,"PRZED MULIt"
c print *,"Processor",myrank," computed Usccorr"
C
C 12/1/95 Multi-body terms
C 01/27/2015 added by adasko
C the energy component below is energy transfer into lipid environment
C based on partition function
+C print *,"przed lipidami"
if (wliptran.gt.0) then
- call Eliptransfer
+ call Eliptransfer(eliptran)
endif
+C print *,"za lipidami"
#ifdef TIMING
time_enecalc=time_enecalc+MPI_Wtime()-time00
#endif
energia(17)=estr
energia(20)=Uconst+Uconst_back
energia(21)=esccor
- energia(22)=eliptrans
+ energia(22)=eliptran
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"
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
- energia(22)=eliptrans
+ eliptran=energia(22)
#ifdef SPLITELE
etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1
& +wang*ebe+wtor*etors+wscloc*escloc
#ifdef MPI
include 'mpif.h'
#endif
- double precision gradbufc(3,maxres),gradbufx(3,maxres),
- & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,maxres)
+ double precision gradbufc(3,0:maxres),gradbufx(3,0:maxres),
+ & glocbuf(4*maxres),gradbufc_sum(3,0:maxres),gloc_scbuf(3,0:maxres)
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
& wstrain*ghpbc(j,i)
+ & +wliptran*gliptranc(j,i)
+
enddo
enddo
#else
& wcorr6*gradcorr6_long(j,i)+
& wturn6*gcorr6_turn_long(j,i)+
& wstrain*ghpbc(j,i)
+ & +wliptran*gliptranc(j,i)
enddo
enddo
#endif
estr=energia(17)
Uconst=energia(20)
esccor=energia(21)
+ eliptran=energia(22)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp,
& estr,wbond,ebe,wang,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor,
& edihcnstr,ebr*nss,
- & Uconst,etot
+ & Uconst,eliptran,wliptran,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST= ',1pE16.6,' (Constraint energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,
& ecorr,wcorr,
& ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3,
& eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr,
- & ebr*nss,Uconst,etot
+ & ebr*nss,Uconst,eliptran,wliptran,etot
10 format (/'Virtual-chain energies:'//
& 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/
& 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/
& 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
& 'UCONST=',1pE16.6,' (Constraint energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#endif
return
c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
+C have you changed here?
e1=fac*fac*aa(itypi,itypj)
e2=fac*bb(itypi,itypj)
evdwij=e1+e2
rij=1.0D0/r_inv_ij
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
+C have you changed here?
e1=fac*fac*aa(itypi,itypj)
e2=fac*bb(itypi,itypj)
evdwij=e_augm+e1+e2
call sc_angular
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
+C have you changed here?
fac=(rrij*sigsq)**expon2
e1=fac*fac*aa(itypi,itypj)
e2=fac*bb(itypi,itypj)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
+C here to start with
+C if (c(i,3).gt.
e1=fac*fac*aa(itypi,itypj)
e2=fac*bb(itypi,itypj)
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
gg(3)=zj*fac
C Calculate angular part of the gradient.
call sc_grad
+ endif
ENDIF ! dyn_ss
enddo ! j
enddo ! iint
C
C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition
do i=iturn3_start,iturn3_end
+ if (i.le.1) cycle
+C write(iout,*) "tu jest i",i
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+2).eq.ntyp1
& .or. itype(i+3).eq.ntyp1
num_cont_hb(i)=num_conti
enddo
do i=iturn4_start,iturn4_end
+ if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+3).eq.ntyp1
& .or. itype(i+4).eq.ntyp1
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
do i=iatel_s,iatel_e
+ if (i.le.1) cycle
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+2).eq.ntyp1
& .or. itype(i-1).eq.ntyp1
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
-c write (iout,*) i,j,itype(i),itype(j)
+C write (iout,*) i,j
+ if (j.le.1) cycle
if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1
& .or.itype(j+2).eq.ntyp1
& .or.itype(j-1).eq.ntyp1
c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
- if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
- & iabs(itype(jjj)).eq.1) then
+C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+C & iabs(itype(jjj)).eq.1) then
cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
if (.not.dyn_ss .and. i.le.nss) then
C 15/02/13 CC dynamic SSbond - additional check
if (ii.gt.nres
& .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
->>>>>>> f5379d3246c4bd95e946c4d35d4a1c13e329c4cb
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
endif
end
CCC----------------------------------------------
subroutine Eliptransfer(eliptran)
+ implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
include 'COMMON.SBRIDGE'
+C print *,"wchodze"
C structure of box:
C water
C--bordliptop-- buffore starts
C--buflipbot--- lipid ends buffore starts
C--bordlipbot--buffore ends
eliptran=0.0
- do i=1,nres
+ do i=ilip_start,ilip_end
+ if (itype(i).eq.ntyp1) cycle
+
+ positi=(mod((c(3,i)+c(3,i+1)),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+C print *,i
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
+ if ((positi.gt.bordlipbot)
+ &.and.(positi.lt.bordliptop)) then
C the energy transfer exist
- if (mod(c(3,i),boxzsize).lt.buflipbot) then
+ if (positi.lt.buflipbot) then
C what fraction I am in
fracinbuf=1.0d0-
- & ((mod(c(3,i),boxzsize)-bordlipbot)/lipbufthick)
+ & ((positi-bordlipbot)/lipbufthick)
C lipbufthick is thickenes of lipid buffore
- ssslip=sscale(fracinbuf)
+ sslip=sscalelip(fracinbuf)
ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
- eliptran=eliptran+sslip
- gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0
+C gliptranc(3,i-2)=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)
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
ssgradlip=sscagradlip(fracinbuf)/lipbufthick
- eliptran=eliptran+sslip
- gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran
+ eliptran=eliptran+sslip*pepliptran
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
+C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
print *, "doing sscalefor top part"
else
- eliptran=eliptran+1.0d0
+ eliptran=eliptran+pepliptran
print *,"I am in true lipid"
endif
C else
C eliptran=elpitran+0.0 ! I am in water
endif
enddo
+C print *, "nic nie bylo w lipidzie?"
C now multiply all by the peptide group transfer factor
- eliptran=eliptran*pepliptran
+C eliptran=eliptran*pepliptran
C now the same for side chains
- do i=1,nres
+ do i=ilip_start,ilip_end
+ if (itype(i).eq.ntyp1) cycle
+ positi=(mod(c(3,i+nres),boxzsize))
+ if (positi.le.0) positi=positi+boxzsize
+C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
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 respos=mod(c(3,i+nres),boxzsize)
+ if ((positi.gt.bordlipbot)
+ & .and.(positi.lt.bordliptop)) then
C the energy transfer exist
- if (mod(c(3,i+nres),boxzsize).lt.buflipbot) then
+ if (positi.lt.buflipbot) then
fracinbuf=1.0d0-
- & ((mod(c(3,i+nres),boxzsize)-bordlipbot)/lipbufthick)
+ & ((positi-bordlipbot)/lipbufthick)
C lipbufthick is thickenes of lipid buffore
- ssslip=sscale(fracinbuf)
+ sslip=sscalelip(fracinbuf)
ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
eliptran=eliptran+sslip*liptranene(itype(i))
- gliptranx(3,i)=gliptranx(3,i)+ssgradlip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))/2.0d0
+ gliptranc(3,i-1)=
+ &+ssgradlip*liptranene(itype(i))
print *,"doing sccale for lower part"
- elseif (mod(c(3,i+nres),boxzsize).gt.bufliptop) then
+ elseif (positi.gt.bufliptop) then
fracinbuf=1.0d0-
- &((bordliptop-mod(c(3,i+nres),boxzsize))/lipbufthick)
- ssslip=sscale(fracinbuf)
+ &((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(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"
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))/2.0d0
+ gliptranc(3,i-1)=
+ &+ssgradlip*liptranene(itype(i))
+ print *, "doing sscalefor top part",sslip,fracinbuf
else
eliptran=eliptran+liptranene(itype(i))
print *,"I am in true lipid"
endif
+ endif ! if in lipid or buffor
C else
C eliptran=elpitran+0.0 ! I am in water
enddo
-
+ return
+ end
projects / unres.git / blobdiff