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,
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ 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
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ 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
+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
+ 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
C Calculate angle-dependent terms of energy and contributions to their
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,
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij/(1.0-sss)*(-sssgrad)/sigmaii(itypi,itypj)*rij
c fac=0.0d0
C Calculate the radial part of the gradient
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
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ 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
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
alf1=alp(itypi)
alf2=alp(itypj)
alf12=0.5D0*(alf1+alf2)
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+ 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
+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
+ 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.gt.0.0d0) then
C Calculate angle-dependent terms of energy and contributions to their
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,
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij/sss*sssgrad/sigmaii(itypi,itypj)*rij
c fac=0.0d0
C Calculate the radial part of the gradient
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)))
C
do i=iturn3_start,iturn3_end
if (itype(i).eq.ntyp1.or. itype(i+1).eq.ntyp1
- & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1) cycle
+ & .or. itype(i+2).eq.ntyp1 .or. itype(i+3).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
+ & .or. itype(i+4).eq.ntyp1
+ & ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,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=0
call eelecij_scale(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
do i=iturn4_start,iturn4_end
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
& .or. itype(i+3).eq.ntyp1
- & .or. itype(i+4).eq.ntyp1) cycle
+ & .or. itype(i+4).eq.ntyp1
+ & .or. itype(i+5).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
+ & ) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,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)
c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3
c
do i=iatel_s,iatel_e
- if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
+ if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1
+ & .or. itype(i+2).eq.ntyp1
+ & .or. itype(i-1).eq.ntyp1
+ &) cycle
dxi=dc(1,i)
dyi=dc(2,i)
dzi=dc(3,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
c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
- if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1) cycle
+ if (itype(j).eq.ntyp1 .or. itype(j+1).eq.ntyp1
+ & .or.itype(j+2).eq.ntyp1
+ & .or.itype(j-1).eq.ntyp1
+ &) cycle
call eelecij_scale(i,j,ees,evdw1,eel_loc)
enddo ! j
num_cont_hb(i)=num_conti
& 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)
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-xmedi
- yj=c(2,j)+0.5D0*dyj-ymedi
- zj=c(3,j)+0.5D0*dzj-zmedi
+ 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_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)
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
cgrad gelc(l,k)=gelc(l,k)+ggg(l)
cgrad enddo
cgrad enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- 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)
c do k=1,3
c ghalf=0.5D0*ggg(k)
c gvdwpp(k,i)=gvdwpp(k,i)+ghalf
cgrad enddo
cgrad enddo
c 9/28/08 AL Gradient compotents will be summed only at the end
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+C ggg(1)=facvdw*xj
+C ggg(2)=facvdw*yj
+C ggg(3)=facvdw*zj
+ 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)
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)
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=0
c write (iout,*) 'i',i,' ielstart',ielstart_vdw(i),
c & ' ielend',ielend_vdw(i)
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-xmedi
- yj=c(2,j)+0.5D0*dyj-ymedi
- zj=c(3,j)+0.5D0*dzj-zmedi
+ 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_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)
sss=sscale(rij/rpp(iteli,itelj))
+ sssgrad=sscagrad(rij/rpp(iteli,itelj))
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
- ggg(2)=facvdw*yj
- 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)
+C ggg(1)=facvdw*xj
+C ggg(2)=facvdw*yj
+C ggg(3)=facvdw*zj
do k=1,3
gvdwpp(k,j)=gvdwpp(k,j)+ggg(k)
gvdwpp(k,i)=gvdwpp(k,i)-ggg(k)
dimension ggg(3)
evdw2=0.0D0
evdw2_14=0.0d0
-cd print '(a)','Enter ESCP'
+CD print '(a)','Enter ESCP KURWA'
cd write (iout,*) 'iatscp_s=',iatscp_s,' iatscp_e=',iatscp_e
do i=iatscp_s,iatscp_e
if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) cycle
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
do iint=1,nscp_gr(i)
do j=iscpstart(i,iint),iscpend(i,iint)
c yj=c(2,nres+j)-yi
c zj=c(3,nres+j)-zi
C Uncomment following three lines for Ca-p interactions
- xj=c(1,j)-xi
- yj=c(2,j)-yi
- zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ 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
+
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)))
if (sss.lt.1.0d0) then
-
fac=rrij**expon2
e1=fac*fac*aad(itypj,iteli)
e2=fac*bad(itypj,iteli)
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)
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
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
do iint=1,nscp_gr(i)
c yj=c(2,nres+j)-yi
c zj=c(3,nres+j)-zi
C Uncomment following three lines for Ca-p interactions
- xj=c(1,j)-xi
- yj=c(2,j)-yi
- zj=c(3,j)-zi
+ xj=c(1,j)
+ yj=c(2,j)
+ zj=c(3,j)
+ 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
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)))
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)
ggg(1)=xj*fac
ggg(2)=yj*fac
ggg(3)=zj*fac
projects / unres.git / blobdiff