C-----------------------------------------------------------------------
+ double precision function sscalelip(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+C if(r.lt.r_cut-rlamb) then
+C sscale=1.0d0
+C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+C gamm=(r-(r_cut-rlamb))/rlamb
+ sscalelip=1.0d0+r*r*(2*r-3.0d0)
+C else
+C sscale=0d0
+C endif
+ return
+ end
+C-----------------------------------------------------------------------
+ double precision function sscagradlip(r)
+ double precision r,gamm
+ include "COMMON.SPLITELE"
+C if(r.lt.r_cut-rlamb) then
+C sscagrad=0.0d0
+C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+C gamm=(r-(r_cut-rlamb))/rlamb
+ sscagradlip=r*(6*r-6.0d0)
+C else
+C sscagrad=0.0d0
+C endif
+ return
+ end
+
+C-----------------------------------------------------------------------
double precision function sscale(r)
double precision r,gamm
include "COMMON.SPLITELE"
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
evdw=evdw+(1.0d0-sss)*evdwij
C
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
evdw=evdw+sss*evdwij
C
if (sss.lt.1.0d0) then
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
if (sss.gt.0.0d0) then
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*(1.0d0-sss)
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*sss
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-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
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*(1.0d0-sss)
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,chi1,chi2,chip1,chip2,
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*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)
endif
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-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
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*sss
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,chi1,chi2,chip1,chip2,
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=ssgradlipi*evdwij
+ gg_lipj(3)=ssgradlipj*evdwij
C Calculate angular part of the gradient.
call sc_grad_scale(sss)
endif
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+(evdwij+e_augm)*sss
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
enddo
c write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
& +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
& +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
C Calculate the components of the gradient in DC and X
C
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
enddo
return
end
cd enddo
cd call check_vecgrad
cd stop
+C print *,"WCHODZE3"
if (icheckgrad.eq.1) then
do i=1,nres-1
fac=1.0d0/dsqrt(scalar(dc(1,i),dc(1,i)))
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) xmedi=xmedi+boxxsize
+ ymedi=mod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=mod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=num_cont_hb(i)
call eelecij_scale(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1)
& 0.0d0,0.0d0,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 (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
xj_safe=xj
yj_safe=yj
zj_safe=zj
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
+ 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
double precision scal_el /0.5d0/
#endif
evdw1=0.0D0
+C print *,"WCHODZE"
c write (iout,*) "iatel_s_vdw",iatel_s_vdw,
c & " iatel_e_vdw",iatel_e_vdw
call flush(iout)
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
- dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ dist_init=(xj-xmedi)**2+(yj-ymedi)**2+(zj-zmedi)**2
xj_safe=xj
yj_safe=yj
zj_safe=zj
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
+ 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