TABLE OF CONTENTS
ABINIT/xchelu [ Functions ]
NAME
xchelu
FUNCTION
Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input rho.
NOTES
Hedin-Lundqvist exchange and correlation (xc)-- L. Hedin and B.I. Lundqvist, J. Phys. C. 4, 2064 (1971).
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, LG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
npt=number of real space points on which density is provided order=gives the maximal derivative of Exc computed. rspts(npt)=Wigner-Seitz radii at each point
OUTPUT
exc(npt)=exchange-correlation energy density (hartree) vxc(npt)=xc potential (d($\rho$*exc)/d($\rho$)) (hartree) if(order>1) dvxc(npt)=derivative d(vxc)/d($\rho$) (hartree*bohr^3)
PARENTS
drivexc
CHILDREN
SOURCE
37 #if defined HAVE_CONFIG_H 38 #include "config.h" 39 #endif 40 41 #include "abi_common.h" 42 43 44 subroutine xchelu(exc,npt,order,rspts,vxc,dvxc) ! dvxc is optional 45 46 use defs_basis 47 use m_profiling_abi 48 use m_errors 49 50 !This section has been created automatically by the script Abilint (TD). 51 !Do not modify the following lines by hand. 52 #undef ABI_FUNC 53 #define ABI_FUNC 'xchelu' 54 !End of the abilint section 55 56 implicit none 57 58 !Arguments ------------------------------------ 59 !scalars 60 integer,intent(in) :: npt,order 61 !arrays 62 real(dp),intent(in) :: rspts(npt) 63 real(dp),intent(out) :: exc(npt),vxc(npt) 64 real(dp),intent(out),optional :: dvxc(npt) 65 66 !Local variables------------------------------- 67 !aa and cc are H-L fitting parameters A and C (C in hartree) 68 !rs = (3/(4 Pi))**(1/3) * rho(r)**(-1/3). 69 !scalars 70 integer :: ipt 71 real(dp),parameter :: aa=21_dp,c1_21=one/21_dp,c4_9=4.0_dp/9.0_dp,cc=0.0225_dp 72 real(dp) :: dfac,efac,rs,rsm1,vfac,xx 73 character(len=500) :: message 74 75 ! ************************************************************************* 76 77 !Checks the values of order 78 if(order<0 .or. order>2)then 79 write(message, '(a,a,a,i0)' )& 80 & 'With Hedin-Lundqvist xc functional, the only',ch10,& 81 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 82 MSG_BUG(message) 83 end if 84 85 !Compute vfac=(3/(2*Pi))^(2/3) 86 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 87 !Compute efac=(3/4)*vfac 88 efac=0.75_dp*vfac 89 !Compute dfac=(4*Pi/9)*vfac 90 dfac=(4.0_dp*pi/9.0_dp)*vfac 91 !separate cases with respect to order 92 if (order==2) then 93 ! Loop over grid points 94 do ipt=1,npt 95 rs=rspts(ipt) 96 rsm1=one/rs 97 ! compute energy density exc (hartree) 98 xx=rs*c1_21 99 exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+& 100 & half*xx-xx*xx-third) - efac*rsm1 101 ! compute xc potential d(rho*exc)/d(rho) (hartree) 102 vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1 103 ! compute d(vxc)/d(rho) (hartree*bohr^3) 104 dvxc(ipt)=-(rs**2)*((c4_9*pi)*cc*rs/(one+xx) + dfac) 105 end do 106 else 107 ! Loop over grid points 108 do ipt=1,npt 109 rs=rspts(ipt) 110 rsm1=one/rs 111 ! compute energy density exc (hartree) 112 xx=rs*c1_21 113 exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+& 114 & half*xx-xx*xx-third) - efac*rsm1 115 ! compute xc potential d(rho*exc)/d(rho) (hartree) 116 vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1 117 end do 118 end if 119 ! 120 end subroutine xchelu