integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
& sigij,r0ij,rcut,sss1,sssgrad1,sqrij
- double precision sscale,sscagrad
+ double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
+ & faclip
+ double precision sscale,sscagrad,sscagradlip,sscalelip
double precision boxshift
+ double precision gg_lipi(3),gg_lipj(3)
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
c do i=iatsc_s,iatsc_e
do ikont=g_listscsc_start,g_listscsc_end
i=newcontlisti(ikont)
yi=c(2,nres+i)
zi=c(3,nres+i)
call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
C
C Calculate SC interaction energy.
C
yj=c(2,nres+j)
zj=c(3,nres+j)
call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ 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
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
if (sss.lt.1.0d0) then
rrij=1.0D0/rij
fac=rrij**expon2
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=e1+e2
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=(sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
- gvdwc(k,i)=gvdwc(k,i)-gg(k)
- gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
enddo
endif
c enddo ! j
integer i,j,k,itypi,itypj,itypi1,num_conti,iint,ikont
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
& sigij,r0ij,rcut,sqrij,sss1,sssgrad1
- double precision sscale,sscagrad
+ double precision sscale,sscagrad,sscagradlip,sscalelip
+ double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
+ & faclip
double precision boxshift
+ double precision gg_lipi(3),gg_lipj(3)
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
c do i=iatsc_s,iatsc_e
do ikont=g_listscsc_start,g_listscsc_end
i=newcontlisti(ikont)
yi=c(2,nres+i)
zi=c(3,nres+i)
call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
C Change 12/1/95
num_conti=0
C
yj=c(2,nres+j)
zj=c(3,nres+j)
call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ 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
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=e1+e2
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=(sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
- gvdwc(k,i)=gvdwc(k,i)-gg(k)
- gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
enddo
endif
c enddo ! j
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
& fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
logical scheck
- double precision sscale,sscagrad
+ double precision sscale,sscagrad,sscagradlip,sscalelip
+ double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
+ & faclip
double precision boxshift
+ double precision gg_lipi(3),gg_lipj(3)
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
c do i=iatsc_s,iatsc_e
do ikont=g_listscsc_start,g_listscsc_end
i=newcontlisti(ikont)
yi=c(2,nres+i)
zi=c(3,nres+i)
call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
C
C Calculate SC interaction energy.
C
yj=c(2,nres+j)
zj=c(3,nres+j)
call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ 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
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
& sscagrad(rij/sigma(itypi,itypj),r_cut_respa)
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=e_augm+e1+e2
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=(sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
- gvdwc(k,i)=gvdwc(k,i)-gg(k)
- gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
enddo
endif
c enddo ! j
double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
& fac_augm,e_augm,r_inv_ij,r_shift_inv,sss1,sssgrad1
logical scheck
- double precision sscale,sscagrad
+ double precision sscale,sscagrad,sscagradlip,sscalelip
+ double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
+ & faclip
double precision boxshift
+ double precision gg_lipi(3),gg_lipj(3)
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
c do i=iatsc_s,iatsc_e
do ikont=g_listscsc_start,g_listscsc_end
i=newcontlisti(ikont)
yi=c(2,nres+i)
zi=c(3,nres+i)
call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
C
C Calculate SC interaction energy.
C
yj=c(2,nres+j)
zj=c(3,nres+j)
call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ 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
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
if (sss.gt.0.0d0) then
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=e_augm+e1+e2
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=(sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
- gvdwc(k,i)=gvdwc(k,i)-gg(k)
- gvdwc(k,j)=gvdwc(k,j)+gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
enddo
endif
c enddo ! j
integer itypi,itypj,itypi1,iint,ind,ikont
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
double precision sss1,sssgrad1
- double precision sscale,sscagrad
+ double precision sscale,sscagrad,sscagradlip,sscalelip
+ double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
+ & faclip
double precision boxshift
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
- evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
c if (icall.eq.0) then
c lprn=.true.
c else
yi=c(2,nres+i)
zi=c(3,nres+i)
call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
yj=c(2,nres+j)
zj=c(3,nres+j)
call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ 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
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
C Calculate the angular part of the gradient and sum add the contributions
C to the appropriate components of the Cartesian gradient.
call sc_grad_scale((1.0d0-sss)*sss1)
double precision evdw
integer itypi,itypj,itypi1,iint,ind,ikont
double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
- double precision sscale,sscagrad
+ double precision sscale,sscagrad,sscagradlip,sscalelip
+ double precision fracinbuf,sslipi,sslipj,ssgradlipj,ssgradlipi,
+ & faclip
double precision boxshift
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
c print *,'Entering EBP nnt=',nnt,' nct=',nct,' expon=',expon
c if (icall.eq.0) then
c lprn=.true.
yi=c(2,nres+i)
zi=c(3,nres+i)
call to_box(xi,yi,zi)
+ call lipid_layer(xi,yi,zi,sslipi,ssgradlipi)
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
yj=c(2,nres+j)
zj=c(3,nres+j)
call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ 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
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=(sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj))))/expon
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
+ do k=1,3
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
+ gvdwc(k,i)=gvdwc(k,i)-gg(k)+gg_lipi(k)
+ gvdwc(k,j)=gvdwc(k,j)+gg(k)+gg_lipj(k)
+ enddo
C Calculate the angular part of the gradient and sum add the contributions
C to the appropriate components of the Cartesian gradient.
call sc_grad_scale(sss)
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
- evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
lprn=.false.
c if (icall.eq.0) lprn=.false.
ind=0
call to_box(xj,yj,zj)
call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +aa_aq(itypi,itypj)*(2.0d0-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
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
& evdwij
endif
- if (energy_dec) write (iout,'(a6,2i5,4f10.5)')
- & 'evdw',i,j,rij,sss,sss1,evdwij
+ if (energy_dec) write (iout,'(a,2i5,5f10.5,e15.5)')
+ & 'r sss evdw',i,j,1.0d0/rij,sss1,sss,sslipi,sslipj,evdwij
C Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
- gg_lipi(3)=ssgradlipi*evdwij
- gg_lipj(3)=ssgradlipj*evdwij
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
+
C Calculate angular part of the gradient.
call sc_grad_scale((1.0d0-sss)*sss1)
endif
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
- evdw=0.0D0
+ gg_lipi=0.0D0
+ gg_lipj=0.0d0
lprn=.false.
c if (icall.eq.0) lprn=.false.
ind=0
call to_box(xj,yj,zj)
call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
- & +aa_aq(itypi,itypj)*(2.0d0-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
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+c write (iout,*) "aa bb",aa_lip(itypi,itypj),
+c & bb_lip(itypi,itypj),aa_aq(itypi,itypj),
+c & bb_aq(itypi,itypj),aa,bb
+c write (iout,*) (sslipi+sslipj)/2.0d0,
+c & (2.0d0-sslipi-sslipj)/2.0d0
xj=boxshift(xj-xi,boxxsize)
yj=boxshift(yj-yi,boxysize)
zj=boxshift(zj-zi,boxzsize)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
& evdwij
endif
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'evdw',i,j,evdwij
+ if (energy_dec) write (iout,'(a,2i5,4f10.5,e15.5)')
+ & 'r sss evdw',i,j,1.0d0/rij,sss,sslipi,sslipj,evdwij
C Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
- gg_lipi(3)=ssgradlipi*evdwij
- gg_lipj(3)=ssgradlipj*evdwij
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
+c write (iout,*) "gglip",i,j,gg_lipi,gg_lipj
C Calculate angular part of the gradient.
call sc_grad_scale(sss)
endif
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
double precision sss1,sssgrad1
evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
lprn=.false.
c if (icall.eq.0) lprn=.true.
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ & *(eps3rt*eps3rt)*sss1*(1.0d0-sss)/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
C Calculate angular part of the gradient.
call sc_grad_scale((1.0d0-sss)*sss1)
endif
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
double precision boxshift
evdw=0.0D0
+ gg_lipi=0.0d0
+ gg_lipj=0.0d0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
lprn=.false.
c if (icall.eq.0) lprn=.true.
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
+ faclip=fac
e1=fac*fac*aa
e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ & *(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ & +faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
C Calculate angular part of the gradient.
call sc_grad_scale(sss)
endif
C
C Calculate the components of the gradient in DC and X
C
+c write (iout,*) "scgrad gglip",i,j,gg_lipi,gg_lipj
do l=1,3
gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
common /locel/ a_temp,agg,aggi,aggi1,aggj,aggj1,a22,a23,a32,a33,
& dxi,dyi,dzi,dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
& num_conti,j1,j2
+ double precision sslipi,sslipj,ssgradlipi,ssgradlipj
+ common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
double precision scal_el /1.0d0/
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
num_conti=0
call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
#ifdef FOURBODY
num_conti=num_cont_hb(i)
#endif
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
#ifdef FOURBODY
num_conti=num_cont_hb(i)
double precision sscale,sscagrad
double precision scalar
double precision boxshift
-
+ double precision sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij,
+ & faclipij2
+ common /lipcalc/ sslipi,sslipj,ssgradlipi,ssgradlipj,faclipij
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
double precision scal_el /1.0d0/
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ faclipij=(sslipi+sslipj)/2.0d0*lipscale+1.0d0
+ faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
xj=boxshift(xj-xmedi,boxxsize)
yj=boxshift(yj-ymedi,boxysize)
zj=boxshift(zj-zmedi,boxzsize)
eesij=el1+el2
C 12/26/95 - for the evaluation of multi-body H-bonding interactions
ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
- ees=ees+eesij*sss1
- evdw1=evdw1+evdwij*(1.0d0-sss)*sss1
+ ees=ees+eesij*sss1*faclipij2
+ evdw1=evdw1+evdwij*(1.0d0-sss)*sss1*faclipij2
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,
cd & xmedi,ymedi,zmedi,xj,yj,zj
if (energy_dec) then
- write (iout,'(a6,2i5,0pf7.3,2f7.3)')
- & 'evdw1',i,j,evdwij,sss,sss1
- write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
+ write (iout,'(a6,2i5,0pf7.3,2i5,e11.3,5f10.5)')
+ & 'evdw1',i,j,evdwij,iteli,itelj,aaa,sss,sss1,sssgrad,sssgrad1,rij
+ write (iout,'(a6,2i5,0pf7.3,6f8.5)') 'ees',i,j,eesij,
+ & fac_shield(i),fac_shield(j),sslipi,sslipj,faclipij,faclipij2
endif
C
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
- aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij
+c old aux=(facel+sssgrad1*(1.0d0-sss)*eesij*rmij)*faclipij2
+ aux=(facel+sssgrad1*eesij*rmij)*faclipij2
c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
ggg(1)=aux*xj
ggg(2)=aux*yj
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+ gelc_long(3,j)=gelc_long(3,j)+
+ & ssgradlipj*eesij/2.0d0*lipscale**2*sss1
+ gelc_long(3,i)=gelc_long(3,i)+
+ & ssgradlipi*eesij/2.0d0*lipscale**2*sss1
c gelc_long(3,i)=gelc_long(3,i)+
c ssgradlipi*eesij/2.0d0*lipscale**2*sss1
*
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
- facvdw=facvdw+
- & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij
-c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ facvdw=(facvdw+
+ &(-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij)
+ & *faclipij2
ggg(1)=facvdw*xj
ggg(2)=facvdw*yj
ggg(3)=facvdw*zj
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)+
+ & sss1*(1.0d0-sss)*ssgradlipj*evdwij/2.0d0*lipscale**2
+ gvdwpp(3,i)=gvdwpp(3,i)+
+ & sss1*(1.0d0-sss)*ssgradlipi*evdwij/2.0d0*lipscale**2
*
* Loop over residues i+1 thru j-1.
*
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
- aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
- & *eesij*rmij
+ aux=(fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
+ & *eesij*rmij)*faclipij2
c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
ggg(1)=aux*xj
ggg(2)=aux*yj
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
enddo
+ gvdwpp(3,j)=gvdwpp(3,j)+
+ & sss1*ssgradlipj*evdwij/2.0d0*lipscale**2
+ gvdwpp(3,i)=gvdwpp(3,i)+
+ & sss1*ssgradlipi*evdwij/2.0d0*lipscale**2
#endif
*
* Angular part
cd & (dcosg(k),k=1,3)
do k=1,3
ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1
- & *fac_shield(i)**2*fac_shield(j)**2
+ & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
enddo
gelc(k,i)=gelc(k,i)
& +((ecosa*(dc_norm(k,j)-cosa*dc_norm(k,i))
& +ecosb*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1))*sss1
- & *fac_shield(i)**2*fac_shield(j)**2
+ & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
c & *((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))*sss1
- & *fac_shield(i)**2*fac_shield(j)**2
+ & *fac_shield(i)**2*fac_shield(j)**2*faclipij2
c & *((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)
C fac_shield(j)=0.6
endif
eel_loc_ij=eel_loc_ij
- & *fac_shield(i)*fac_shield(j)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
eel_loc=eel_loc+eel_loc_ij
if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
& +a23*gmuij1(2)
& +a32*gmuij1(3)
& +a33*gmuij1(4))
- & *fac_shield(i)*fac_shield(j)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c write(iout,*) "derivative over thatai"
c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
c & a33*gmuij1(4)
& +a33*gmuij2(4)
gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
& geel_loc_ij*wel_loc
- & *fac_shield(i)*fac_shield(j)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c Derivative over j residue
geel_loc_ji=a22*gmuji1(1)
gloc(nphi+j,icg)=gloc(nphi+j,icg)+
& geel_loc_ji*wel_loc
- & *fac_shield(i)*fac_shield(j)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
geel_loc_ji=
& +a22*gmuji2(1)
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)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
#endif
cC Paral derivatives in virtual-bond dihedral angles gamma
if (i.gt.1)
& gel_loc_loc(i-1)=gel_loc_loc(i-1)+
& (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)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c & *fac_shield(i)*fac_shield(j)
c & *((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)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c & *fac_shield(i)*fac_shield(j)
c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
do l=1,3
ggg(l)=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)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c & *fac_shield(i)*fac_shield(j)
c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
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/faclipij
+ gel_loc_long(3,i)=gel_loc_long(3,i)+
+ & ssgradlipi*eel_loc_ij/2.0d0*lipscale/faclipij
cgrad do k=i+1,j2
cgrad do l=1,3
cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
do l=1,3
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)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c & *((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)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c & *((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)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c & *((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)*sss1
+ & *fac_shield(i)*fac_shield(j)*sss1*faclipij
c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
enddo
double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
& dist_temp, dist_init,sss_grad
double precision sscale,sscagrad
+ double precision sslipi,ssgradlipi,sslipj,ssgradlipj
double precision boxshift
integer ikont
+ double precision faclipij2
evdw1=0.0D0
C print *,"WCHODZE"
c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
call to_box(xmedi,ymedi,zmedi)
+ call lipid_layer(xmedi,ymedi,zmedi,sslipi,ssgradlipi)
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
c & ' ielend',ielend_vdw(i)
yj=c(2,j)+0.5D0*dyj
zj=c(3,j)+0.5D0*dzj
call to_box(xj,yj,zj)
+ call lipid_layer(xj,yj,zj,sslipj,ssgradlipj)
+ faclipij2=(sslipi+sslipj)/2.0d0*lipscale**2+1.0d0
xj=boxshift(xj-xmedi,boxxsize)
yj=boxshift(yj-ymedi,boxysize)
zj=boxshift(zj-zmedi,boxzsize)
if (energy_dec) then
write (iout,'(a6,2i5,0pf7.3,f7.3)') 'evdw1',i,j,evdwij,sss
endif
- evdw1=evdw1+evdwij*sss
+ evdw1=evdw1+evdwij*sss*faclipij2
if (energy_dec) write (iout,'(a10,2i5,0pf7.3)')
& 'evdw1_sum',i,j,evdw1
C
C Calculate contributions to the Cartesian gradient.
C
- facvdw=-6*rrmij*(ev1+evdwij)*sss
- ggg(1)=facvdw*xj+sssgrad*rmij*evdwij*xj/rpp(iteli,itelj)
- ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj/rpp(iteli,itelj)
- ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj/rpp(iteli,itelj)
+ facvdw=(-6*rrmij*(ev1+evdwij)*sss+sssgrad*rmij*evdwij/
+ & rpp(iteli,itelj))*faclipij2
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
C ggg(1)=facvdw*xj
C ggg(2)=facvdw*yj
C ggg(3)=facvdw*zj
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
endif
c enddo ! j
enddo ! i
projects / unres.git / blobdiff