Cleanup: function exc_t rather useless

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@11086 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2014-07-15 14:50:37 +00:00
parent c4aacc31b3
commit 903200ce58
7 changed files with 31 additions and 69 deletions

View File

@ -792,7 +792,6 @@ CONTAINS
REAL(dp), PARAMETER :: rho_eq_0(ndmx) = ZERO ! ccharge=0 when nlcc=.f.
!
REAL(dp) :: &
exc_t, & ! exchange-correlation function
eh, exc, edc, & ! hartree, xc and double counting energies
eloc, & ! local energy
rhovtot(ndmx), & ! total valence charge
@ -851,12 +850,12 @@ CONTAINS
rh(is) = vcharge_(i,is)/pawset_%grid%r2(i)/FPI
ENDDO
IF (nlcc_) rhc = ccharge_(i)/pawset_%grid%r2(i)/FPI
CALL vxc_t(rh,rhc,lsd,vxcr)
CALL vxc_t(lsd,rh,rhc,exc,vxcr)
vxc(i,1:nspin_)=vxcr(1:nspin_)
IF (nlcc_) THEN
aux(i)=exc_t(rh,rhc,lsd) * (rhovtot(i)+ccharge_(i))
aux(i)=exc * (rhovtot(i)+ccharge_(i))
ELSE
aux(i)=exc_t(rh,rhc,lsd) * rhovtot(i)
aux(i)=exc * rhovtot(i)
END IF
END DO
IF (dft_is_gradient()) THEN

View File

@ -26,7 +26,7 @@ subroutine elsdps
use funct, only : dft_is_gradient
implicit none
real(DP) :: &
exc_t, & ! exchange-correlation function
excc, vxcc(2), & ! exch-corr energy from core charge
int_0_inf_dr, & ! the integral function
rh0(2), & ! the charge in a given point
rhc, & ! core charge in a given point
@ -67,7 +67,8 @@ subroutine elsdps
rh0(2)=0.0_DP
do i=1,grid%mesh
rhc= rhoc(i)/grid%r2(i)/fpi
exccc(i) = exc_t(rh0,rhc,lsd)*rhoc(i)
call vxc_t(lsd,rh0,rhc,excc,vxcc)
exccc(i) = excc*rhoc(i)
enddo
if (dft_is_gradient()) then
allocate(rho_aux(ndmx,2), stat=ierr)

View File

@ -22,7 +22,7 @@ subroutine elsdps_paw( )
use funct, only : dft_is_gradient
implicit none
real(DP) :: &
exc_t, & ! exchange-correlation function
excc, vxcc(2), & ! exch-corr energy from core charge
int_0_inf_dr, & ! the integral function
rh0(2), & ! the charge in a given point
rhc, & ! core charge in a given point
@ -52,7 +52,8 @@ subroutine elsdps_paw( )
rh0(2)=0.0_DP
do i=1,grid%mesh
rhc= rhoc(i)/grid%r2(i)/fpi
exccc(i) = exc_t(rh0,rhc,lsd)*rhoc(i)
call vxc_t(lsd,rh0,rhc,excc,vxcc)
exccc(i) = excc*rhoc(i)
enddo
if (dft_is_gradient()) then
allocate(rho_aux(ndmx,2), stat=ierr)

View File

@ -23,7 +23,7 @@ subroutine new_potential &
logical :: nlcc, gga, oep, meta
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)
real(DP):: zed,enne,rh(2),rhc, exc_t
real(DP):: zed,enne,rh(2),rhc, excp
real(DP),allocatable:: vgc(:,:), egc(:), rhotot(:)
! real(DP),allocatable:: vx(:,:)
real(DP),allocatable:: dchi0(:,:)
@ -66,8 +66,8 @@ subroutine new_potential &
vxcp(:)=0.0_dp
exc(i) =0.0_dp
else
call vxc_t(rh,rhc,lsd,vxcp)
exc(i)=exc_t(rh,rhc,lsd)
call vxc_t(lsd,rh,rhc,excp,vxcp)
exc(i)=excp
endif
do is=1,nspin
vxc(i,is)=vxcp(is)

View File

@ -20,12 +20,12 @@ subroutine sic_correction(n,vhn1,vhn2,egc)
use radial_grids, only: hartree
implicit none
integer :: n
real(DP):: vhn1(ndmx),vhn2(ndmx), egc(ndmx), exc_t
real(DP):: vhn1(ndmx),vhn2(ndmx), egc(ndmx)
REAL(dp) :: & ! compatibility with metaGGA - not yet used
tau(ndmx) = 0.0_dp, vtau(ndmx) = 0.0_dp
!
integer :: i
real(DP):: rh(2), rhc, vxcp(2)
real(DP):: rh(2), rhc, vxcp(2), excp
real(DP):: vgc(ndmx,2), egc0(ndmx), rhotot(ndmx,2)
logical :: gga
@ -58,9 +58,9 @@ subroutine sic_correction(n,vhn1,vhn2,egc)
vhn1(i) = e2*vhn1(i)
rh(1) = rhotot(i,1)/grid%r2(i)/fpi
if (nlcc) rhc = rhoc(i)/grid%r2(i)/fpi
call vxc_t(rh,rhc,lsd,vxcp)
call vxc_t(lsd,rh,rhc,excp,vxcp)
vhn2(i)= vhn1(i)+vxcp(1)
egc(i)= exc_t(rh,rhc,lsd)*rhotot(i,1)
egc(i)= excp*rhotot(i,1)
end do
if (.not.gga) return

