TABLE OF CONTENTS


ABINIT/xcwign [ Functions ]

[ Top ] [ Functions ]

NAME

 xcwign

FUNCTION

 Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input $\rho$.
 Wigner exchange and correlation (xc)--see e.g. David Pines,
 Elementary Excitations in Solids, p. 94, NY 1964.
 Expression is exc=-(0.44)/(rs+7.8)-efac/rs (hartree), efac below.
 rs = $(3/(4\pi))^{1/3}* \rho (r)^{-1/3}$.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR)
 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)=corresponding Wigner-Seitz radii, precomputed

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 xcwign(exc,npt,order,rspts,vxc,& !Mandatory arguments
 45 &                dvxc)                           !Optional arguments
 46 
 47  use defs_basis
 48  use m_errors
 49  use m_profiling_abi
 50 
 51 !This section has been created automatically by the script Abilint (TD).
 52 !Do not modify the following lines by hand.
 53 #undef ABI_FUNC
 54 #define ABI_FUNC 'xcwign'
 55 !End of the abilint section
 56 
 57  implicit none
 58 
 59 !Arguments ------------------------------------
 60 !scalars
 61  integer,intent(in) :: npt,order
 62 !arrays
 63  real(dp),intent(in) :: rspts(npt)
 64  real(dp),intent(out) :: exc(npt),vxc(npt)
 65  real(dp),intent(out),optional :: dvxc(npt)
 66 
 67 !Local variables-------------------------------
 68 !c1 and c2 are the Wigner parameters in hartree and bohr resp.
 69 !scalars
 70  integer :: ipt
 71  real(dp),parameter :: c1=0.44_dp,c2=7.8_dp,c4_3=4.0_dp/3.0_dp
 72  real(dp),parameter :: c8_27=8.0_dp/27.0_dp
 73  real(dp) :: dfac,efac,rs,rsc2m1,rsm1,vfac,vxcnum
 74  character(len=500) :: message
 75 
 76 ! *************************************************************************
 77 
 78 !Checks the values of order
 79  if(order<0 .or. order>2)then
 80    write(message, '(a,a,a,i0)' )&
 81 &   'With Wigner xc functional, the only',ch10,&
 82 &   'allowed values for order are 0, 1 or 2, while it is found to be',order
 83    MSG_BUG(message)
 84  end if
 85 
 86 !Checks the compatibility between the order and the presence of the optional arguments
 87  if(order <= 1 .and. present(dvxc))then
 88    write(message, '(a,a,a,i3)' )&
 89 &   'The order chosen does not need the presence',ch10,&
 90 &   'of the vector dvxc, that is needed only with order=2 , while we have',order
 91    MSG_BUG(message)
 92  end if
 93 
 94 !Compute vfac=(3/(2*Pi))^(2/3)
 95  vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
 96 !Compute efac=(3/4)*vfac
 97  efac=0.75_dp*vfac
 98 !Compute dfac=(4*Pi/9)*vfac
 99  dfac=(4.0_dp*pi/9.0_dp)*vfac
100 
101 !separate cases with respect to order
102  if (order==2) then
103    
104 !  Loop over grid points
105    do ipt=1,npt
106      rs=rspts(ipt)
107      rsm1=1.0_dp/rs
108      rsc2m1=1.0_dp/(rs+c2)
109 !    compute energy density (hartree)
110      exc(ipt)=-c1*rsc2m1-efac*rsm1
111      vxcnum=-(c4_3*rs+c2)*c1
112 !    compute potential (hartree)
113      vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1
114 !    compute d(vxc)/d(rho) (hartree*bohr^3)
115      dvxc(ipt)=-(c8_27*pi)*(c1*rs**4)*(rs+rs+c2)*rsc2m1**3-dfac*rs**2
116    end do
117  else
118    
119 !  Loop over grid points
120    do ipt=1,npt
121      rs=rspts(ipt)
122      rsm1=1.0_dp/rs
123      rsc2m1=1.0_dp/(rs+c2)
124 !    compute energy density (hartree)
125      exc(ipt)=-c1*rsc2m1-efac*rsm1
126      vxcnum=-(c4_3*rs+c2)*c1
127 !    compute potential (hartree)
128      vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1
129    end do
130    
131  end if
132 !
133 end subroutine xcwign