C 21/5/07 Calculate local sicdechain correlation energy
C
call eback_sc_corr(esccor)
+
+ if (wliptran.gt.0) then
+ call Eliptransfer(eliptran)
+ endif
+
C
C 12/1/95 Multi-body terms
C
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+ & +wliptran*eliptran
#else
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2
& +welec*fact(1)*(ees+evdw1)
& +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6
& +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d
& +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr
+ & +wliptran*eliptran
#endif
energia(0)=etot
energia(1)=evdw
energia(20)=edihcnstr
energia(21)=evdw_t
energia(24)=ethetacnstr
+ energia(22)=eliptran
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
& wcorr6*fact(5)*gradcorr6(j,i)+
& wturn6*fact(5)*gcorr6_turn(j,i)+
& wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(j,i)
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*fact(2)*gsccorx(j,i)
+ & +wliptran*gliptranx(j,i)
enddo
#else
do i=1,nct
& wcorr6*fact(5)*gradcorr6(j,i)+
& wturn6*fact(5)*gcorr6_turn(j,i)+
& wsccor*fact(2)*gsccorc(j,i)
+ & +wliptran*gliptranc(j,i)
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*fact(1)*gsccorx(j,i)
+ & +wliptran*gliptranx(j,i)
enddo
#endif
enddo
edihcnstr=energia(20)
estr=energia(18)
ethetacnstr=energia(24)
+ eliptran=energia(22)
#ifdef SPLITELE
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1,
& wvdwpp,
& ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5),
& eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2),
& eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5),
- & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot
+ & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,
+ & 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)'/
& 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/
+ & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/
& 'ETOT= ',1pE16.6,' (total)')
#else
write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),estr,wbond,
& ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2),
& eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3),
& eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor,
- & edihcnstr,ethetacnstr,ebr*nss,etot
+ & edihcnstr,ethetacnstr,ebr*nss,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)'/
& 'ETHETC= ',1pE16.6,' (valence angle constraints)'/
& 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic 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
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
ij=icant(itypi,itypj)
c ROZNICA z cluster
cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
rij=1.0D0/r_inv_ij
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
ij=icant(itypi,itypj)
eneps_temp(1,ij)=eneps_temp(1,ij)+(e1+a_augm)
cd & bb(itypi,itypj),augm(itypi,itypj),epsi,sigm,
cd & sigma(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
eneps_temp(1,ij)=eneps_temp(1,ij)+e1*aux
& /dabs(eps(itypi,itypj))
eneps_temp(2,ij)=eneps_temp(2,ij)+e2*aux/eps(itypi,itypj)
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij
else
evdw_t=evdw_t+evdwij
endif
if (calc_grad) then
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),15(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,chi1,chi2,chip1,chip2,
if (yi.lt.0) yi=yi+boxysize
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
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
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+C write(iout,*),aa,aa_lip(itypi,itypj),aa_aq(itypi,itypj)
C checking the distance
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
if (sss.le.0.0) cycle
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
+
call sc_angular
sigsq=1.0D0/sigsq
sig=sig0ij*dsqrt(sigsq)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
- if (bb(itypi,itypj).gt.0) then
+ if (bb.gt.0) then
evdw=evdw+evdwij*sss
else
evdw_t=evdw_t+evdwij*sss
c & " ij",ij," eneps",aux*e1/dabs(eps(itypi,itypj)),
c & aux*e2/eps(itypi,itypj)
c if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
#ifdef DEBUG
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
evdwij=evdwij*eps2rt*eps3rt
- if (bb(itypi,itypj).gt.0.0d0) then
+ if (bb.gt.0.0d0) then
evdw=evdw+evdwij+e_augm
else
evdw_t=evdw_t+evdwij+e_augm
do k=1,2
mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1)
enddo
+C write (iout,*) 'mumu',i,b1(1,iti),Ub2(1,i-2)
+
C Vectors and matrices dependent on a single virtual-bond dihedral.
call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1))
call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2))
eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
& +a33*muij(4)
c write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
-c write (iout,'(a6,2i5,0pf7.3)')
-c & 'eelloc',i,j,eel_loc_ij
+C write (iout,'(a6,2i5,0pf7.3)')
+C & 'eelloc',i,j,eel_loc_ij
+C write(iout,*) 'muije=',i,j,muij(1),muij(2),muij(3),muij(4)
c write (iout,*) a22,muij(1),a23,muij(2),a32,muij(3)
eel_loc=eel_loc+eel_loc_ij
C Partial derivatives in virtual-bond dihedral angles gamma
enddo
endif
else if (j.eq.i+3 .and. itype(i+2).ne.ntyp1) then
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+C changes suggested by Ana to avoid out of bounds
+ & .or.((i+5).gt.nres)
+ & .or.((i-1).le.0)
+C end of changes suggested by Ana
+ & .or. itype(i+3).eq.ntyp1
+ & .or. itype(i+4).eq.ntyp1
+ & .or. itype(i+5).eq.ntyp1
+ & .or. itype(i).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1) goto 178
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Fourth-order contributions
gcorr4_turn(l,j1)=gcorr4_turn(l,j1)-(s1+s2+s3)
enddo
endif
+ 178 continue
endif
return
end
return
end
crc-------------------------------------------------
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+ subroutine Eliptransfer(eliptran)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.GEO'
+ include 'COMMON.VAR'
+ include 'COMMON.LOCAL'
+ include 'COMMON.CHAIN'
+ include 'COMMON.DERIV'
+ include 'COMMON.INTERACT'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CALC'
+ include 'COMMON.CONTROL'
+ include 'COMMON.SPLITELE'
+ include 'COMMON.SBRIDGE'
+C this is done by Adasko
+C print *,"wchodze"
+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 do i=1,1
+ if (itype(i).eq.ntyp1) cycle
+
+ positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),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 ((positi.gt.bordlipbot)
+ &.and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ 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
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=sscagradlip(fracinbuf)/lipbufthick
+ 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
+C print *, "doing sscalefor top part"
+C print *,i,sslip,fracinbuf,ssgradlip
+ else
+ eliptran=eliptran+pepliptran
+C 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
+C eliptran=eliptran*pepliptran
+C now the same for side chains
+CV do i=1,1
+ do i=1,nres
+ 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
+C respos=mod(c(3,i+nres),boxzsize)
+C print *,positi,bordlipbot,buflipbot
+ if ((positi.gt.bordlipbot)
+ & .and.(positi.lt.bordliptop)) then
+C the energy transfer exist
+ if (positi.lt.buflipbot) then
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslip=sscalelip(fracinbuf)
+ ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
+ eliptran=eliptran+sslip*liptranene(itype(i))
+ gliptranx(3,i)=gliptranx(3,i)
+ &+ssgradlip*liptranene(itype(i))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *,"doing sccale for lower part"
+ elseif (positi.gt.bufliptop) then
+ fracinbuf=1.0d0-
+ &((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))
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))
+C print *, "doing sscalefor top part",sslip,fracinbuf
+ else
+ eliptran=eliptran+liptranene(itype(i))
+C 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
+
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
SUBROUTINE MATVEC2(A1,V1,V2)
implicit real*8 (a-h,o-z)
include 'DIMENSIONS'
return
end
C-----------------------------------------------------------------------
+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
+ sscalelip=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
+ sscagradlip=r*(6*r-6.0d0)
+C else
+C sscagrad=0.0d0
+C endif
+ return
+ end
+
+C-----------------------------------------------------------------------
projects / unres.git / blobdiff