C-----------------------------------------------------------------------
double precision function sscalelip(r)
+ implicit none
double precision r,gamm
include "COMMON.SPLITELE"
C if(r.lt.r_cut-rlamb) then
end
C-----------------------------------------------------------------------
double precision function sscagradlip(r)
+ implicit none
double precision r,gamm
include "COMMON.SPLITELE"
C if(r.lt.r_cut-rlamb) then
end
C-----------------------------------------------------------------------
- double precision function sscale(r)
- double precision r,gamm
+ double precision function sscale(r,r_cut)
+ implicit none
+ double precision r,r_cut,gamm
include "COMMON.SPLITELE"
if(r.lt.r_cut-rlamb) then
sscale=1.0d0
return
end
C-----------------------------------------------------------------------
-C-----------------------------------------------------------------------
- double precision function sscagrad(r)
- double precision r,gamm
+ double precision function sscagrad(r,r_cut)
+ implicit none
+ double precision r,r_cut,gamm
include "COMMON.SPLITELE"
if(r.lt.r_cut-rlamb) then
sscagrad=0.0d0
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the LJ potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
- parameter (accur=1.0d-10)
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.LOCAL'
include 'COMMON.SBRIDGE'
include 'COMMON.NAMES'
include 'COMMON.IOUNITS'
- include 'COMMON.CONTACTS'
- dimension gg(3)
+ include "COMMON.SPLITELE"
+c include 'COMMON.CONTACTS'
+ double precision gg(3)
+ double precision evdw,evdwij
+ integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+ double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
+ & sigij,r0ij,rcut,sss1,sssgrad1,sqrij
+ double precision sscale,sscagrad
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
cd & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
rij=xj*xj+yj*yj+zj*zj
- sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
+ sqrij=dsqrt(rrij)
+ eps0ij=eps(itypi,itypj)
+ sss1=sscale(sqrij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(sqrij,r_cut_int)
+ sssgrad=
+ & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
+ sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
if (sss.lt.1.0d0) then
rrij=1.0D0/rij
- eps0ij=eps(itypi,itypj)
fac=rrij**expon2
e1=fac*fac*aa
e2=fac*bb
evdwij=e1+e2
- evdw=evdw+(1.0d0-sss)*evdwij
+ evdw=evdw+(1.0d0-sss)*sss1*evdwij/sqrij/expon
C
C Calculate the components of the gradient in DC and X
C
- fac=-rrij*(e1+evdwij)*(1.0d0-sss)
+ fac=-rrij*(e1+evdwij)*(1.0d0-sss)*sss1
+ & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
+ & +(1.0d0-sss)*sssgrad1)/sqrij
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the LJ potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
- parameter (accur=1.0d-10)
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.LOCAL'
include 'COMMON.SBRIDGE'
include 'COMMON.NAMES'
include 'COMMON.IOUNITS'
- include 'COMMON.CONTACTS'
- dimension gg(3)
+ include "COMMON.SPLITELE"
+c include 'COMMON.CONTACTS'
+ double precision gg(3)
+ double precision evdw,evdwij
+ integer i,j,k,itypi,itypj,itypi1,num_conti,iint
+ double precision xi,yi,zi,xj,yj,zj,rij,eps0ij,fac,e1,e2,rrij,
+ & sigij,r0ij,rcut,sqrij,sss1,sssgrad1
+ double precision sscale,sscagrad
c write(iout,*)'Entering ELJ nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
cd write (iout,*) 'i=',i,' iint=',iint,' istart=',istart(i,iint),
cd & 'iend=',iend(i,iint)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
zj=c(3,nres+j)-zi
C Change 12/1/95 to calculate four-body interactions
rij=xj*xj+yj*yj+zj*zj
- sss=sscale(dsqrt(rij)/sigma(itypi,itypj))
+ sqrij=dsqrt(rij)
+ sss=sscale(sqrij/sigma(itypi,itypj),r_cut_respa)
if (sss.gt.0.0d0) then
+ sssgrad=
+ & sscagrad(sqrij/sigma(itypi,itypj),r_cut_respa)
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
C
C Calculate the components of the gradient in DC and X
C
- fac=-rrij*(e1+evdwij)*sss
+ fac=-rrij*(e1+evdwij)*sss+evdwij*sssgrad/sqrij/expon
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the LJK potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
- dimension gg(3)
+ include "COMMON.SPLITELE"
+ double precision gg(3)
+ double precision evdw,evdwij
+ integer i,j,k,itypi,itypj,itypi1,iint
+ 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
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
e_augm=augm(itypi,itypj)*fac_augm
r_inv_ij=dsqrt(rrij)
rij=1.0D0/r_inv_ij
- sss=sscale(rij/sigma(itypi,itypj))
+ sss1=sscale(rij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sssgrad1=sscagrad(rij,r_cut_int)
+ sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
if (sss.lt.1.0d0) then
+ sssgrad=
+ & 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
e1=fac*fac*aa
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)
- evdw=evdw+(1.0d0-sss)*evdwij
+ evdw=evdw+(1.0d0-sss)*sss1*evdwij
C
C Calculate the components of the gradient in DC and X
C
fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
- fac=fac*(1.0d0-sss)
+ fac=fac*(1.0d0-sss)*sss1
+ & +evdwij*(-sss1*sssgrad/sigma(itypi,itypj)
+ & +(1.0d0-sss)*sssgrad1)*r_inv_ij/expon
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.NAMES'
- dimension gg(3)
+ include "COMMON.SPLITELE"
+ double precision gg(3)
+ double precision evdw,evdwij
+ integer i,j,k,itypi,itypj,itypi1,iint
+ 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
c print *,'Entering ELJK nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
C
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
xj=c(1,nres+j)-xi
yj=c(2,nres+j)-yi
e_augm=augm(itypi,itypj)*fac_augm
r_inv_ij=dsqrt(rrij)
rij=1.0D0/r_inv_ij
- sss=sscale(rij/sigma(itypi,itypj))
+ sss=sscale(rij/sigma(itypi,itypj),r_cut_respa)
if (sss.gt.0.0d0) then
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
C Calculate the components of the gradient in DC and X
C
fac=-2.0D0*rrij*e_augm-r_inv_ij*r_shift_inv*(e1+e1+e2)
+ & +evdwij*sssgrad/sigma(itypi,itypj)*r_inv_ij/expon
fac=fac*sss
gg(1)=xj*fac
gg(2)=yj*fac
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the Berne-Pechukas potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include "COMMON.SPLITELE"
+ integer icall
common /srutu/ icall
+ double precision evdw
+ integer itypi,itypj,itypi1,iint,ind
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
+ double precision sss1,sssgrad1
+ double precision sscale,sscagrad
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
c endif
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
-
+ sss1=sscale(1.0d0/rij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
if (sss.lt.1.0d0) then
-
+ sssgrad=
+ & sscagrad(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
+ sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
C Calculate the angle-dependent terms of energy & contributions to derivatives.
call sc_angular
C Calculate whole angle-dependent part of epsilon and contributions
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*(1.0d0-sss)
+ evdw=evdw+evdwij*(1.0d0-sss)*sss1
if (lprn) then
sigm=dabs(aa/bb)**(1.0D0/6.0D0)
epsi=bb**2/aa
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)
sigder=fac/sigsq
- fac=rrij*fac
+ fac=(fac+evdwij*(sss1/(1.0d0-sss)*sssgrad/
+ & sigmaii(itypi,itypj)+(1.0d0-sss)/sss1*sssgrad1))*rij
C Calculate radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
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)
+ call sc_grad_scale((1.0d0-sss)*sss1)
endif
enddo ! j
enddo ! iint
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include "COMMON.SPLITELE"
+ integer icall
common /srutu/ icall
+ double precision evdw
+ integer itypi,itypj,itypi1,iint,ind
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
+ double precision sscale,sscagrad
c double precision rrsave(maxdim)
logical lprn
evdw=0.0D0
c endif
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
if (sss.gt.0.0d0) then
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)
sigder=fac/sigsq
- fac=rrij*fac
+ fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rrij
C Calculate radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the Gay-Berne potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
+ include "COMMON.SPLITELE"
logical lprn
integer xshift,yshift,zshift
+ double precision evdw
+ integer itypi,itypj,itypi1,iint,ind
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
+ double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
+ & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+ double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+ double precision subchap,sss1,sssgrad1
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
c if (icall.eq.0) lprn=.false.
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
xj=c(1,nres+j)
yj=c(2,nres+j)
zj=c(3,nres+j)
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- 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
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ 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
+ if (zj.lt.buflipbot) then
C what fraction I am in
- fracinbuf=1.0d0-
- & ((positi-bordlipbot)/lipbufthick)
+ fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-positi)/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
-
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
- dxj=dc_norm(1,nres+j)
- dyj=dc_norm(2,nres+j)
- dzj=dc_norm(3,nres+j)
- rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
- rij=dsqrt(rrij)
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
- sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
- if (sss.lt.1.0d0) then
-
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ else if (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/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
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if (dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+ dxj=dc_norm(1,nres+j)
+ dyj=dc_norm(2,nres+j)
+ 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)
+ if (sss1.eq.0.0d0) cycle
+ 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.
+ sssgrad=
+ & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
+ sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
call sc_angular
sigsq=1.0D0/sigsq
sig=sig0ij*dsqrt(sigsq)
c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*(1.0d0-sss)
+ evdw=evdw+evdwij*(1.0d0-sss)*sss1
if (lprn) then
sigm=dabs(aa/bb)**(1.0D0/6.0D0)
epsi=bb**2/aa
& evdwij
endif
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'evdw',i,j,evdwij
+ if (energy_dec) write (iout,'(a6,2i5,4f10.5)')
+ & 'evdw',i,j,rij,sss,sss1,evdwij
C Calculate gradient components.
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
- fac=rij*fac
- fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
+ fac=(fac+evdwij*(-sss1*sssgrad/(1.0d0-sss)
+ & /sigmaii(itypi,itypj)+(1.0d0-sss)*sssgrad1/sss1))*rij
c fac=0.0d0
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg_lipi(3)=ssgradlipi*evdwij
gg_lipj(3)=ssgradlipj*evdwij
C Calculate angular part of the gradient.
- call sc_grad_scale(1.0d0-sss)
+ call sc_grad_scale((1.0d0-sss)*sss1)
endif
enddo ! j
enddo ! iint
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the Gay-Berne potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
include 'COMMON.CONTROL'
+ include "COMMON.SPLITELE"
logical lprn
integer xshift,yshift,zshift
+ double precision evdw
+ integer itypi,itypj,itypi1,iint,ind
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,xi,yi,zi
+ double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
+ & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+ double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+ double precision subchap
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
c if (icall.eq.0) lprn=.false.
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
xj=c(1,nres+j)
yj=c(2,nres+j)
zj=c(3,nres+j)
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- 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
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ 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
+ if (zj.lt.buflipbot) then
C what fraction I am in
- fracinbuf=1.0d0-
- & ((positi-bordlipbot)/lipbufthick)
+ fracinbuf=1.0d0-((zi-bordlipbot)/lipbufthick)
C lipbufthick is thickenes of lipid buffore
- sslipj=sscalelip(fracinbuf)
- ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
- elseif (zi.gt.bufliptop) then
- fracinbuf=1.0d0-((bordliptop-positi)/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
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- subchap=1
- endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/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
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
- sssgrad=sscagrad((1.0d0/rij)/sigmaii(itypi,itypj))
+ 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
C Calculate angle-dependent terms of energy and contributions to their
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
- fac=rij*fac
- fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
+ fac=(fac+evdwij*sssgrad/sss/sigmaii(itypi,itypj))*rij
c fac=0.0d0
C Calculate the radial part of the gradient
gg(1)=xj*fac
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include "COMMON.SPLITELE"
+ integer icall
common /srutu/ icall
logical lprn
+ integer itypi,itypj,itypi1,iint,ind
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
+ & xi,yi,zi,fac_augm,e_augm
+ double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
+ & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+ double precision dist,sscale,sscagrad,sscagradlip,sscalelip
+ double precision sss1,sssgrad1
evdw=0.0D0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
c if (icall.eq.0) lprn=.true.
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss1=sscale(1.0d0/rij,r_cut_int)
+ if (sss1.eq.0.0d0) cycle
if (sss.lt.1.0d0) then
+ sssgrad=
+ & sscagrad((1.0d0/rij)/sigmaii(itypi,itypj),r_cut_respa)
+ sssgrad1=sscagrad(1.0d0/rij,r_cut_int)
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
fac_augm=rrij**expon
e_augm=augm(itypi,itypj)*fac_augm
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
+ evdw=evdw+(evdwij+e_augm)*sss1*(1.0d0-sss)
if (lprn) then
sigm=dabs(aa/bb)**(1.0D0/6.0D0)
epsi=bb**2/aa
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac-2*expon*rrij*e_augm
+ fac=fac+(evdwij+e_augm)*
+ & (-sss1*sssgrad/(1.0d0-sss)/sigmaii(itypi,itypj)
+ & +(1.0d0-sss)*sssgrad1/sss1)*rij
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
C Calculate angular part of the gradient.
- call sc_grad_scale(1.0d0-sss)
+ call sc_grad_scale((1.0d0-sss)*sss1)
endif
enddo ! j
enddo ! iint
C This subroutine calculates the interaction energy of nonbonded side chains
C assuming the Gay-Berne-Vorobjev potential of interaction.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
include 'COMMON.IOUNITS'
include 'COMMON.CALC'
+ include "COMMON.SPLITELE"
+ integer icall
common /srutu/ icall
logical lprn
+ integer itypi,itypj,itypi1,iint,ind
+ double precision eps0ij,epsi,sigm,fac,e1,e2,rrij,r0ij,
+ & xi,yi,zi,fac_augm,e_augm
+ double precision evdw
+ double precision fracinbuf,sslipi,evdwij_przed_tri,sig0ij,
+ & sslipj,ssgradlipj,ssgradlipi,dist_init,xj_safe,yj_safe,zj_safe,
+ & xj_temp,yj_temp,zj_temp,dist_temp,sig,rij_shift,faclip
+ double precision dist,sscale,sscagrad,sscagradlip,sscalelip
evdw=0.0D0
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
c if (icall.eq.0) lprn=.true.
ind=0
do i=iatsc_s,iatsc_e
- itypi=itype(i)
+ itypi=iabs(itype(i))
if (itypi.eq.ntyp1) cycle
- itypi1=itype(i+1)
+ itypi1=iabs(itype(i+1))
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
do iint=1,nint_gr(i)
do j=istart(i,iint),iend(i,iint)
ind=ind+1
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
c dscj_inv=dsc_inv(itypj)
dscj_inv=vbld_inv(j+nres)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
- sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)),r_cut_respa)
if (sss.gt.0.0d0) then
e1=e1*eps1*eps2rt**2*eps3rt**2
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
- fac=rij*fac-2*expon*rrij*e_augm
+ fac=rij*fac-2*expon*rrij*e_augm+
+ & (evdwij+e_augm)*sssgrad/sigmaii(itypi,itypj)/sss*rij
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
include 'COMMON.DERIV'
include 'COMMON.CALC'
include 'COMMON.IOUNITS'
+ include "COMMON.SPLITELE"
double precision dcosom1(3),dcosom2(3)
double precision scalfac
eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
+#ifdef FOURBODY
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+#endif
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
include 'COMMON.SHIELD'
+ include "COMMON.SPLITELE"
dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
eello_turn3=0.0d0
eello_turn4=0.0d0
ind=0
+#ifdef FOURBODY
do i=1,nres
num_cont_hb(i)=0
enddo
+#endif
cd print '(a)','Enter EELEC'
cd write (iout,*) 'iatel_s=',iatel_s,' iatel_e=',iatel_e
do i=1,nres
num_conti=0
call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo
do i=iturn4_start,iturn4_end
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
if (ymedi.lt.0) ymedi=ymedi+boxysize
zmedi=mod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
+#ifdef FOURBODY
num_conti=num_cont_hb(i)
+#endif
call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
& call eturn4(i,eello_turn4)
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo ! i
c
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
if (ymedi.lt.0) ymedi=ymedi+boxysize
zmedi=mod(zmedi,boxzsize)
if (zmedi.lt.0) zmedi=zmedi+boxzsize
-c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
+c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend
+#ifdef FOURBODY
num_conti=num_cont_hb(i)
+#endif
do j=ielstart(i),ielend(i)
if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
C & .or.itype(j+2).eq.ntyp1
&) cycle
call eelecij_scale(i,j,ees,evdw1,eel_loc)
enddo ! j
+#ifdef FOURBODY
num_cont_hb(i)=num_conti
+#endif
enddo ! i
c write (iout,*) "Number of loop steps in EELEC:",ind
cd do i=1,nres
end
C-------------------------------------------------------------------------------
subroutine eelecij_scale(i,j,ees,evdw1,eel_loc)
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
#ifdef MPI
include "mpif.h"
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
+#ifdef FOURBODY
include 'COMMON.CONTACTS'
+ include 'COMMON.CONTMAT'
+#endif
+ include 'COMMON.CORRMAT'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
include 'COMMON.TIME1'
include 'COMMON.SHIELD'
+ include "COMMON.SPLITELE"
integer xshift,yshift,zshift
- dimension ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
+ double precision ggg(3),gggp(3),gggm(3),erij(3),dcosb(3),dcosg(3),
& erder(3,3),uryg(3,3),urzg(3,3),vryg(3,3),vrzg(3,3)
double precision acipa(2,2),agg(3,4),aggi(3,4),aggi1(3,4),
& aggj(3,4),aggj1(3,4),a_temp(2,2),muij(4),gmuij1(4),gmuji1(4),
& gmuij2(4),gmuji2(4)
+ integer j1,j2,num_conti
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
+ integer k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap,ind,itypi,itypj
+ integer ilist,iresshield
+ double precision rlocshield
+ double precision ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
+ double precision ees,evdw1,eel_loc,aaa,bbb,ael3i
+ double precision dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,
+ & rij,r3ij,r6ij,cosa,cosb,cosg,fac,ev1,ev2,fac3,fac4,
+ & evdwij,el1,el2,eesij,ees0ij,facvdw,facel,fac1,ecosa,
+ & ecosb,ecosg,ury,urz,vry,vrz,facr,a22der,a23der,a32der,
+ & a33der,eel_loc_ij,cosa4,wij,cosbg1,cosbg2,ees0pij,
+ & ees0pij1,ees0mij,ees0mij1,fac3p,ees0mijp,ees0pijp,
+ & ecosa1,ecosb1,ecosg1,ecosa2,ecosb2,ecosg2,ecosap,ecosbp,
+ & ecosgp,ecosam,ecosbm,ecosgm,ghalf,geel_loc_ij,geel_loc_ji,
+ & dxi,dyi,dzi,a22,a23,a32,a33
+ double precision dist_init,xmedi,ymedi,zmedi,xj_safe,yj_safe,
+ & zj_safe,xj_temp,yj_temp,zj_temp,dist_temp,dx_normi,dy_normi,
+ & dz_normi,aux
+ double precision sss1,sssgrad1
+ double precision sscale,sscagrad
+ double precision scalar
+
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
double precision scal_el /1.0d0/
c time00=MPI_Wtime()
cd write (iout,*) "eelecij",i,j
C print *,"WCHODZE2"
- ind=ind+1
- iteli=itel(i)
- itelj=itel(j)
- if (j.eq.i+2 .and. itelj.eq.2) iteli=2
- aaa=app(iteli,itelj)
- bbb=bpp(iteli,itelj)
- ael6i=ael6(iteli,itelj)
- ael3i=ael3(iteli,itelj)
- dxj=dc(1,j)
- dyj=dc(2,j)
- dzj=dc(3,j)
- dx_normj=dc_norm(1,j)
- dy_normj=dc_norm(2,j)
- dz_normj=dc_norm(3,j)
- xj=c(1,j)+0.5D0*dxj
- yj=c(2,j)+0.5D0*dyj
- zj=c(3,j)+0.5D0*dzj
- xj=mod(xj,boxxsize)
- if (xj.lt.0) xj=xj+boxxsize
- yj=mod(yj,boxysize)
- if (yj.lt.0) yj=yj+boxysize
- zj=mod(zj,boxzsize)
- if (zj.lt.0) zj=zj+boxzsize
+ ind=ind+1
+ iteli=itel(i)
+ itelj=itel(j)
+ if (j.eq.i+2 .and. itelj.eq.2) iteli=2
+ aaa=app(iteli,itelj)
+ bbb=bpp(iteli,itelj)
+ ael6i=ael6(iteli,itelj)
+ ael3i=ael3(iteli,itelj)
+ dxj=dc(1,j)
+ dyj=dc(2,j)
+ dzj=dc(3,j)
+ dx_normj=dc_norm(1,j)
+ dy_normj=dc_norm(2,j)
+ dz_normj=dc_norm(3,j)
+ xj=c(1,j)+0.5D0*dxj
+ yj=c(2,j)+0.5D0*dyj
+ zj=c(3,j)+0.5D0*dzj
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
xj_safe=xj
yj_safe=yj
zj_temp=zj
isubchap=1
endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
- rij=xj*xj+yj*yj+zj*zj
- rrmij=1.0D0/rij
- rij=dsqrt(rij)
- rmij=1.0D0/rij
+ rij=xj*xj+yj*yj+zj*zj
+ rrmij=1.0D0/rij
+ rij=dsqrt(rij)
+ rmij=1.0D0/rij
c For extracting the short-range part of Evdwpp
- sss=sscale(rij/rpp(iteli,itelj))
- sssgrad=sscagrad(rij/rpp(iteli,itelj))
- r3ij=rrmij*rmij
- r6ij=r3ij*r3ij
- cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
- cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
- cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
- fac=cosa-3.0D0*cosb*cosg
- ev1=aaa*r6ij*r6ij
+ sss1=sscale(rij,r_cut_int)
+ if (sss1.eq.0.0d0) return
+ sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
+ sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
+ sssgrad1=sscagrad(rij,r_cut_int)
+ r3ij=rrmij*rmij
+ r6ij=r3ij*r3ij
+ cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
+ cosb=(xj*dx_normi+yj*dy_normi+zj*dz_normi)*rmij
+ cosg=(xj*dx_normj+yj*dy_normj+zj*dz_normj)*rmij
+ fac=cosa-3.0D0*cosb*cosg
+ ev1=aaa*r6ij*r6ij
c 4/26/02 - AL scaling down 1,4 repulsive VDW interactions
- if (j.eq.i+2) ev1=scal_el*ev1
- ev2=bbb*r6ij
- fac3=ael6i*r6ij
- fac4=ael3i*r3ij
- evdwij=ev1+ev2
- if (shield_mode.eq.0) then
- fac_shield(i)=1.0
- fac_shield(j)=1.0
- endif
- el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
- el2=fac4*fac
- el1=el1*fac_shield(i)**2*fac_shield(j)**2
- el2=el2*fac_shield(i)**2*fac_shield(j)**2
- eesij=el1+el2
+ if (j.eq.i+2) ev1=scal_el*ev1
+ ev2=bbb*r6ij
+ fac3=ael6i*r6ij
+ fac4=ael3i*r3ij
+ evdwij=ev1+ev2
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+ endif
+ el1=fac3*(4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg))
+ el2=fac4*fac
+ el1=el1*fac_shield(i)**2*fac_shield(j)**2
+ el2=el2*fac_shield(i)**2*fac_shield(j)**2
+ 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
- evdw1=evdw1+evdwij*(1.0d0-sss)
+ ees0ij=4.0D0+fac*fac-3.0D0*(cosb*cosb+cosg*cosg)
+ ees=ees+eesij*sss1
+ evdw1=evdw1+evdwij*(1.0d0-sss)*sss1
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,f7.3)') 'evdw1',i,j,evdwij,sss
- write (iout,'(a6,2i5,0pf7.3)') 'ees',i,j,eesij
- endif
+ 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
+ endif
C
C Calculate contributions to the Cartesian gradient.
C
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)
- facel=-3*rrmij*(el1+eesij)
- fac1=fac
- erij(1)=xj*rmij
- erij(2)=yj*rmij
- erij(3)=zj*rmij
+ facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
+c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ facel=-3*rrmij*(el1+eesij)*sss1
+c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ fac1=fac
+ erij(1)=xj*rmij
+ erij(2)=yj*rmij
+ erij(3)=zj*rmij
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
- ggg(1)=facel*xj
- ggg(2)=facel*yj
- ggg(3)=facel*zj
- if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ aux=facel+sssgrad1*(1.0d0-sss)*eesij*rmij
+c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(1)=aux*xj
+ ggg(2)=aux*yj
+ ggg(3)=aux*zj
+c ggg(1)=facel*xj
+c ggg(2)=facel*yj
+c ggg(3)=facel*zj
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
& (shield_mode.gt.0)) then
C print *,i,j
- do ilist=1,ishield_list(i)
- iresshield=shield_list(ilist,i)
- do k=1,3
- rlocshield=grad_shield_side(k,ilist,i)*eesij/fac_shield(i)
- & *2.0
+ do ilist=1,ishield_list(i)
+ iresshield=shield_list(ilist,i)
+ do k=1,3
+ rlocshield=grad_shield_side(k,ilist,i)*eesij*sss1
+ & /fac_shield(i)*2.0*sss1
gshieldx(k,iresshield)=gshieldx(k,iresshield)+
& rlocshield
- & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)*2.0
+ & +grad_shield_loc(k,ilist,i)*eesij*sss1/fac_shield(i)*2.0
+ & *sss1
gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
C gshieldc_loc(k,iresshield)=gshieldc_loc(k,iresshield)
C & +grad_shield_loc(k,ilist,i)*eesij/fac_shield(i)
C
C enddo
C endif
- enddo
enddo
- do ilist=1,ishield_list(j)
- iresshield=shield_list(ilist,j)
- do k=1,3
+ enddo
+ do ilist=1,ishield_list(j)
+ iresshield=shield_list(ilist,j)
+ do k=1,3
rlocshield=grad_shield_side(k,ilist,j)*eesij/fac_shield(j)
- & *2.0
+ & *2.0*sss1
gshieldx(k,iresshield)=gshieldx(k,iresshield)+
- & rlocshield
- & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0
+ & rlocshield
+ & +grad_shield_loc(k,ilist,j)*eesij/fac_shield(j)*2.0*sss1
gshieldc(k,iresshield-1)=gshieldc(k,iresshield-1)+rlocshield
- enddo
enddo
+ enddo
- do k=1,3
- gshieldc(k,i)=gshieldc(k,i)+
- & grad_shield(k,i)*eesij/fac_shield(i)*2.0
- gshieldc(k,j)=gshieldc(k,j)+
- & grad_shield(k,j)*eesij/fac_shield(j)*2.0
- gshieldc(k,i-1)=gshieldc(k,i-1)+
- & grad_shield(k,i)*eesij/fac_shield(i)*2.0
- gshieldc(k,j-1)=gshieldc(k,j-1)+
- & grad_shield(k,j)*eesij/fac_shield(j)*2.0
+ do k=1,3
+ gshieldc(k,i)=gshieldc(k,i)+
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
+ gshieldc(k,j)=gshieldc(k,j)+
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
+ gshieldc(k,i-1)=gshieldc(k,i-1)+
+ & grad_shield(k,i)*eesij/fac_shield(i)*2.0*sss1
+ gshieldc(k,j-1)=gshieldc(k,j-1)+
+ & grad_shield(k,j)*eesij/fac_shield(j)*2.0*sss1
- enddo
- endif
+ enddo
+ endif
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,j)=gelc(k,j)+ghalf
c enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
- do k=1,3
- gelc_long(k,j)=gelc_long(k,j)+ggg(k)
- gelc_long(k,i)=gelc_long(k,i)-ggg(k)
- enddo
+ do k=1,3
+ gelc_long(k,j)=gelc_long(k,j)+ggg(k)
+ gelc_long(k,i)=gelc_long(k,i)-ggg(k)
+ enddo
+c gelc_long(3,i)=gelc_long(3,i)+
+c ssgradlipi*eesij/2.0d0*lipscale**2*sss1
*
* Loop over residues i+1 thru j-1.
*
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
- 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=facvdw+
+ & (-sss1*sssgrad/rpp(iteli,itelj)+(1.0d0-sss)*sssgrad1)*rmij*evdwij
+c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
c gvdwpp(k,j)=gvdwpp(k,j)+ghalf
c enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
- do k=1,3
- gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
- gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
- enddo
+ do k=1,3
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ enddo
*
* Loop over residues i+1 thru j-1.
*
cgrad enddo
cgrad enddo
#else
- facvdw=ev1+evdwij*(1.0d0-sss)
- facel=el1+eesij
- fac1=fac
- fac=-3*rrmij*(facvdw+facvdw+facel)
- erij(1)=xj*rmij
- erij(2)=yj*rmij
- erij(3)=zj*rmij
+ facvdw=-6*rrmij*(ev1+evdwij)*(1.0d0-sss)*sss1
+c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ facel=-3*rrmij*(el1+eesij)*sss1
+c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+
+c facvdw=ev1+evdwij*(1.0d0-sss)*sss1
+c facel=el1+eesij
+ fac1=fac
+ fac=-3*rrmij*(facvdw+facvdw+facel)
+ erij(1)=xj*rmij
+ erij(2)=yj*rmij
+ erij(3)=zj*rmij
*
* Radial derivatives. First process both termini of the fragment (i,j)
*
- ggg(1)=fac*xj
- ggg(2)=fac*yj
- ggg(3)=fac*zj
+ aux=fac+(sssgrad1*(1.0d0-sss)-sssgrad*sss1/rpp(iteli,itelj))
+ & *eesij*rmij
+c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
+ ggg(1)=aux*xj
+ ggg(2)=aux*yj
+ ggg(3)=axu*zj
+c ggg(1)=fac*xj
+c ggg(2)=fac*yj
+c ggg(3)=fac*zj
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
c gelc(k,j)=gelc(k,j)+ghalf
c enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
- do k=1,3
- gelc_long(k,j)=gelc(k,j)+ggg(k)
- gelc_long(k,i)=gelc(k,i)-ggg(k)
- enddo
+ do k=1,3
+ gelc_long(k,j)=gelc(k,j)+ggg(k)
+ gelc_long(k,i)=gelc(k,i)-ggg(k)
+ enddo
*
* Loop over residues i+1 thru j-1.
*
C ggg(1)=facvdw*xj
C ggg(2)=facvdw*yj
C ggg(3)=facvdw*zj
- 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)
- do k=1,3
- gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
- gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
- enddo
+ facvdw=facvdw
+ & (-sssgrad*sss1/rpp(iteli,itelj)+sssgrad1*(1.0d0-sss))*rmij*evdwij
+ ggg(1)=facvdw*xj
+ ggg(2)=facvdw*yj
+ ggg(3)=facvdw*zj
+ do k=1,3
+ gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
+ gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
+ enddo
#endif
*
* Angular part
*
- ecosa=2.0D0*fac3*fac1+fac4
- fac4=-3.0D0*fac4
- fac3=-6.0D0*fac3
- ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
- ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
- do k=1,3
- dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
- dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
- enddo
+ ecosa=2.0D0*fac3*fac1+fac4
+ fac4=-3.0D0*fac4
+ fac3=-6.0D0*fac3
+ ecosb=(fac3*(fac1*cosg+cosb)+cosg*fac4)
+ ecosg=(fac3*(fac1*cosb+cosg)+cosb*fac4)
+ do k=1,3
+ dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
+ dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
+ enddo
cd print '(2i3,2(3(1pd14.5),3x))',i,j,(dcosb(k),k=1,3),
cd & (dcosg(k),k=1,3)
- do k=1,3
- ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*
- & fac_shield(i)**2*fac_shield(j)**2
+ do k=1,3
+ ggg(k)=(ecosb*dcosb(k)+ecosg*dcosg(k))*sss1
+ & *fac_shield(i)**2*fac_shield(j)**2
+c & *((sslipi+sslipj)/2.0d0*lipscale**2+1.0d0)
- enddo
+ enddo
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gelc(k,i)=gelc(k,i)+ghalf
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
- do k=1,3
- 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))
- & *fac_shield(i)**2*fac_shield(j)**2
+ do k=1,3
+ 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
+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))
- & *fac_shield(i)**2*fac_shield(j)**2
- gelc_long(k,j)=gelc_long(k,j)+ggg(k)
- gelc_long(k,i)=gelc_long(k,i)-ggg(k)
- enddo
- IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
- & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
- & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
+ 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
+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)
+ enddo
+ IF (wel_loc.gt.0.0d0 .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0
+ & .or. wcorr6.gt.0.0d0 .or. wturn3.gt.0.0d0
+ & .or. wturn4.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
C
C 9/25/99 Mixed third-order local-electrostatic terms. The local-interaction
C energy of a peptide unit is assumed in the form of a second-order
C Macromolecules, 1974, 7, 797-806 for definition). This correlation terms
C are computed for EVERY pair of non-contiguous peptide groups.
C
- if (j.lt.nres-1) then
- j1=j+1
- j2=j-1
- else
- j1=j-1
- j2=j-2
- endif
- kkk=0
- do k=1,2
- do l=1,2
- kkk=kkk+1
- muij(kkk)=mu(k,i)*mu(l,j)
+ if (j.lt.nres-1) then
+ j1=j+1
+ j2=j-1
+ else
+ j1=j-1
+ j2=j-2
+ endif
+ kkk=0
+ do k=1,2
+ do l=1,2
+ kkk=kkk+1
+ muij(kkk)=mu(k,i)*mu(l,j)
#ifdef NEWCORR
- gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
+ gmuij1(kkk)=gtb1(k,i+1)*mu(l,j)
c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j)
- gmuij2(kkk)=gUb2(k,i)*mu(l,j)
- gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
+ gmuij2(kkk)=gUb2(k,i)*mu(l,j)
+ gmuji1(kkk)=mu(k,i)*gtb1(l,j+1)
c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i)
- gmuji2(kkk)=mu(k,i)*gUb2(l,j)
+ gmuji2(kkk)=mu(k,i)*gUb2(l,j)
#endif
- enddo
- enddo
+ enddo
+ enddo
cd write (iout,*) 'EELEC: i',i,' j',j
cd write (iout,*) 'j',j,' j1',j1,' j2',j2
cd write(iout,*) 'muij',muij
- ury=scalar(uy(1,i),erij)
- urz=scalar(uz(1,i),erij)
- vry=scalar(uy(1,j),erij)
- vrz=scalar(uz(1,j),erij)
- a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
- a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
- a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
- a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
- fac=dsqrt(-ael6i)*r3ij
- a22=a22*fac
- a23=a23*fac
- a32=a32*fac
- a33=a33*fac
+ ury=scalar(uy(1,i),erij)
+ urz=scalar(uz(1,i),erij)
+ vry=scalar(uy(1,j),erij)
+ vrz=scalar(uz(1,j),erij)
+ a22=scalar(uy(1,i),uy(1,j))-3*ury*vry
+ a23=scalar(uy(1,i),uz(1,j))-3*ury*vrz
+ a32=scalar(uz(1,i),uy(1,j))-3*urz*vry
+ a33=scalar(uz(1,i),uz(1,j))-3*urz*vrz
+ fac=dsqrt(-ael6i)*r3ij
+ a22=a22*fac
+ a23=a23*fac
+ a32=a32*fac
+ a33=a33*fac
cd write (iout,'(4i5,4f10.5)')
cd & i,itortyp(itype(i)),j,itortyp(itype(j)),a22,a23,a32,a33
cd write (iout,'(6f10.5)') (muij(k),k=1,4),fac,eel_loc_ij
cd write (iout,'(9f10.5/)')
cd & fac22,a22,fac23,a23,fac32,a32,fac33,a33,eel_loc_ij
C Derivatives of the elements of A in virtual-bond vectors
- call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
- do k=1,3
- uryg(k,1)=scalar(erder(1,k),uy(1,i))
- uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
- uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
- urzg(k,1)=scalar(erder(1,k),uz(1,i))
- urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
- urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
- vryg(k,1)=scalar(erder(1,k),uy(1,j))
- vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
- vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
- vrzg(k,1)=scalar(erder(1,k),uz(1,j))
- vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
- vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
- enddo
+ call unormderiv(erij(1),unmat(1,1),rmij,erder(1,1))
+ do k=1,3
+ uryg(k,1)=scalar(erder(1,k),uy(1,i))
+ uryg(k,2)=scalar(uygrad(1,k,1,i),erij(1))
+ uryg(k,3)=scalar(uygrad(1,k,2,i),erij(1))
+ urzg(k,1)=scalar(erder(1,k),uz(1,i))
+ urzg(k,2)=scalar(uzgrad(1,k,1,i),erij(1))
+ urzg(k,3)=scalar(uzgrad(1,k,2,i),erij(1))
+ vryg(k,1)=scalar(erder(1,k),uy(1,j))
+ vryg(k,2)=scalar(uygrad(1,k,1,j),erij(1))
+ vryg(k,3)=scalar(uygrad(1,k,2,j),erij(1))
+ vrzg(k,1)=scalar(erder(1,k),uz(1,j))
+ vrzg(k,2)=scalar(uzgrad(1,k,1,j),erij(1))
+ vrzg(k,3)=scalar(uzgrad(1,k,2,j),erij(1))
+ enddo
C Compute radial contributions to the gradient
- facr=-3.0d0*rrmij
- a22der=a22*facr
- a23der=a23*facr
- a32der=a32*facr
- a33der=a33*facr
- agg(1,1)=a22der*xj
- agg(2,1)=a22der*yj
- agg(3,1)=a22der*zj
- agg(1,2)=a23der*xj
- agg(2,2)=a23der*yj
- agg(3,2)=a23der*zj
- agg(1,3)=a32der*xj
- agg(2,3)=a32der*yj
- agg(3,3)=a32der*zj
- agg(1,4)=a33der*xj
- agg(2,4)=a33der*yj
- agg(3,4)=a33der*zj
+ facr=-3.0d0*rrmij
+ a22der=a22*facr
+ a23der=a23*facr
+ a32der=a32*facr
+ a33der=a33*facr
+ agg(1,1)=a22der*xj
+ agg(2,1)=a22der*yj
+ agg(3,1)=a22der*zj
+ agg(1,2)=a23der*xj
+ agg(2,2)=a23der*yj
+ agg(3,2)=a23der*zj
+ agg(1,3)=a32der*xj
+ agg(2,3)=a32der*yj
+ agg(3,3)=a32der*zj
+ agg(1,4)=a33der*xj
+ agg(2,4)=a33der*yj
+ agg(3,4)=a33der*zj
C Add the contributions coming from er
- fac3=-3.0d0*fac
- do k=1,3
- agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
- agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
- agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
- agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
- enddo
- do k=1,3
+ fac3=-3.0d0*fac
+ do k=1,3
+ agg(k,1)=agg(k,1)+fac3*(uryg(k,1)*vry+vryg(k,1)*ury)
+ agg(k,2)=agg(k,2)+fac3*(uryg(k,1)*vrz+vrzg(k,1)*ury)
+ agg(k,3)=agg(k,3)+fac3*(urzg(k,1)*vry+vryg(k,1)*urz)
+ agg(k,4)=agg(k,4)+fac3*(urzg(k,1)*vrz+vrzg(k,1)*urz)
+ enddo
+ do k=1,3
C Derivatives in DC(i)
cgrad ghalf1=0.5d0*agg(k,1)
cgrad ghalf2=0.5d0*agg(k,2)
cgrad ghalf3=0.5d0*agg(k,3)
cgrad ghalf4=0.5d0*agg(k,4)
- aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
- & -3.0d0*uryg(k,2)*vry)!+ghalf1
- aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
- & -3.0d0*uryg(k,2)*vrz)!+ghalf2
- aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
- & -3.0d0*urzg(k,2)*vry)!+ghalf3
- aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
- & -3.0d0*urzg(k,2)*vrz)!+ghalf4
+ aggi(k,1)=fac*(scalar(uygrad(1,k,1,i),uy(1,j))
+ & -3.0d0*uryg(k,2)*vry)!+ghalf1
+ aggi(k,2)=fac*(scalar(uygrad(1,k,1,i),uz(1,j))
+ & -3.0d0*uryg(k,2)*vrz)!+ghalf2
+ aggi(k,3)=fac*(scalar(uzgrad(1,k,1,i),uy(1,j))
+ & -3.0d0*urzg(k,2)*vry)!+ghalf3
+ aggi(k,4)=fac*(scalar(uzgrad(1,k,1,i),uz(1,j))
+ & -3.0d0*urzg(k,2)*vrz)!+ghalf4
C Derivatives in DC(i+1)
- aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
- & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
- aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
- & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
- aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
- & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
- aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
- & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
+ aggi1(k,1)=fac*(scalar(uygrad(1,k,2,i),uy(1,j))
+ & -3.0d0*uryg(k,3)*vry)!+agg(k,1)
+ aggi1(k,2)=fac*(scalar(uygrad(1,k,2,i),uz(1,j))
+ & -3.0d0*uryg(k,3)*vrz)!+agg(k,2)
+ aggi1(k,3)=fac*(scalar(uzgrad(1,k,2,i),uy(1,j))
+ & -3.0d0*urzg(k,3)*vry)!+agg(k,3)
+ aggi1(k,4)=fac*(scalar(uzgrad(1,k,2,i),uz(1,j))
+ & -3.0d0*urzg(k,3)*vrz)!+agg(k,4)
C Derivatives in DC(j)
- aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
- & -3.0d0*vryg(k,2)*ury)!+ghalf1
- aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
- & -3.0d0*vrzg(k,2)*ury)!+ghalf2
- aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
- & -3.0d0*vryg(k,2)*urz)!+ghalf3
- aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
- & -3.0d0*vrzg(k,2)*urz)!+ghalf4
+ aggj(k,1)=fac*(scalar(uygrad(1,k,1,j),uy(1,i))
+ & -3.0d0*vryg(k,2)*ury)!+ghalf1
+ aggj(k,2)=fac*(scalar(uzgrad(1,k,1,j),uy(1,i))
+ & -3.0d0*vrzg(k,2)*ury)!+ghalf2
+ aggj(k,3)=fac*(scalar(uygrad(1,k,1,j),uz(1,i))
+ & -3.0d0*vryg(k,2)*urz)!+ghalf3
+ aggj(k,4)=fac*(scalar(uzgrad(1,k,1,j),uz(1,i))
+ & -3.0d0*vrzg(k,2)*urz)!+ghalf4
C Derivatives in DC(j+1) or DC(nres-1)
- aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
- & -3.0d0*vryg(k,3)*ury)
- aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
- & -3.0d0*vrzg(k,3)*ury)
- aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
- & -3.0d0*vryg(k,3)*urz)
- aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
- & -3.0d0*vrzg(k,3)*urz)
+ aggj1(k,1)=fac*(scalar(uygrad(1,k,2,j),uy(1,i))
+ & -3.0d0*vryg(k,3)*ury)
+ aggj1(k,2)=fac*(scalar(uzgrad(1,k,2,j),uy(1,i))
+ & -3.0d0*vrzg(k,3)*ury)
+ aggj1(k,3)=fac*(scalar(uygrad(1,k,2,j),uz(1,i))
+ & -3.0d0*vryg(k,3)*urz)
+ aggj1(k,4)=fac*(scalar(uzgrad(1,k,2,j),uz(1,i))
+ & -3.0d0*vrzg(k,3)*urz)
cgrad if (j.eq.nres-1 .and. i.lt.j-2) then
cgrad do l=1,4
cgrad aggj1(k,l)=aggj1(k,l)+agg(k,l)
cgrad enddo
cgrad endif
+ enddo
+ acipa(1,1)=a22
+ acipa(1,2)=a23
+ acipa(2,1)=a32
+ acipa(2,2)=a33
+ a22=-a22
+ a23=-a23
+ do l=1,2
+ do k=1,3
+ agg(k,l)=-agg(k,l)
+ aggi(k,l)=-aggi(k,l)
+ aggi1(k,l)=-aggi1(k,l)
+ aggj(k,l)=-aggj(k,l)
+ aggj1(k,l)=-aggj1(k,l)
enddo
- acipa(1,1)=a22
- acipa(1,2)=a23
- acipa(2,1)=a32
- acipa(2,2)=a33
+ enddo
+ if (j.lt.nres-1) then
a22=-a22
- a23=-a23
- do l=1,2
+ a32=-a32
+ do l=1,3,2
do k=1,3
agg(k,l)=-agg(k,l)
aggi(k,l)=-aggi(k,l)
aggj1(k,l)=-aggj1(k,l)
enddo
enddo
- if (j.lt.nres-1) then
- a22=-a22
- a32=-a32
- do l=1,3,2
- do k=1,3
- agg(k,l)=-agg(k,l)
- aggi(k,l)=-aggi(k,l)
- aggi1(k,l)=-aggi1(k,l)
- aggj(k,l)=-aggj(k,l)
- aggj1(k,l)=-aggj1(k,l)
- enddo
+ else
+ a22=-a22
+ a23=-a23
+ a32=-a32
+ a33=-a33
+ do l=1,4
+ do k=1,3
+ agg(k,l)=-agg(k,l)
+ aggi(k,l)=-aggi(k,l)
+ aggi1(k,l)=-aggi1(k,l)
+ aggj(k,l)=-aggj(k,l)
+ aggj1(k,l)=-aggj1(k,l)
enddo
- else
- a22=-a22
- a23=-a23
- a32=-a32
- a33=-a33
- do l=1,4
- do k=1,3
- agg(k,l)=-agg(k,l)
- aggi(k,l)=-aggi(k,l)
- aggi1(k,l)=-aggi1(k,l)
- aggj(k,l)=-aggj(k,l)
- aggj1(k,l)=-aggj1(k,l)
- enddo
- enddo
- endif
- ENDIF ! WCORR
- IF (wel_loc.gt.0.0d0) THEN
+ enddo
+ endif
+ ENDIF ! WCORR
+ IF (wel_loc.gt.0.0d0) THEN
C Contribution to the local-electrostatic energy coming from the i-j pair
- eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
- & +a33*muij(4)
-cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
+ eel_loc_ij=a22*muij(1)+a23*muij(2)+a32*muij(3)
+ & +a33*muij(4)
+cd write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij
- if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
- & 'eelloc',i,j,eel_loc_ij
+ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)')
+ & 'eelloc',i,j,eel_loc_ij
- if (shield_mode.eq.0) then
- fac_shield(i)=1.0
- fac_shield(j)=1.0
-C else
-C fac_shield(i)=0.4
-C fac_shield(j)=0.6
- endif
- eel_loc_ij=eel_loc_ij
- & *fac_shield(i)*fac_shield(j)
- eel_loc=eel_loc+eel_loc_ij
+ if (shield_mode.eq.0) then
+ fac_shield(i)=1.0
+ fac_shield(j)=1.0
+C else
+C fac_shield(i)=0.4
+C fac_shield(j)=0.6
+ endif
+ eel_loc_ij=eel_loc_ij
+ & *fac_shield(i)*fac_shield(j)*sss1
+ eel_loc=eel_loc+eel_loc_ij
- if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
+ if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and.
& (shield_mode.gt.0)) then
C print *,i,j
C & *2.0
gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
& rlocshield
- & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
+ & +grad_shield_loc(k,ilist,i)*eel_loc_ij/fac_shield(i)
gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
& +rlocshield
enddo
C & *2.0
gshieldx_ll(k,iresshield)=gshieldx_ll(k,iresshield)+
& rlocshield
- & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
+ & +grad_shield_loc(k,ilist,j)*eel_loc_ij/fac_shield(j)
gshieldc_ll(k,iresshield-1)=gshieldc_ll(k,iresshield-1)
& +rlocshield
& grad_shield(k,i)*eel_loc_ij/fac_shield(i)
gshieldc_ll(k,j-1)=gshieldc_ll(k,j-1)+
& grad_shield(k,j)*eel_loc_ij/fac_shield(j)
- enddo
- endif
+ enddo
+ endif
#ifdef NEWCORR
- geel_loc_ij=(a22*gmuij1(1)
+ geel_loc_ij=(a22*gmuij1(1)
& +a23*gmuij1(2)
& +a32*gmuij1(3)
& +a33*gmuij1(4))
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss1
c write(iout,*) "derivative over thatai"
c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3),
c & a33*gmuij1(4)
- gloc(nphi+i,icg)=gloc(nphi+i,icg)+
+ gloc(nphi+i,icg)=gloc(nphi+i,icg)+
& geel_loc_ij*wel_loc
c write(iout,*) "derivative over thatai-1"
c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3),
c & a33*gmuij2(4)
- geel_loc_ij=
+ geel_loc_ij=
& a22*gmuij2(1)
& +a23*gmuij2(2)
& +a32*gmuij2(3)
& +a33*gmuij2(4)
- gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
+ gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+
& geel_loc_ij*wel_loc
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss1
c Derivative over j residue
- geel_loc_ji=a22*gmuji1(1)
+ geel_loc_ji=a22*gmuji1(1)
& +a23*gmuji1(2)
& +a32*gmuji1(3)
& +a33*gmuji1(4)
gloc(nphi+j,icg)=gloc(nphi+j,icg)+
& geel_loc_ji*wel_loc
- & *fac_shield(i)*fac_shield(j)
+ & *fac_shield(i)*fac_shield(j)*sss1
- geel_loc_ji=
+ geel_loc_ji=
& +a22*gmuji2(1)
& +a23*gmuji2(2)
& +a32*gmuji2(3)
c write(iout,*) "derivative over thataj-1"
c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3),
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)
+ gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+
+ & geel_loc_ji*wel_loc
+ & *fac_shield(i)*fac_shield(j)*sss1
#endif
-cC Partial derivatives in virtual-bond dihedral angles gamma
- if (i.gt.1)
+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)
+ & (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
+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)
+
+ 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
+c & *fac_shield(i)*fac_shield(j)
+c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
C Derivatives of eello in DC(i+1) thru DC(j-1) or DC(nres-2)
- do l=1,3
- 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)
+ aux=eel_loc_ij/sss1*sssgrad1*rmij
+ ggg(1)=aux*xj
+ ggg(2)=aux*yj
+ ggg(3)=aux*zj
+ 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
+c & *fac_shield(i)*fac_shield(j)
+c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
- gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
- gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
-cgrad ghalf=0.5d0*ggg(l)
-cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
-cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
- enddo
+ gel_loc_long(l,j)=gel_loc_long(l,j)+ggg(l)
+ gel_loc_long(l,i)=gel_loc_long(l,i)-ggg(l)
+cgrad ghalf=0.5d0*ggg(l)
+cgrad gel_loc(l,i)=gel_loc(l,i)+ghalf
+cgrad gel_loc(l,j)=gel_loc(l,j)+ghalf
+ enddo
cgrad do k=i+1,j2
cgrad do l=1,3
cgrad gel_loc(l,k)=gel_loc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
C Remaining derivatives of eello
- 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)
+c gel_loc_long(3,j)=gel_loc_long(3,j)+ &
+c ssgradlipj*eel_loc_ij/2.0d0*lipscale/ &
+c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
+c
+c gel_loc_long(3,i)=gel_loc_long(3,i)+ &
+c ssgradlipi*eel_loc_ij/2.0d0*lipscale/ &
+c ((sslipi+sslipj)/2.0d0*lipscale+1.0d0)*sss_ele_cut
- 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)
+ 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
+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)
+ 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
+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)
+ 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
+c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
- enddo
- ENDIF
+ 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
+c & *((sslipi+sslipj)/2.0d0*lipscale+1.0d0)
+
+ enddo
+ ENDIF
+#ifdef FOURBODY
C Change 12/26/95 to calculate four-body contributions to H-bonding energy
c if (j.gt.i+1 .and. num_conti.le.maxconts) then
- if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
- & .and. num_conti.le.maxconts) then
+ if (wcorr+wcorr4+wcorr5+wcorr6.gt.0.0d0
+ & .and. num_conti.le.maxconts) then
c write (iout,*) i,j," entered corr"
C
C Calculate the contact function. The ith column of the array JCONT will
C contain the numbers of atoms that make contacts with the atom I (of numbers
C greater than I). The arrays FACONT and GACONT will contain the values of
C the contact function and its derivative.
-c r0ij=1.02D0*rpp(iteli,itelj)
-c r0ij=1.11D0*rpp(iteli,itelj)
- r0ij=2.20D0*rpp(iteli,itelj)
-c r0ij=1.55D0*rpp(iteli,itelj)
- call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
- if (fcont.gt.0.0D0) then
- num_conti=num_conti+1
- if (num_conti.gt.maxconts) then
- write (iout,*) 'WARNING - max. # of contacts exceeded;',
- & ' will skip next contacts for this conf.'
- else
- jcont_hb(num_conti,i)=j
-cd write (iout,*) "i",i," j",j," num_conti",num_conti,
-cd & " jcont_hb",jcont_hb(num_conti,i)
- IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
- & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
+c r0ij=1.02D0*rpp(iteli,itelj)
+c r0ij=1.11D0*rpp(iteli,itelj)
+ r0ij=2.20D0*rpp(iteli,itelj)
+c r0ij=1.55D0*rpp(iteli,itelj)
+ call gcont(rij,r0ij,1.0D0,0.2d0*r0ij,fcont,fprimcont)
+ if (fcont.gt.0.0D0) then
+ num_conti=num_conti+1
+ if (num_conti.gt.maxconts) then
+ write (iout,*) 'WARNING - max. # of contacts exceeded;',
+ & ' will skip next contacts for this conf.'
+ else
+ jcont_hb(num_conti,i)=j
+cd write (iout,*) "i",i," j",j," num_conti",num_conti,
+cd " jcont_hb",jcont_hb(num_conti,i)
+ IF (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.
+ & wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0) THEN
C 9/30/99 (AL) - store components necessary to evaluate higher-order loc-el
C terms.
- d_cont(num_conti,i)=rij
+ d_cont(num_conti,i)=rij
cd write (2,'(3e15.5)') rij,r0ij+0.2d0*r0ij,rij
C --- Electrostatic-interaction matrix ---
- a_chuj(1,1,num_conti,i)=a22
- a_chuj(1,2,num_conti,i)=a23
- a_chuj(2,1,num_conti,i)=a32
- a_chuj(2,2,num_conti,i)=a33
+ a_chuj(1,1,num_conti,i)=a22
+ a_chuj(1,2,num_conti,i)=a23
+ a_chuj(2,1,num_conti,i)=a32
+ a_chuj(2,2,num_conti,i)=a33
C --- Gradient of rij
- do kkk=1,3
- grij_hb_cont(kkk,num_conti,i)=erij(kkk)
- enddo
- kkll=0
- do k=1,2
- do l=1,2
- kkll=kkll+1
- do m=1,3
- a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
- a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
- a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
- a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
- a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
- enddo
+ do kkk=1,3
+ grij_hb_cont(kkk,num_conti,i)=erij(kkk)
+ enddo
+ kkll=0
+ do k=1,2
+ do l=1,2
+ kkll=kkll+1
+ do m=1,3
+ a_chuj_der(k,l,m,1,num_conti,i)=agg(m,kkll)
+ a_chuj_der(k,l,m,2,num_conti,i)=aggi(m,kkll)
+ a_chuj_der(k,l,m,3,num_conti,i)=aggi1(m,kkll)
+ a_chuj_der(k,l,m,4,num_conti,i)=aggj(m,kkll)
+ a_chuj_der(k,l,m,5,num_conti,i)=aggj1(m,kkll)
enddo
enddo
- ENDIF
- IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
+ enddo
+ ENDIF
+ IF (wcorr4.eq.0.0d0 .and. wcorr.gt.0.0d0) THEN
C Calculate contact energies
- cosa4=4.0D0*cosa
- wij=cosa-3.0D0*cosb*cosg
- cosbg1=cosb+cosg
- cosbg2=cosb-cosg
-c fac3=dsqrt(-ael6i)/r0ij**3
- fac3=dsqrt(-ael6i)*r3ij
-c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
- ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
- if (ees0tmp.gt.0) then
- ees0pij=dsqrt(ees0tmp)
- else
- ees0pij=0
- endif
-c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
- ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
- if (ees0tmp.gt.0) then
- ees0mij=dsqrt(ees0tmp)
- else
- ees0mij=0
- endif
-c ees0mij=0.0D0
- if (shield_mode.eq.0) then
+ cosa4=4.0D0*cosa
+ wij=cosa-3.0D0*cosb*cosg
+ cosbg1=cosb+cosg
+ cosbg2=cosb-cosg
+c fac3=dsqrt(-ael6i)/r0ij**3
+ fac3=dsqrt(-ael6i)*r3ij
+c ees0pij=dsqrt(4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1)
+ ees0tmp=4.0D0+cosa4+wij*wij-3.0D0*cosbg1*cosbg1
+ if (ees0tmp.gt.0) then
+ ees0pij=dsqrt(ees0tmp)
+ else
+ ees0pij=0
+ endif
+c ees0mij=dsqrt(4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2)
+ ees0tmp=4.0D0-cosa4+wij*wij-3.0D0*cosbg2*cosbg2
+ if (ees0tmp.gt.0) then
+ ees0mij=dsqrt(ees0tmp)
+ else
+ ees0mij=0
+ endif
+c ees0mij=0.0D0
+ if (shield_mode.eq.0) then
fac_shield(i)=1.0d0
fac_shield(j)=1.0d0
- else
+ else
ees0plist(num_conti,i)=j
C fac_shield(i)=0.4d0
C fac_shield(j)=0.6d0
- endif
- ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
- & *fac_shield(i)*fac_shield(j)
- ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
- & *fac_shield(i)*fac_shield(j)
+ endif
+ ees0p(num_conti,i)=0.5D0*fac3*(ees0pij+ees0mij)
+ & *fac_shield(i)*fac_shield(j)*sss1
+ ees0m(num_conti,i)=0.5D0*fac3*(ees0pij-ees0mij)
+ & *fac_shield(i)*fac_shield(j)*sss1
C Diagnostics. Comment out or remove after debugging!
c ees0p(num_conti,i)=0.5D0*fac3*ees0pij
c write (iout,*) 'i=',i,' j=',j,' rij=',rij,' r0ij=',r0ij,
c & ' ees0ij=',ees0p(num_conti,i),ees0m(num_conti,i),' fcont=',fcont
C Angular derivatives of the contact function
- ees0pij1=fac3/ees0pij
- ees0mij1=fac3/ees0mij
- fac3p=-3.0D0*fac3*rrmij
- ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
- ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
-c ees0mij1=0.0D0
- ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
- ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
- ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
- ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
- ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
- ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
- ecosap=ecosa1+ecosa2
- ecosbp=ecosb1+ecosb2
- ecosgp=ecosg1+ecosg2
- ecosam=ecosa1-ecosa2
- ecosbm=ecosb1-ecosb2
- ecosgm=ecosg1-ecosg2
+ ees0pij1=fac3/ees0pij
+ ees0mij1=fac3/ees0mij
+ fac3p=-3.0D0*fac3*rrmij
+ ees0pijp=0.5D0*fac3p*(ees0pij+ees0mij)
+ ees0mijp=0.5D0*fac3p*(ees0pij-ees0mij)
+c ees0mij1=0.0D0
+ ecosa1= ees0pij1*( 1.0D0+0.5D0*wij)
+ ecosb1=-1.5D0*ees0pij1*(wij*cosg+cosbg1)
+ ecosg1=-1.5D0*ees0pij1*(wij*cosb+cosbg1)
+ ecosa2= ees0mij1*(-1.0D0+0.5D0*wij)
+ ecosb2=-1.5D0*ees0mij1*(wij*cosg+cosbg2)
+ ecosg2=-1.5D0*ees0mij1*(wij*cosb-cosbg2)
+ ecosap=ecosa1+ecosa2
+ ecosbp=ecosb1+ecosb2
+ ecosgp=ecosg1+ecosg2
+ ecosam=ecosa1-ecosa2
+ ecosbm=ecosb1-ecosb2
+ ecosgm=ecosg1-ecosg2
C Diagnostics
c ecosap=ecosa1
c ecosbp=ecosb1
c ecosbm=0.0D0
c ecosgm=0.0D0
C End diagnostics
- facont_hb(num_conti,i)=fcont
- fprimcont=fprimcont/rij
+ facont_hb(num_conti,i)=fcont
+ fprimcont=fprimcont/rij
cd facont_hb(num_conti,i)=1.0D0
C Following line is for diagnostics.
cd fprimcont=0.0D0
- do k=1,3
- dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
- dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
- enddo
- do k=1,3
- gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
- gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
- enddo
- gggp(1)=gggp(1)+ees0pijp*xj
- gggp(2)=gggp(2)+ees0pijp*yj
- gggp(3)=gggp(3)+ees0pijp*zj
- gggm(1)=gggm(1)+ees0mijp*xj
- gggm(2)=gggm(2)+ees0mijp*yj
- gggm(3)=gggm(3)+ees0mijp*zj
+ do k=1,3
+ dcosb(k)=rmij*(dc_norm(k,i)-erij(k)*cosb)
+ dcosg(k)=rmij*(dc_norm(k,j)-erij(k)*cosg)
+ enddo
+ do k=1,3
+ gggp(k)=ecosbp*dcosb(k)+ecosgp*dcosg(k)
+ gggm(k)=ecosbm*dcosb(k)+ecosgm*dcosg(k)
+ enddo
+ gggp(1)=gggp(1)+ees0pijp*xj
+ & +ees0p(num_conti,i)/sss1*rmij*xj*sssgrad1
+ gggp(2)=gggp(2)+ees0pijp*yj
+ & +ees0p(num_conti,i)/sss1*rmij*yj*sssgrad1
+ gggp(3)=gggp(3)+ees0pijp*zj
+ & +ees0p(num_conti,i)/sss1*rmij*zj*sssgrad1
+ gggm(1)=gggm(1)+ees0mijp*xj
+ & +ees0m(num_conti,i)/sss1*rmij*xj*sssgrad1
+ gggm(2)=gggm(2)+ees0mijp*yj
+ & +ees0m(num_conti,i)/sss1*rmij*yj*sssgrad1
+ gggm(3)=gggm(3)+ees0mijp*zj
+ & +ees0m(num_conti,i)/sss1*rmij*zj*sssgrad1
C Derivatives due to the contact function
- gacont_hbr(1,num_conti,i)=fprimcont*xj
- gacont_hbr(2,num_conti,i)=fprimcont*yj
- gacont_hbr(3,num_conti,i)=fprimcont*zj
- do k=1,3
+ gacont_hbr(1,num_conti,i)=fprimcont*xj
+ gacont_hbr(2,num_conti,i)=fprimcont*yj
+ gacont_hbr(3,num_conti,i)=fprimcont*zj
+ do k=1,3
c
c 10/24/08 cgrad and ! comments indicate the parts of the code removed
c following the change of gradient-summation algorithm.
c
cgrad ghalfp=0.5D0*gggp(k)
cgrad ghalfm=0.5D0*gggm(k)
- gacontp_hb1(k,num_conti,i)=!ghalfp
- & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
- & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
- & *fac_shield(i)*fac_shield(j)
+ gacontp_hb1(k,num_conti,i)=!ghalfp
+ & +(ecosap*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosbp*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ & *sss1*fac_shield(i)*fac_shield(j)
- gacontp_hb2(k,num_conti,i)=!ghalfp
- & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
- & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- & *fac_shield(i)*fac_shield(j)
+ gacontp_hb2(k,num_conti,i)=!ghalfp
+ & +(ecosap*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosgp*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ & *sss1*fac_shield(i)*fac_shield(j)
- gacontp_hb3(k,num_conti,i)=gggp(k)
- & *fac_shield(i)*fac_shield(j)
+ gacontp_hb3(k,num_conti,i)=gggp(k)
+ & *sss1*fac_shield(i)*fac_shield(j)
- gacontm_hb1(k,num_conti,i)=!ghalfm
- & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
- & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
- & *fac_shield(i)*fac_shield(j)
+ gacontm_hb1(k,num_conti,i)=!ghalfm
+ & +(ecosam*(dc_norm(k,j)-cosa*dc_norm(k,i))
+ & + ecosbm*(erij(k)-cosb*dc_norm(k,i)))*vbld_inv(i+1)
+ & *sss1*fac_shield(i)*fac_shield(j)
- gacontm_hb2(k,num_conti,i)=!ghalfm
- & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
- & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
- & *fac_shield(i)*fac_shield(j)
+ gacontm_hb2(k,num_conti,i)=!ghalfm
+ & +(ecosam*(dc_norm(k,i)-cosa*dc_norm(k,j))
+ & + ecosgm*(erij(k)-cosg*dc_norm(k,j)))*vbld_inv(j+1)
+ & *sss1*fac_shield(i)*fac_shield(j)
- gacontm_hb3(k,num_conti,i)=gggm(k)
- & *fac_shield(i)*fac_shield(j)
+ gacontm_hb3(k,num_conti,i)=gggm(k)
+ & *sss1*fac_shield(i)*fac_shield(j)
- enddo
- ENDIF ! wcorr
- endif ! num_conti.le.maxconts
- endif ! fcont.gt.0
- endif ! j.gt.i+1
- if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
- do k=1,4
- do l=1,3
- ghalf=0.5d0*agg(l,k)
- aggi(l,k)=aggi(l,k)+ghalf
- aggi1(l,k)=aggi1(l,k)+agg(l,k)
- aggj(l,k)=aggj(l,k)+ghalf
enddo
+ ENDIF ! wcorr
+ endif ! num_conti.le.maxconts
+ endif ! fcont.gt.0
+ endif ! j.gt.i+1
+#endif
+ if (wturn3.gt.0.0d0 .or. wturn4.gt.0.0d0) then
+ do k=1,4
+ do l=1,3
+ ghalf=0.5d0*agg(l,k)
+ aggi(l,k)=aggi(l,k)+ghalf
+ aggi1(l,k)=aggi1(l,k)+agg(l,k)
+ aggj(l,k)=aggj(l,k)+ghalf
+ enddo
+ enddo
+ if (j.eq.nres-1 .and. i.lt.j-2) then
+ do k=1,4
+ do l=1,3
+ aggj1(l,k)=aggj1(l,k)+agg(l,k)
enddo
- if (j.eq.nres-1 .and. i.lt.j-2) then
- do k=1,4
- do l=1,3
- aggj1(l,k)=aggj1(l,k)+agg(l,k)
- enddo
- enddo
- endif
- endif
+ enddo
+ endif
+ endif
c t_eelecij=t_eelecij+MPI_Wtime()-time00
return
end
C
C Compute Evdwpp
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.CONTROL'
include 'COMMON.IOUNITS'
include 'COMMON.CHAIN'
include 'COMMON.DERIV'
include 'COMMON.INTERACT'
- include 'COMMON.CONTACTS'
+c include 'COMMON.CONTACTS'
include 'COMMON.TORSION'
include 'COMMON.VECTORS'
include 'COMMON.FFIELD'
- dimension ggg(3)
+ include "COMMON.SPLITELE"
+ double precision ggg(3)
integer xshift,yshift,zshift
c 4/26/02 - AL scaling factor for 1,4 repulsive VDW interactions
#ifdef MOMENT
double precision scal_el /0.5d0/
#endif
c write (iout,*) "evdwpp_short"
+ integer i,j,k,iteli,itelj,num_conti,ind,isubchap
+ double precision dxi,dyi,dzi,dxj,dyj,dzj,aaa,bbb
+ double precision xj,yj,zj,rij,rrmij,r3ij,r6ij,evdw1,
+ & dx_normi,dy_normi,dz_normi,xmedi,ymedi,zmedi,
+ & dx_normj,dy_normj,dz_normj,rmij,ev1,ev2,evdwij,facvdw
+ double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
+ & dist_temp, dist_init,sss_grad
+ double precision sscale,sscagrad
evdw1=0.0D0
C print *,"WCHODZE"
c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
- xmedi=mod(xmedi,boxxsize)
- if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
- ymedi=mod(ymedi,boxysize)
- if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
- zmedi=mod(zmedi,boxzsize)
- if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
+ xmedi=mod(xmedi,boxxsize)
+ if (xmedi.lt.0.0d0) xmedi=xmedi+boxxsize
+ ymedi=mod(ymedi,boxysize)
+ if (ymedi.lt.0.0d0) ymedi=ymedi+boxysize
+ zmedi=mod(zmedi,boxzsize)
+ if (zmedi.lt.0.0d0) zmedi=zmedi+boxzsize
num_conti=0
c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
c & ' ielend',ielend_vdw(i)
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- isubchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
- xj=xj_safe+xshift*boxxsize
- yj=yj_safe+yshift*boxysize
- zj=zj_safe+zshift*boxzsize
- dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
- if(dist_temp.lt.dist_init) then
- dist_init=dist_temp
- xj_temp=xj
- yj_temp=yj
- zj_temp=zj
- isubchap=1
- endif
- enddo
- enddo
- enddo
- if (isubchap.eq.1) then
- xj=xj_temp-xmedi
- yj=yj_temp-ymedi
- zj=zj_temp-zmedi
- else
- xj=xj_safe-xmedi
- yj=yj_safe-ymedi
- zj=zj_safe-zmedi
- endif
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ isubchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ isubchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (isubchap.eq.1) then
+ xj=xj_temp-xmedi
+ yj=yj_temp-ymedi
+ zj=zj_temp-zmedi
+ else
+ xj=xj_safe-xmedi
+ yj=yj_safe-ymedi
+ zj=zj_safe-zmedi
+ endif
rij=xj*xj+yj*yj+zj*zj
rrmij=1.0D0/rij
rij=dsqrt(rij)
c sss=sscale(rij/rpp(iteli,itelj))
c sssgrad=sscagrad(rij/rpp(iteli,itelj))
- sss=sscale(rij)
- sssgrad=sscagrad(rij)
+ sss=sscale(rij/rpp(iteli,itelj),r_cut_respa)
+ sssgrad=sscagrad(rij/rpp(iteli,itelj),r_cut_respa)
if (sss.gt.0.0d0) then
rmij=1.0D0/rij
r3ij=rrmij*rmij
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)
+ 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)
C ggg(1)=facvdw*xj
C ggg(2)=facvdw*yj
C ggg(3)=facvdw*zj
include 'COMMON.FFIELD'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
+ include "COMMON.SPLITELE"
logical lprint_short
common /shortcheck/ lprint_short
- dimension ggg(3)
+ double precision ggg(3)
integer xshift,yshift,zshift
+ integer i,iint,j,k,iteli,itypj,subchap
+ double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
+ & fac,e1,e2,rij
+ double precision evdw2,evdw2_14,evdwij
+ double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
+ & dist_temp, dist_init
+ double precision sscale,sscagrad
if (energy_dec) write (iout,*) "escp_long:",r_cut,rlamb
evdw2=0.0D0
evdw2_14=0.0d0
xi=0.5D0*(c(1,i)+c(1,i+1))
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
if (itypj.eq.ntyp1) cycle
C Uncomment following three lines for SC-p interactions
c xj=c(1,nres+j)-xi
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
c end correction
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
xj=xj_safe+xshift*boxxsize
yj=yj_safe+yshift*boxysize
zj=zj_safe+zshift*boxzsize
zj_temp=zj
subchap=1
endif
- enddo
- enddo
- enddo
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
- sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
- sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
+ sss1=sscale(1.0d0/(dsqrt(rrij)),r_cut_int)
+ if (sss1.eq.0) cycle
+ sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
+ sssgrad=
+ & sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
+ sssgrad1=sscagrad(1.0d0/dsqrt(rrij),r_cut_int)
if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
& " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
if (sss.lt.1.0d0) then
if (iabs(j-i) .le. 2) then
e1=scal14*e1
e2=scal14*e2
- evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)
+ evdw2_14=evdw2_14+(e1+e2)*(1.0d0-sss)*sss1
endif
evdwij=e1+e2
- evdw2=evdw2+evdwij*(1.0d0-sss)
+ evdw2=evdw2+evdwij*(1.0d0-sss)*sss1
if (energy_dec) write (iout,'(a6,2i5,2(0pf7.3))')
& 'evdw2',i,j,sss,evdwij
C
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
C
- fac=-(evdwij+e1)*rrij*(1.0d0-sss)
- fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
+ fac=-(evdwij+e1)*rrij*(1.0d0-sss)*sss1
+ fac=fac+evdwij*dsqrt(rrij)*(-sssgrad/rscp(itypj,iteli)
+ & +sssgrad1)/expon
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
C peptide-group centers and side chains and its gradient in virtual-bond and
C side-chain vectors.
C
- implicit real*8 (a-h,o-z)
+ implicit none
include 'DIMENSIONS'
include 'COMMON.GEO'
include 'COMMON.VAR'
include 'COMMON.FFIELD'
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
+ include "COMMON.SPLITELE"
integer xshift,yshift,zshift
logical lprint_short
common /shortcheck/ lprint_short
- dimension ggg(3)
+ integer i,iint,j,k,iteli,itypj,subchap
+ double precision xi,yi,zi,xj,yj,zj,rrij,sss1,sssgrad1,
+ & fac,e1,e2,rij
+ double precision evdw2,evdw2_14,evdwij
+ double precision xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,
+ & dist_temp, dist_init
+ double precision ggg(3)
+ double precision sscale,sscagrad
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,rlamb
+ if (energy_dec) write (iout,*) "escp_short:",r_cut_int,rlamb
do i=iatscp_s,iatscp_e
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))
yi=0.5D0*(c(2,i)+c(2,i+1))
zi=0.5D0*(c(3,i)+c(3,i+1))
- xi=mod(xi,boxxsize)
- if (xi.lt.0) xi=xi+boxxsize
- yi=mod(yi,boxysize)
- if (yi.lt.0) yi=yi+boxysize
- zi=mod(zi,boxzsize)
- if (zi.lt.0) zi=zi+boxzsize
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
c if (lprint_short)
c & write (iout,*) "i",i," itype",itype(i),itype(i+1),
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
- itypj=itype(j)
+ itypj=iabs(itype(j))
c if (lprint_short)
c & write (iout,*) "j",j," itypj",itypj
if (itypj.eq.ntyp1) cycle
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
c end correction
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
c if (lprint_short) then
c write (iout,*) i,j,xi,yi,zi,xj,yj,zj
c write (iout,*) "dist_init",dsqrt(dist_init)
c endif
- xj_safe=xj
- yj_safe=yj
- zj_safe=zj
- subchap=0
- do xshift=-1,1
- do yshift=-1,1
- do zshift=-1,1
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
xj=xj_safe+xshift*boxxsize
yj=yj_safe+yshift*boxysize
zj=zj_safe+zshift*boxzsize
zj_temp=zj
subchap=1
endif
- enddo
- enddo
- enddo
+ enddo
+ enddo
+ enddo
c if (lprint_short) write (iout,*) "dist_temp",dsqrt(dist_temp)
- if (subchap.eq.1) then
- xj=xj_temp-xi
- yj=yj_temp-yi
- zj=zj_temp-zi
- else
- xj=xj_safe-xi
- yj=yj_safe-yi
- zj=zj_safe-zi
- endif
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
c sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
c sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)))
- sss=sscale(1.0d0/(dsqrt(rrij)))
- sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
+ sss=sscale(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),r_cut_respa)
+ sssgrad=sscagrad(1.0d0/(dsqrt(rrij)*rscp(itypj,iteli)),
+ & r_cut_respa)
if (energy_dec) write (iout,*) "rrij",1.0d0/dsqrt(rrij),
& " rscp",rscp(itypj,iteli)," subchap",subchap," sss",sss
c if (lprint_short) write (iout,*) "rij",1.0/dsqrt(rrij),
C Calculate contributions to the gradient in the virtual-bond and SC vectors.
C
fac=-(evdwij+e1)*rrij*sss
- fac=fac+(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
+ fac=fac+evdwij*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)/expon
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
projects / unres.git / blobdiff