mirror of https://gitlab.com/QEF/q-e.git
A tentative to improve the GGA instability close to the origin.
git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@4144 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
parent
99e6e79569
commit
eb7580e5f0
|
@ -87,7 +87,7 @@ subroutine c6_dft (mesh, zed, grid)
|
|||
end if
|
||||
|
||||
rhoc1=0.d0
|
||||
call new_potential(ndmx,mesh,grid,zed,vxt,lsd,.false.,latt,enne,rhoc1,rho,vh,vnew)
|
||||
call new_potential(ndmx,mesh,grid,zed,vxt,lsd,.false.,latt,enne,rhoc1,rho,vh,vnew,0)
|
||||
error = 0.d0
|
||||
do i=1,mesh
|
||||
error = error + abs( vpot(i,1)-vnew(i,1) ) * grid%r2(i) * grid%dx
|
||||
|
|
|
@ -82,7 +82,7 @@ subroutine descreening
|
|||
call chargeps(rhos,phits,nwfts,llts,jjts,octs,iwork)
|
||||
|
||||
call new_potential(ndmx,grid%mesh,grid,0.0_dp,vxt,lsd,nlcc,latt,enne,&
|
||||
rhoc,rhos,vh,vaux)
|
||||
rhoc,rhos,vh,vaux,1)
|
||||
|
||||
do n=1,grid%mesh
|
||||
vpstot(n,1)=vpsloc(n)
|
||||
|
|
|
@ -44,7 +44,7 @@ subroutine elsd (mesh,zed,grid,rho,zeta,vxt,vh,nlcc,nwf,enl,ll,lsd,&
|
|||
|
||||
rhoc=0.0_DP
|
||||
if (mesh.ne.grid%mesh) call errore('elsd','mesh dimension is not as expected',1)
|
||||
call vxcgc(ndm,mesh,nspin,grid%r,grid%r2,rho,rhoc,vgc,egc)
|
||||
call vxcgc(ndm,mesh,nspin,grid%r,grid%r2,rho,rhoc,vgc,egc,0)
|
||||
|
||||
rhc=0.0_DP
|
||||
do i=1,mesh
|
||||
|
|
|
@ -62,10 +62,11 @@ implicit none
|
|||
egc=0.0_DP
|
||||
if (gga.and.nlcc) then
|
||||
f1=0.0_DP
|
||||
call vxcgc(ndmx,grid%mesh,nspin,grid%r,grid%r2,f1,rhoc,vgc,egcc)
|
||||
call vxcgc(ndmx,grid%mesh,nspin,grid%r,grid%r2,f1,rhoc,vgc,egcc,1)
|
||||
endif
|
||||
|
||||
if (gga) call vxcgc(ndmx,grid%mesh,nspin,grid%r,grid%r2,rhos,rhoc,vgc,egc)
|
||||
if (gga) call vxcgc(ndmx,grid%mesh,nspin,grid%r,grid%r2,rhos, &
|
||||
rhoc,vgc,egc,1)
|
||||
|
||||
rh0(1)=0.0_DP
|
||||
rh0(2)=0.0_DP
|
||||
|
|
|
@ -8,7 +8,7 @@
|
|||
!
|
||||
!---------------------------------------------------------------
|
||||
subroutine new_potential &
|
||||
(ndm,mesh,grid,zed,vxt,lsd,nlcc,latt,enne,rhoc,rho,vh,vnew)
|
||||
(ndm,mesh,grid,zed,vxt,lsd,nlcc,latt,enne,rhoc,rho,vh,vnew,iflag)
|
||||
!---------------------------------------------------------------
|
||||
! set up the selfconsistent atomic potential
|
||||
!
|
||||
|
@ -19,6 +19,7 @@ subroutine new_potential &
|
|||
use ld1inc, only : nwf, vx
|
||||
implicit none
|
||||
type(radial_grid_type),intent(in):: grid
|
||||
integer, intent(in) :: iflag
|
||||
logical :: nlcc, gga, oep
|
||||
integer :: ndm,mesh,lsd,latt,i,is,nu, nspin, ierr
|
||||
real(DP):: rho(ndm,2),vxcp(2),vnew(ndm,2),vxt(ndm),vh(ndm), rhoc(ndm)
|
||||
|
@ -27,6 +28,7 @@ subroutine new_potential &
|
|||
! real(DP),allocatable:: vx(:,:)
|
||||
real(DP),allocatable:: dchi0(:,:)
|
||||
|
||||
|
||||
if (mesh.ne.grid%mesh) call errore('new_potential','mesh dimension is not as expected',1)
|
||||
gga=dft_is_gradient()
|
||||
oep=get_iexch().eq.4
|
||||
|
@ -70,7 +72,7 @@ subroutine new_potential &
|
|||
allocate(egc(ndm),stat=ierr)
|
||||
call errore('new_potential','allocating vgc and egc',ierr)
|
||||
|
||||
call vxcgc(ndm,mesh,nspin,grid%r,grid%r2,rho,rhoc,vgc,egc)
|
||||
call vxcgc(ndm,mesh,nspin,grid%r,grid%r2,rho,rhoc,vgc,egc,iflag)
|
||||
do is=1,nspin
|
||||
do i=1,mesh
|
||||
vnew(i,is)=vnew(i,is)+vgc(i,is)
|
||||
|
|
|
@ -75,7 +75,7 @@ subroutine scf
|
|||
! calculate new potential
|
||||
!
|
||||
call new_potential(ndmx,grid%mesh,grid,zed,vxt,&
|
||||
lsd,.false.,latt,enne,rhoc1,rho,vh,vnew)
|
||||
lsd,.false.,latt,enne,rhoc1,rho,vh,vnew,0)
|
||||
!
|
||||
! calculate SIC correction potential (if present)
|
||||
!
|
||||
|
@ -97,7 +97,7 @@ subroutine scf
|
|||
call vpack(grid%mesh,ndmx,nspin,vnew,vpot,1)
|
||||
call dmixp(grid%mesh*nspin,vnew,vpot,beta,tr2,iter,id,eps0,conv)
|
||||
call vpack(grid%mesh,ndmx,nspin,vnew,vpot,-1)
|
||||
! write(6,*) iter, eps0
|
||||
! write(6,*) iter, eps0
|
||||
!
|
||||
if (conv) then
|
||||
if (nerr /= 0) call infomsg ('scf','errors in KS equations')
|
||||
|
|
|
@ -124,23 +124,29 @@ subroutine set_rho_core
|
|||
write(stdout,110) grid%r(ik), a, b
|
||||
110 format (5x, ' r < ',f4.2,' : rho core = a sin(br)/r', &
|
||||
' a=',f7.2,' b=',f7.2/)
|
||||
1100 continue
|
||||
if (file_core .ne. ' ') then
|
||||
write(stdout,*) '***Writing file ',trim(file_core),' ***'
|
||||
write(stdout,'(6x, "***Writing file ",a, " ***")') trim(file_core)
|
||||
if (ionode) &
|
||||
open(unit=26,file=file_core, status='unknown', iostat=ios, err=300 )
|
||||
300 call mp_bcast(ios, ionode_id)
|
||||
call errore('set_rho_core','opening file '//file_core,abs(ios))
|
||||
if (ionode) then
|
||||
do n=1,grid%mesh
|
||||
write(26,'(4f20.10)') grid%r(n),rhoc(n),rhov(n),rhoco(n)
|
||||
enddo
|
||||
if (totrho>1.d-6) then
|
||||
do n=1,grid%mesh
|
||||
write(26,'(4f20.10)') grid%r(n),rhoc(n),rhov(n),rhoco(n)
|
||||
enddo
|
||||
else
|
||||
do n=1,grid%mesh
|
||||
write(26,'(2f20.10)') grid%r(n),rhov(n)
|
||||
enddo
|
||||
endif
|
||||
close(26)
|
||||
endif
|
||||
endif
|
||||
totrho = int_0_inf_dr(rhoc,grid,grid%mesh,2)
|
||||
write(stdout,'(13x,''integrated core pseudo-charge : '',f6.2)') totrho
|
||||
write(stdout,'(6x,''Integrated core pseudo-charge : '',f6.2)') totrho
|
||||
if (.not.nlcc) rhoc(1:grid%mesh) = 0.0_dp
|
||||
1100 continue
|
||||
deallocate (rhoco, rhov)
|
||||
return
|
||||
end subroutine set_rho_core
|
||||
|
|
|
@ -64,7 +64,7 @@ subroutine sic_correction(n,vhn1,vhn2,egc)
|
|||
! add gradient-correction terms to exchange-correlation potential
|
||||
!
|
||||
egc0=egc
|
||||
call vxcgc(ndmx,grid%mesh,nspin,grid%r,grid%r2,rhotot,rhoc,vgc,egc)
|
||||
call vxcgc(ndmx,grid%mesh,nspin,grid%r,grid%r2,rhotot,rhoc,vgc,egc,1)
|
||||
do i=1,grid%mesh
|
||||
vhn2(i)=vhn2(i)+vgc(i,1)
|
||||
egc(i)=egc(i)*grid%r2(i)*fpi+egc0(i)
|
||||
|
|
216
atomic/vxcgc.f90
216
atomic/vxcgc.f90
|
@ -7,11 +7,11 @@
|
|||
!
|
||||
!
|
||||
!---------------------------------------------------------------
|
||||
subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc)
|
||||
subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc,iflag)
|
||||
!---------------------------------------------------------------
|
||||
!
|
||||
!
|
||||
! This routine compute the exchange and correlation potential and
|
||||
! This routine computes the exchange and correlation potential and
|
||||
! energy to be added to the local density, to have the first
|
||||
! gradient correction.
|
||||
! In input the density is rho(r) (multiplied by 4*pi*r2).
|
||||
|
@ -21,25 +21,26 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc)
|
|||
use kinds, only : DP
|
||||
use constants, only : fpi
|
||||
use funct, only : gcxc, gcx_spin, gcc_spin
|
||||
use ld1inc, only : grid
|
||||
implicit none
|
||||
integer :: ndm,mesh,nspin,ndm1
|
||||
real(DP) :: r(mesh), r2(mesh), rho(ndm,2), rhoc(ndm), &
|
||||
vgc(ndm,2), egc(ndm)
|
||||
integer, intent(in) :: ndm,mesh,nspin,iflag
|
||||
real(DP), intent(in) :: r(mesh), r2(mesh), rho(ndm,2), rhoc(ndm)
|
||||
real(DP), intent(out):: vgc(ndm,2), egc(ndm)
|
||||
|
||||
integer :: i, is, ierr
|
||||
real(DP) :: sx,sc,v1x,v2x,v1c,v2c
|
||||
real(DP) :: v1xup, v1xdw, v2xup, v2xdw, v1cup, v1cdw
|
||||
real(DP) :: segno, arho, grho2(2)
|
||||
real(DP) :: rh, zeta, grh2
|
||||
real(DP) :: segno, arho
|
||||
real(DP) :: rh, zeta, grh2, grho2(2)
|
||||
real(DP),parameter :: eps=1.e-12_dp
|
||||
|
||||
real(DP), pointer :: grho(:,:), h(:,:), dh(:)
|
||||
real(DP), allocatable :: grho(:,:), h(:,:), dh(:), rhoaux(:,:)
|
||||
!
|
||||
! First compute the charge and the charge gradient, assumed
|
||||
! to have spherical symmetry. The gradient is the derivative of
|
||||
! the charge with respect to the modulus of r. The last point is
|
||||
! assumed to have zero gradient as happens in an atom.
|
||||
! the charge with respect to the modulus of r.
|
||||
!
|
||||
allocate(rhoaux(mesh,2),stat=ierr)
|
||||
allocate(grho(mesh,2),stat=ierr)
|
||||
allocate(h(mesh,2),stat=ierr)
|
||||
allocate(dh(mesh),stat=ierr)
|
||||
|
@ -49,21 +50,9 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc)
|
|||
|
||||
do is=1,nspin
|
||||
do i=1, mesh
|
||||
rho(i,is)=(rho(i,is)+rhoc(i)/nspin)/fpi/r2(i)
|
||||
rhoaux(i,is)=(rho(i,is)+rhoc(i)/nspin)/fpi/r2(i)
|
||||
enddo
|
||||
do i=2, mesh-1
|
||||
grho(i,is)=( (r(i+1)-r(i))**2*(rho(i-1,is)-rho(i,is)) &
|
||||
-(r(i-1)-r(i))**2*(rho(i+1,is)-rho(i,is)) ) &
|
||||
/((r(i+1)-r(i))*(r(i-1)-r(i))*(r(i+1)-r(i-1)))
|
||||
enddo
|
||||
grho(mesh,is)=0.0_dp
|
||||
!
|
||||
! The gradient in the first point is a linear interpolation of the
|
||||
! gradient at point 2 and 3. The final result is not really sensitive to
|
||||
! the value of these derivatives.
|
||||
!
|
||||
grho(1,is)=grho(2,is)+(grho(3,is)-grho(2,is)) &
|
||||
*(r(1)-r(2))/(r(3)-r(2))
|
||||
call gradient(rhoaux(1,is),grho(1,is),r,mesh,iflag)
|
||||
enddo
|
||||
|
||||
if (nspin.eq.1) then
|
||||
|
@ -71,8 +60,8 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc)
|
|||
! GGA case
|
||||
!
|
||||
do i=1,mesh
|
||||
arho=abs(rho(i,1))
|
||||
segno=sign(1.0_dp,rho(i,1))
|
||||
arho=abs(rhoaux(i,1))
|
||||
segno=sign(1.0_dp,rhoaux(i,1))
|
||||
if (arho.gt.eps.and.abs(grho(i,1)).gt.eps) then
|
||||
call gcxc(arho,grho(i,1)**2,sx,sc,v1x,v2x,v1c,v2c)
|
||||
egc(i)=(sx+sc)*segno
|
||||
|
@ -97,21 +86,17 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc)
|
|||
! or gradients are zero or negative must
|
||||
! be detected within the gcxc_spin routine
|
||||
!
|
||||
! call gcxc_spin(rho(i,1),rho(i,2),grho(i,1),grho(i,2), &
|
||||
! sx,sc,v1xup,v1xdw,v2xup,v2xdw, &
|
||||
! v1cup,v1cdw,v2c)
|
||||
!
|
||||
! spin-polarised case
|
||||
!
|
||||
do is = 1, nspin
|
||||
grho2(is)=grho(i,is)**2
|
||||
enddo
|
||||
|
||||
call gcx_spin (rho(i, 1), rho(i, 2), grho2(1), grho2(2), &
|
||||
call gcx_spin (rhoaux(i, 1), rhoaux(i, 2), grho2(1), grho2(2), &
|
||||
sx, v1xup, v1xdw, v2xup, v2xdw)
|
||||
rh = rho(i, 1) + rho(i, 2)
|
||||
rh = rhoaux(i, 1) + rhoaux(i, 2)
|
||||
if (rh.gt.eps) then
|
||||
zeta = (rho (i, 1) - rho (i, 2) ) / rh
|
||||
zeta = (rhoaux (i, 1) - rhoaux (i, 2) ) / rh
|
||||
grh2 = (grho (i, 1) + grho (i, 2) ) **2
|
||||
call gcc_spin (rh, zeta, grh2, sc, v1cup, v1cdw, v2c)
|
||||
else
|
||||
|
@ -136,15 +121,7 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc)
|
|||
! and correlation potential.
|
||||
!
|
||||
do is=1,nspin
|
||||
do i=2,mesh-1
|
||||
dh(i)=( (r(i+1)-r(i))**2*(h(i-1,is)-h(i,is)) &
|
||||
-(r(i-1)-r(i))**2*(h(i+1,is)-h(i,is)) ) &
|
||||
/( (r(i+1)-r(i))*(r(i-1)-r(i))*(r(i+1)-r(i-1)) )
|
||||
enddo
|
||||
|
||||
dh(1)=dh(2)+(dh(3)-dh(2)) &
|
||||
*(r(1)-r(2))/(r(3)-r(2))
|
||||
dh(mesh)=0.0_dp
|
||||
call gradient(h(1,is),dh,r,mesh,iflag)
|
||||
!
|
||||
! Finally we compute the total exchange and correlation energy and
|
||||
! potential. We put the original values on the charge and multiply
|
||||
|
@ -152,7 +129,6 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc)
|
|||
|
||||
do i=1, mesh
|
||||
vgc(i,is)=vgc(i,is)-dh(i)/r2(i)
|
||||
rho(i,is)=rho(i,is)*fpi*r2(i)-rhoc(i)/nspin
|
||||
vgc(i,is)=2.0_dp*vgc(i,is)
|
||||
if (is.eq.1) egc(i)=2.0_dp*egc(i)
|
||||
! if (is.eq.1.and.i.lt.4) write(6,'(3f20.12)') &
|
||||
|
@ -163,6 +139,160 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc)
|
|||
deallocate(dh)
|
||||
deallocate(h)
|
||||
deallocate(grho)
|
||||
deallocate(rhoaux)
|
||||
|
||||
return
|
||||
end subroutine vxcgc
|
||||
|
||||
subroutine gradient(f,gf,r,mesh,iflag)
|
||||
!
|
||||
! This subroutine calculates the derivative with respect to r of a
|
||||
! radial function defined on the mesh r. If iflag=0 it uses all mesh
|
||||
! points. If iflag=1 it uses only a coarse grained mesh close to the
|
||||
! origin, to avoid large errors in the derivative when the function
|
||||
! is too smooth.
|
||||
!
|
||||
use kinds, only : DP
|
||||
use ld1inc, only: zed
|
||||
use radial_grids, only : series
|
||||
implicit none
|
||||
integer, intent(in) :: mesh, iflag
|
||||
real(DP), intent(in) :: f(mesh), r(mesh)
|
||||
real(DP), intent(out) :: gf(mesh)
|
||||
|
||||
integer :: i,j,k,imin,npoint
|
||||
real(DP) :: delta, b(5), faux(6), raux(6)
|
||||
!
|
||||
! This formula is used in the all-electron case.
|
||||
!
|
||||
if (iflag==0) then
|
||||
do i=2, mesh-1
|
||||
gf(i)=( (r(i+1)-r(i))**2*(f(i-1)-f(i)) &
|
||||
-(r(i-1)-r(i))**2*(f(i+1)-f(i)) ) &
|
||||
/((r(i+1)-r(i))*(r(i-1)-r(i))*(r(i+1)-r(i-1)))
|
||||
enddo
|
||||
gf(mesh)=0.0_dp
|
||||
!
|
||||
! The gradient in the first point is a linear interpolation of the
|
||||
! gradient at point 2 and 3.
|
||||
!
|
||||
gf(1) = gf(2) + (gf(3)-gf(2)) * (r(1)-r(2)) / (r(3)-r(2))
|
||||
return
|
||||
endif
|
||||
!
|
||||
! If the input function is slowly changing (as the pseudocharge),
|
||||
! the previous formula is affected by numerical errors close to the
|
||||
! origin where the r points are too close one to the other. Therefore
|
||||
! we calculate the gradient on a coarser mesh. This gradient is often
|
||||
! more accurate but still does not remove all instabilities observed
|
||||
! with the GGA.
|
||||
! At larger r the distances between points become larger than delta
|
||||
! and this formula coincides with the previous one.
|
||||
! (ADC 08/2007)
|
||||
!
|
||||
|
||||
delta=0.00001_dp
|
||||
|
||||
imin=1
|
||||
points: do i=2, mesh
|
||||
do j=i+1,mesh
|
||||
if (r(j)>r(i)+delta) then
|
||||
do k=i-1,1,-1
|
||||
if (r(k)<r(i)-delta) then
|
||||
gf(i)=( (r(j)-r(i))**2*(f(k)-f(i)) &
|
||||
-(r(k)-r(i))**2*(f(j)-f(i)) ) &
|
||||
/((r(j)-r(i))*(r(k)-r(i))*(r(j)-r(k)))
|
||||
cycle points
|
||||
endif
|
||||
enddo
|
||||
!
|
||||
! if the code arrives here there are not enough points on the left:
|
||||
! r(i)-delta is smaller than r(1).
|
||||
!
|
||||
imin=i
|
||||
cycle points
|
||||
endif
|
||||
enddo
|
||||
!
|
||||
! If the code arrives here there are not enough points on the right.
|
||||
! It should happen only at mesh.
|
||||
! NB: the f function is assumed to be vanishing for large r, so the gradient
|
||||
! in the last points is taken as zero.
|
||||
!
|
||||
gf(i)=0.0_DP
|
||||
enddo points
|
||||
!
|
||||
! In the first imin points the previous formula cannot be
|
||||
! used. We interpolate with a polynomial the points already found
|
||||
! and extrapolate in the points from 1 to imin.
|
||||
! Presently we fit 5 points with a 3th degree polynomial.
|
||||
!
|
||||
npoint=5
|
||||
raux=0.0_DP
|
||||
faux=0.0_DP
|
||||
faux(1)=gf(imin+1)
|
||||
raux(1)=r(imin+1)
|
||||
j=imin+1
|
||||
points_fit: do k=2,npoint
|
||||
do i=j,mesh-1
|
||||
if (r(i)>r(imin+1)+(k-1)*delta) then
|
||||
faux(k)=gf(i)
|
||||
raux(k)=r(i)
|
||||
j=i+1
|
||||
cycle points_fit
|
||||
endif
|
||||
enddo
|
||||
enddo points_fit
|
||||
call fit_pol(raux,faux,npoint,3,b)
|
||||
do i=1,imin
|
||||
gf(i)=b(1)+r(i)*(b(2)+r(i)*(b(3)+r(i)*b(4)))
|
||||
enddo
|
||||
return
|
||||
end subroutine gradient
|
||||
|
||||
subroutine fit_pol(xdata,ydata,n,degree,b)
|
||||
!
|
||||
! This routine finds the coefficients of the least-square polynomial which
|
||||
! interpolates the n input data points.
|
||||
!
|
||||
use kinds, ONLY : DP
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: n, degree
|
||||
real(DP), intent(in) :: xdata(n), ydata(n)
|
||||
real(DP), intent(out) :: b(degree+1)
|
||||
|
||||
integer :: ipiv(degree+1), info, i, j, k
|
||||
real(DP) :: bmat(degree+1,degree+1), amat(degree+1,n)
|
||||
|
||||
amat(1,:)=1.0_DP
|
||||
do i=2,degree+1
|
||||
do j=1,n
|
||||
amat(i,j)=amat(i-1,j)*xdata(j)
|
||||
enddo
|
||||
enddo
|
||||
do i=1,degree+1
|
||||
b(i)=0.0_DP
|
||||
do k=1,n
|
||||
b(i)=b(i)+ydata(k)*xdata(k)**(i-1)
|
||||
enddo
|
||||
enddo
|
||||
do i=1,degree+1
|
||||
do j=1,degree+1
|
||||
bmat(i,j)=0.0_DP
|
||||
do k=1,n
|
||||
bmat(i,j)=bmat(i,j)+amat(i,k)*amat(j,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! This lapack routine solves the linear system that gives the
|
||||
! coefficients of the interpolating polynomial.
|
||||
!
|
||||
call DGESV(degree+1, 1, bmat, degree+1, ipiv, b, degree+1, info)
|
||||
|
||||
if (info.ne.0) call errore('pol_fit','problems with the linear system', &
|
||||
abs(info))
|
||||
return
|
||||
end subroutine fit_pol
|
||||
|
||||
|
|
|
@ -3,7 +3,7 @@
|
|||
title='H',
|
||||
xmin=-6,
|
||||
rel=0,
|
||||
lsd=1,
|
||||
lsd=0,
|
||||
isic=1,
|
||||
iswitch=1,
|
||||
config='1s1'
|
||||
|
|
Loading…
Reference in New Issue