#ifdef MPI
real(kind=8) :: weights_(n_ene) !,time_Bcast,time_Bcastw
+! shielding effect varibles for MPI
+! real(kind=8) fac_shieldbuf(maxres),
+! & grad_shield_locbuf(3,maxcontsshi,-1:maxres),
+! & grad_shield_sidebuf(3,maxcontsshi,-1:maxres),
+! & grad_shieldbuf(3,-1:maxres)
+! integer ishield_listbuf(maxres),
+! &shield_listbuf(maxcontsshi,maxres)
+
! print*,"ETOTAL Processor",fg_rank," absolute rank",myrank,
! & " nfgtasks",nfgtasks
if (nfgtasks.gt.1) then
! include 'COMMON.SBRIDGE'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ integer :: iint,itypi,itypi1,itypj,subchap
real(kind=8) :: rrij,xi,yi,zi,sig,rij_shift,fac,e1,e2,sigm,epsi
real(kind=8) :: evdw,sig0ij
+ real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
integer :: ii
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ xi=dmod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=dmod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=dmod(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=0.0D0
! alf2=0.0D0
! alf12=0.0D0
- 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=dmod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=dmod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=dmod(zj,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_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)
! write(iout,*)"c ", c(1,:), c(2,:), c(3,:)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
+ sss_ele_cut=sscale_ele(1.0d0/(rij*sigma(itypi,itypj)))
+ sss_ele_grad=sscagrad_ele(1.0d0/(rij*sigma(itypi,itypj)))
+! print *,sss_ele_cut,sss_ele_grad,&
+! 1.0d0/(rij),r_cut_ele,rlamb_ele
+ if (sss_ele_cut.le.0.0) cycle
! Calculate angle-dependent terms of energy and contributions to their
! derivatives.
call sc_angular
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,& !d
! " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2," fac",fac !d
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij
+ evdw=evdw+evdwij*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
+! print *,'before fac',fac,rij,evdwij
+ fac=fac+evdwij*sss_ele_grad/sss_ele_cut&
+ /sigma(itypi,itypj)*rij
+! print *,'grad part scale',fac, &
+! evdwij*sss_ele_grad/sss_ele_cut &
+! /sigma(itypi,itypj)*rij
! fac=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+! print *,'before sc_grad', gg(1),gg(2),gg(3)
! Calculate angular part of the gradient.
call sc_grad
ENDIF ! dyn_ss
! if (.not.allocated(gradbx)) allocate(gradbx(3,nres)) !(3,maxres)
do i=ibondp_start,ibondp_end
+ if (itype(i-1).eq.ntyp1 .and. itype(i).eq.ntyp1) cycle
if (itype(i-1).eq.ntyp1 .or. itype(i).eq.ntyp1) then
- estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
- do j=1,3
- gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
- *dc(j,i-1)/vbld(i)
- enddo
- if (energy_dec) write(iout,*) &
- "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+!C estr1=estr1+gnmr1(vbld(i),-1.0d0,distchainmax)
+!C do j=1,3
+!C gradb(j,i-1)=gnmr1prim(vbld(i),-1.0d0,distchainmax) &
+!C *dc(j,i-1)/vbld(i)
+!C enddo
+!C if (energy_dec) write(iout,*) &
+!C "estr1",i,gnmr1(vbld(i),-1.0d0,distchainmax)
+ diff = vbld(i)-vbldpDUM
else
diff = vbld(i)-vbldp0
+ endif
if (energy_dec) write (iout,'(a7,i5,4f7.3)') &
"estr bb",i,vbld(i),vbldp0,diff,AKP*diff*diff
estr=estr+diff*diff
gradb(j,i-1)=AKP*diff*dc(j,i-1)/vbld(i)
enddo
! write (iout,'(i5,3f10.5)') i,(gradb(j,i-1),j=1,3)
- endif
+! endif
enddo
estr=0.5d0*AKP*estr+estr1
!
! " sigder",sigder
! write (iout,*) "eps1_om12",eps1_om12," eps2rt_om12",eps2rt_om12
! write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12
+!C print *,sss_ele_cut,'in sc_grad'
do k=1,3
dcosom1(k)=rij*(dc_norm(k,nres+i)-om1*erij(k))
dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
enddo
do k=1,3
- gg(k)=gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*sss_ele_cut
+!C print *,'gg',k,gg(k)
enddo
! write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
gvdwx(k,i)=gvdwx(k,i)-gg(k) &
+(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
- +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv &
+ *sss_ele_cut
+
gvdwx(k,j)=gvdwx(k,j)+gg(k) &
+(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
- +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv &
+ *sss_ele_cut
+
! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i)) &
! +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j)) &
! call intcartderiv
! call checkintcartgrad
call zerograd
- aincr=1.0D-4
+ aincr=1.0D-6
write(iout,*) 'Calling CHECK_ECARTINT.'
nf=0
icall=0
endif
return
end function sscale
+!!!!!!!!!! PBCSCALE
+ real(kind=8) function sscale_ele(r)
+! include "COMMON.SPLITELE"
+ real(kind=8) :: r,gamm
+ if(r.lt.r_cut_ele-rlamb_ele) then
+ sscale_ele=1.0d0
+ else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+ gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+ sscale_ele=1.0d0+gamm*gamm*(2*gamm-3.0d0)
+ else
+ sscale_ele=0d0
+ endif
+ return
+ end function sscale_ele
+
+ real(kind=8) function sscagrad_ele(r)
+ real(kind=8) :: r,gamm
+! include "COMMON.SPLITELE"
+ if(r.lt.r_cut_ele-rlamb_ele) then
+ sscagrad_ele=0.0d0
+ else if(r.le.r_cut_ele.and.r.ge.r_cut_ele-rlamb_ele) then
+ gamm=(r-(r_cut_ele-rlamb_ele))/rlamb_ele
+ sscagrad_ele=gamm*(6*gamm-6.0d0)/rlamb_ele
+ else
+ sscagrad_ele=0.0d0
+ endif
+ return
+ end function sscagrad_ele
+!!!!!!!!!!!!!!!
!-----------------------------------------------------------------------------
subroutine elj_long(evdw)
!
! include 'COMMON.CONTROL'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ 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) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
+
evdw=0.0D0
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
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)
+! Searching for nearest neighbour
+ 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-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)))
-
+ 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
if (sss.lt.1.0d0) then
! Calculate angle-dependent terms of energy and contributions to their
! write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
! & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
- evdw=evdw+evdwij*(1.0d0-sss)
+ evdw=evdw+evdwij*(1.0d0-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=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
! include 'COMMON.CONTROL'
logical :: lprn
!el local variables
- integer :: iint,itypi,itypi1,itypj
+ 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) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,&
+ dist_temp, dist_init
evdw=0.0D0
!cccc energy_dec=.false.
! print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
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)
+! dsci_inv=dsc_inv(itypi)
+ dsci_inv=vbld_inv(i+nres)
+
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)-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)
+! Searching for nearest neighbour
+ 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-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)))
+ 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
if (sss.gt.0.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=0.0d0
! Calculate the radial part of the gradient
gg(1)=xj*fac
dcosom2(k)=rij*(dc_norm(k,nres+j)-om2*erij(k))
enddo
do k=1,3
- gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac
+ gg(k)=(gg(k)+eom1*dcosom1(k)+eom2*dcosom2(k))*scalfac&
+ *sss_ele_cut
enddo
! write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
gvdwx(k,i)=gvdwx(k,i)-gg(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
+ +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac&
+ *sss_ele_cut
gvdwx(k,j)=gvdwx(k,j)+gg(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
+ +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac&
+ *sss_ele_cut
! write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
! & +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
! write (iout,*)(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
projects / unres4.git / blobdiff