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)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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
c endif
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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
- evdw=0.0D0
c if (icall.eq.0) then
c lprn=.true.
c else
c endif
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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)
& xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
double precision subchap,sss1,sssgrad1
- double precision boxshift
+ double precision boxshift,rij1
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
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ if (energy_dec)
+ & write(2,*) "g_listscsc_start_long,g_listscsc_end_long",
+ & g_listscsc_start_long,g_listscsc_end_long
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss1=sscale(1.0d0/rij,r_cut_int)
+ rij1=1.0d0/rij
+c sss1=sscale(1.0d0/rij,r_cut_int)
+ sss1=sscale(rij1,r_cut_int)
if (sss1.eq.0.0d0) cycle
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
+ rij1=rij1/sigmaii(itypi,itypj)
+ sss=sscale(rij1,r_cut_respa)
+c sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
if (sss.lt.1.0d0) then
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
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
C assuming the Gay-Berne potential of interaction.
C
implicit none
+ include 'mpif.h'
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
include "COMMON.SPLITELE"
+ include 'COMMON.TIME1'
logical lprn
double precision evdw
integer itypi,itypj,itypi1,iint,ind,ikont
& sslipj,ssgradlipj,ssgradlipi,sig,rij_shift,faclip
double precision dist,sscale,sscagrad,sscagradlip,sscalelip
double precision boxshift
+ double precision time01
+c time01=MPI_Wtime()
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
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ if (energy_dec)
+ & write(2,*) "g_listscsc_start_short,g_listscsc_end_short",
+ & g_listscsc_start_short,g_listscsc_end_short
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
- sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
if (sss.gt.0.0d0) then
+ sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
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
c enddo ! j
c enddo ! iint
enddo ! i
+c time_evdw_short=time_evdw_short+MPI_Wtime()-time01
c write (iout,*) "Number of loop steps in EGB:",ind
cccc energy_dec=.false.
return
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
- evdw=0.0D0
lprn=.false.
c if (icall.eq.0) lprn=.true.
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_long,g_listscsc_end_long
+ i=newcontlisti_long(ikont)
+ j=newcontlistj_long(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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
- evdw=0.0D0
lprn=.false.
c if (icall.eq.0) lprn=.true.
ind=0
c do i=iatsc_s,iatsc_e
- do ikont=g_listscsc_start,g_listscsc_end
- i=newcontlisti(ikont)
- j=newcontlistj(ikont)
+ do ikont=g_listscsc_start_short,g_listscsc_end_short
+ i=newcontlisti_short(ikont)
+ j=newcontlistj_short(ikont)
itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
itypi1=iabs(itype(i+1))
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
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
c do i=iatel_s,iatel_e
+ if (energy_dec)
+ & write(iout,*) "g_listpp_start,g_listpp_end",
+ & g_listpp_start,g_listpp_end
do ikont=g_listpp_start,g_listpp_end
i=newcontlistppi(ikont)
j=newcontlistppj(ikont)
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,
c & " iatel_e_vdw",iatel_e_vdw
c call flush(iout)
c do i=iatel_s_vdw,iatel_e_vdw
- do ikont=g_listpp_vdw_start,g_listpp_vdw_end
- i=newcontlistpp_vdwi(ikont)
- j=newcontlistpp_vdwj(ikont)
+ if (energy_dec)
+ & write(iout,*) "g_listpp_vdw_start_short,g_listpp_vdw_end_short",
+ & g_listpp_vdw_start_short,g_listpp_vdw_end_short
+ do ikont=g_listpp_vdw_start_short,g_listpp_vdw_end_short
+ i=newcontlistpp_vdwi_short(ikont)
+ j=newcontlistpp_vdwj_short(ikont)
if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1) cycle
dxi=dc(1,i)
dyi=dc(2,i)
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
c & write (iout,*) 'ESCP_LONG iatscp_s=',iatscp_s,
c & ' iatscp_e=',iatscp_e
c do i=iatscp_s,iatscp_e
- do ikont=g_listscp_start,g_listscp_end
- i=newcontlistscpi(ikont)
- j=newcontlistscpj(ikont)
+ if (energy_dec)
+ & write(iout,*)"g_listscp_start_long,g_listscp_end_long",
+ & g_listscp_start_long,g_listscp_end_long
+ do ikont=g_listscp_start_long,g_listscp_end_long
+ i=newcontlistscpi_long(ikont)
+ j=newcontlistscpj_long(ikont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
double precision ggg(3)
double precision sscale,sscagrad
double precision boxshift
+ integer ikont
evdw2=0.0D0
evdw2_14=0.0d0
cd print '(a)','Enter ESCP'
c if (lprint_short)
c & write (iout,*) 'ESCP_SHORT iatscp_s=',iatscp_s,
c & ' iatscp_e=',iatscp_e
- if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
- do i=iatscp_s,iatscp_e
+c if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
+ if (energy_dec)
+ & write(iout,*) "g_listscp_start_short,g_listscp_end_short",
+ & g_listscp_start_short,g_listscp_end_short
+ do ikont=g_listscp_start_short,g_listscp_end_short
+ i=newcontlistscpi_short(ikont)
+ j=newcontlistscpj_short(ikont)
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
iteli=itel(i)
xi=0.5D0*(c(1,i)+c(1,i+1))
c if (lprint_short)
c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
c & " nscp_gr",nscp_gr(i)
- do iint=1,nscp_gr(i)
-
- do j=iscpstart(i,iint),iscpend(i,iint)
+c do iint=1,nscp_gr(i)
+c
+c do j=iscpstart(i,iint),iscpend(i,iint)
itypj=iabs(itype(j))
c if (lprint_short)
c & write (iout,*) "j",j," itypj",itypj
gvdwc_scp(k,j)=gvdwc_scp(k,j)+ggg(k)
enddo
endif
- enddo
+c enddo
- enddo ! iint
+c enddo ! iint
enddo ! i
do i=1,nct
do j=1,3
projects / unres.git / blobdiff