View File

@ -30,7 +30,7 @@ subroutine v_of_rho_at (rho,rhoc,vh,vxc,exc,excgga,vnew,nlcc,iflag)
!
logical :: gga, oep
integer :: i,is,nu,ierr
real(DP):: vxcp(2),rh(2),rhc, exc_t
real(DP):: vxcp(2),rh(2),rhc, excp
real(DP),allocatable:: vgc(:,:), egc(:), rhotot(:)
real(DP),allocatable:: dchi0(:,:)
@ -59,11 +59,11 @@ subroutine v_of_rho_at (rho,rhoc,vh,vxc,exc,excgga,vnew,nlcc,iflag)
rh(is) = rho(i,is)/grid%r2(i)/fpi
end do
if (nlcc) rhc = rhoc(i)/grid%r2(i)/fpi
call vxc_t(rh,rhc,lsd,vxcp)
call vxc_t(lsd,rh,rhc,excp,vxcp)
do is=1,nspin
vxc(i,is)=vxcp(is)
end do
exc(i)=exc_t(rh,rhc,lsd)
exc(i)=excp
end do
!
! if gga add gga exchange and correlation potential

View File

@ -7,23 +7,24 @@
!
!
!---------------------------------------------------------------
subroutine vxc_t(rho,rhoc,lsd,vxc)
subroutine vxc_t(lsd,rho,rhoc,exc,vxc)
!---------------------------------------------------------------
!
! this function returns the XC potential in LDA or LSDA approximation
! this function returns the XC potential and energy in LDA or
! LSDA approximation
!
use io_global, only : stdout
use kinds, only : DP
use funct, only : xc, xc_spin
implicit none
integer:: lsd
real(DP):: vxc(2), rho(2),rhoc,arho,zeta
real(DP):: vx(2), vc(2), ex, ec
implicit none
integer, intent(in) :: lsd
real(DP), intent(in) :: rho(2), rhoc
real(DP), intent(out):: exc, vxc(2)
real(DP):: arho, zeta, vx(2), vc(2), ex, ec
!
real(DP), parameter :: e2=2.0_dp, eps=1.e-30_dp
vxc(1)=0.0_dp
if (lsd.eq.1) vxc(2)=0.0_dp
exc=0.0_dp
if (lsd.eq.0) then
!
@ -33,11 +34,13 @@ subroutine vxc_t(rho,rhoc,lsd,vxc)
if (arho.gt.eps) then
call xc(arho,ex,ec,vx(1),vc(1))
vxc(1)=e2*(vx(1)+vc(1))
exc =e2*(ex+ec)
endif
else
!
! LSDA case
!
vxc(2)=0.0_dp
arho = abs(rho(1)+rho(2)+rhoc)
if (arho.gt.eps) then
zeta = (rho(1)-rho(2)) / arho
@ -47,54 +50,12 @@ subroutine vxc_t(rho,rhoc,lsd,vxc)
call xc_spin(arho,zeta,ex,ec,vx(1),vx(2),vc(1),vc(2))
vxc(1) = e2*(vx(1)+vc(1))
vxc(2) = e2*(vx(2)+vc(2))
exc = e2*(ex+ec)
endif
endif
return
end subroutine vxc_t
!---------------------------------------------------------------
function exc_t(rho,rhoc,lsd)
!---------------------------------------------------------------
!
use kinds, only : DP
use funct, only : xc, xc_spin
implicit none
integer:: lsd
real(DP) :: exc_t, rho(2),arho,rhot, zeta,rhoc
real(DP) :: ex, ec, vx(2), vc(2)
real(DP),parameter:: e2 =2.0_DP
exc_t=0.0_DP
if(lsd == 0) then
!
! LDA case
!
rhot = rho(1) + rhoc
arho = abs(rhot)
if (arho.gt.1.e-30_DP) then
call xc(arho,ex,ec,vx(1),vc(1))
exc_t=e2*(ex+ec)
endif
else
!
! LSDA case
!
rhot = rho(1)+rho(2)+rhoc
arho = abs(rhot)
if (arho.gt.1.e-30_DP) then
zeta = (rho(1)-rho(2)) / arho
! In atomic this cannot happen, but in PAW zeta can become
! a little larger than 1, or smaller than -1:
if( abs(zeta) > 1._dp) zeta = sign(1._dp, zeta)
call xc_spin(arho,zeta,ex,ec,vx(1),vx(2),vc(1),vc(2))
exc_t=e2*(ex+ec)
endif
endif
return
end function exc_t
!
!---------------------------------------------------------------
subroutine vxcgc ( ndm, mesh, nspin, r, r2, rho, rhoc, vgc, egc, &