TABLE OF CONTENTS


ABINIT/xchelu [ Functions ]

[ Top ] [ 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