if (yi.lt.0) yi=yi+boxysize
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
if (zj.lt.buflipbot) then
C what fraction I am in
fracinbuf=1.0d0-
- & ((positi-bordlipbot)/lipbufthick)
+ & ((zj-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)
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
sslipj=sscalelip(fracinbuf)
ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
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
+ & +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
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
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)=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 Calculate angular part of the gradient.
call sc_grad_scale(1.0d0-sss)
endif
include 'COMMON.CALC'
include 'COMMON.CONTROL'
logical lprn
+ integer xshift,yshift,zshift
evdw=0.0D0
ccccc energy_dec=.false.
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
if (yi.lt.0) yi=yi+boxysize
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
if (zj.lt.buflipbot) then
C what fraction I am in
fracinbuf=1.0d0-
- & ((positi-bordlipbot)/lipbufthick)
+ & ((zj-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)
+ elseif (zj.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
sslipj=sscalelip(fracinbuf)
ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
else
ssgradlipj=0.0
endif
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
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
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)=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 Calculate angular part of the gradient.
call sc_grad_scale(sss)
endif
enddo ! j
enddo ! iint
enddo ! i
-c write (iout,*) "Number of loop steps in EGB:",ind
+C write (iout,*) "Number of loop steps in EGB:",ind
cccc energy_dec=.false.
return
end
enddo
do k=1,3
gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
+ gg_lipi(k)=gg_lipi(k)*scalfac
+ gg_lipj(k)=gg_lipj(k)*scalfac
enddo
c write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
double precision unmat(3,3) /1.0d0,0.0d0,0.0d0,
& 0.0d0,1.0d0,0.0d0,
& 0.0d0,0.0d0,1.0d0/
+ integer xhift,yshift,zshift
c time00=MPI_Wtime()
cd write (iout,*) "eelecij",i,j
C print *,"WCHODZE2"
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))
+ sss=sscale(rij)
+ sssgrad=sscagrad(rij)
r3ij=rrmij*rmij
r6ij=r3ij*r3ij
cosa=dx_normi*dx_normj+dy_normi*dy_normj+dz_normi*dz_normj
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)
+ ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj
+ ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj
+ ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
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)
+ ggg(1)=facvdw*xj-sssgrad*rmij*evdwij*xj
+ ggg(2)=facvdw*yj-sssgrad*rmij*evdwij*yj
+ ggg(3)=facvdw*zj-sssgrad*rmij*evdwij*zj
do k=1,3
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
rij=xj*xj+yj*yj+zj*zj
rrmij=1.0D0/rij
rij=dsqrt(rij)
- sss=sscale(rij/rpp(iteli,itelj))
- sssgrad=sscagrad(rij/rpp(iteli,itelj))
+ sss=sscale(rij)
+ sssgrad=sscagrad(rij)
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
+ ggg(2)=facvdw*yj+sssgrad*rmij*evdwij*yj
+ ggg(3)=facvdw*zj+sssgrad*rmij*evdwij*zj
C ggg(1)=facvdw*xj
C ggg(2)=facvdw*yj
C ggg(3)=facvdw*zj
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
dimension ggg(3)
+ integer xshift,yshift,zshift
evdw2=0.0D0
evdw2_14=0.0d0
CD print '(a)','Enter ESCP KURWA'
xj=c(1,j)
yj=c(2,j)
zj=c(3,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(zi,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
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)))
+ sss=sscale(1.0d0/(dsqrt(rrij)))
+ sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
if (sss.lt.1.0d0) then
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
C
fac=-(evdwij+e1)*rrij*(1.0d0-sss)
- fac=fac-(evdwij)*sssgrad*dsqrt(rrij)/rscp(itypj,iteli)
+ fac=fac-(evdwij)*sssgrad*dsqrt(rrij)
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
include 'COMMON.IOUNITS'
include 'COMMON.CONTROL'
dimension ggg(3)
+ integer xshift,yshift,zshift
evdw2=0.0D0
evdw2_14=0.0d0
cd print '(a)','Enter ESCP'
xj=c(1,j)
yj=c(2,j)
zj=c(3,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(zi,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
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)))
+ sss=sscale(1.0d0/(dsqrt(rrij)))
+ sssgrad=sscagrad(1.0d0/(dsqrt(rrij)))
if (sss.gt.0.0d0) then
fac=rrij**expon2
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)
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac