xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=0
call eelecij(i,i+2,ees,evdw1,eel_loc)
if (wturn3.gt.0.0d0) call eturn3(i,eello_turn3)
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
num_conti=num_cont_hb(i)
call eelecij(i,i+3,ees,evdw1,eel_loc)
if (wturn4.gt.0.0d0 .and. itype(i+2).ne.ntyp1) &
xmedi=c(1,i)+0.5d0*dxi
ymedi=c(2,i)+0.5d0*dyi
zmedi=c(3,i)+0.5d0*dzi
+ xmedi=dmod(xmedi,boxxsize)
+ if (xmedi.lt.0) xmedi=xmedi+boxxsize
+ ymedi=dmod(ymedi,boxysize)
+ if (ymedi.lt.0) ymedi=ymedi+boxysize
+ zmedi=dmod(zmedi,boxzsize)
+ if (zmedi.lt.0) zmedi=zmedi+boxzsize
+
! write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i)
num_conti=num_cont_hb(i)
do j=ielstart(i),ielend(i)
real(kind=8),dimension(2,2) :: acipa !el,a_temp
!el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1
real(kind=8),dimension(4) :: muij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+ integer xshift,yshift,zshift
!el integer :: num_conti,j1,j2
!el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,&
!el dz_normi,xmedi,ymedi,zmedi
0.0d0,0.0d0,1.0d0/),shape(unmat))
! integer :: maxconts=nres/4
!el local variables
- integer :: k,i,j,iteli,itelj,kkk,l,kkll,m
+ integer :: k,i,j,iteli,itelj,kkk,l,kkll,m,isubchap
real(kind=8) :: ael6i,rrmij,rmij,r0ij,fcont,fprimcont,ees0tmp
real(kind=8) :: ees,evdw1,eel_loc,aaa,bbb,ael3i
real(kind=8) :: dxj,dyj,dzj,dx_normj,dy_normj,dz_normj,xj,yj,zj,&
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-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)
+ 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
+!C print *,i,j
+ 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_ele_cut=sscale_ele(rij)
+ sss_ele_grad=sscagrad_ele(rij)
+! print *,sss_ele_cut,sss_ele_grad,&
+! 1.0d0/(rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) go to 128
+
rmij=1.0D0/rij
r3ij=rrmij*rmij
r6ij=r3ij*r3ij
eesij=el1+el2
! 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
+ ees=ees+eesij*sss_ele_cut
+ evdw1=evdw1+evdwij*sss_ele_cut
!d write(iout,'(2(2i3,2x),7(1pd12.4)/2(3(1pd12.4),5x)/)')
!d & iteli,i,itelj,j,aaa,bbb,ael6i,ael3i,
!d & 1.0D0/dsqrt(rrmij),evdwij,eesij,
! Calculate contributions to the Cartesian gradient.
!
#ifdef SPLITELE
- facvdw=-6*rrmij*(ev1+evdwij)
- facel=-3*rrmij*(el1+eesij)
+ facvdw=-6*rrmij*(ev1+evdwij)*sss_ele_cut
+ facel=-3*rrmij*(el1+eesij)*sss_ele_cut
fac1=fac
erij(1)=xj*rmij
erij(2)=yj*rmij
!grad gelc(l,k)=gelc(l,k)+ggg(l)
!grad enddo
!grad enddo
- ggg(1)=facvdw*xj
- ggg(2)=facvdw*yj
- ggg(3)=facvdw*zj
+ ggg(1)=facvdw*xj+sss_ele_grad*rmij*evdwij*xj
+ ggg(2)=facvdw*yj+sss_ele_grad*rmij*evdwij*yj
+ ggg(3)=facvdw*zj+sss_ele_grad*rmij*evdwij*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gvdwpp(k,i)=gvdwpp(k,i)+ghalf
!grad enddo
!grad enddo
#else
- facvdw=ev1+evdwij
- facel=el1+eesij
+ facvdw=(ev1+evdwij)*sss_ele_cut
+ facel=(el1+eesij)*sss_ele_cut
fac1=fac
fac=-3*rrmij*(facvdw+facvdw+facel)
erij(1)=xj*rmij
!
! Radial derivatives. First process both termini of the fragment (i,j)
!
- ggg(1)=fac*xj
- ggg(2)=fac*yj
- ggg(3)=fac*zj
+ ggg(1)=fac*xj+sss_ele_grad*rmij*(eesij+evdwij)*xj
+ ggg(2)=fac*yj+sss_ele_grad*rmij*(eesij+evdwij)*yj
+ ggg(3)=fac*zj+sss_ele_grad*rmij*(eesij+evdwij)*zj
! do k=1,3
! ghalf=0.5D0*ggg(k)
! gelc(k,i)=gelc(k,i)+ghalf
gelc_long(k,j)=gelc_long(k,j)+ggg(k)
gelc_long(k,i)=gelc_long(k,i)-ggg(k)
enddo
+ 128 continue
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
endif
return
end function sscale
+ real(kind=8) function sscale_grad(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut-rlamb) then
+ sscale_grad=0.0d0
+ else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then
+ gamm=(r-(r_cut-rlamb))/rlamb
+ sscale_grad=gamm*(6*gamm-6.0d0)/rlamb
+ else
+ sscale_grad=0d0
+ endif
+ return
+ end function sscale_grad
+
!!!!!!!!!! PBCSCALE
real(kind=8) function sscale_ele(r)
! include "COMMON.SPLITELE"
!el local variables
integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig,sig0ij,rij_shift
- real(kind=8) :: sss,e1,e2,evdw
+ real(kind=8) :: sss,e1,e2,evdw,sss_grad
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
dist_temp, dist_init
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
if (sss_ele_cut.le.0.0) cycle
if (sss.lt.1.0d0) then
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
- fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij-sss_grad/(1.0-sss)*rij &
+ /sigmaii(itypi,itypj))
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
!el local variables
integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,fac,sigm,epsi,sig0ij,sig
- real(kind=8) :: sss,e1,e2,evdw,rij_shift
+ real(kind=8) :: sss,e1,e2,evdw,rij_shift,sss_grad
real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
dist_temp, dist_init
evdw=0.0D0
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
sss=sscale(1.0d0/(rij*sigmaii(itypi,itypj)))
+ sss_grad=sscale_grad(1.0d0/(rij*sigmaii(itypi,itypj)))
sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
if (sss_ele_cut.le.0.0) cycle
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*sss
+ evdw=evdw+evdwij*sss*sss_ele_cut
if (lprn) then
sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
- fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
- /sigma(itypi,itypj)*rij
+ fac=fac+evdwij*(sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij+sss_grad/sss*rij &
+ /sigmaii(itypi,itypj))
! fac=0.0d0
! Calculate the radial part of the gradient
projects / unres4.git / blobdiff