TABLE OF CONTENTS


ABINIT/m_xclda [ Modules ]

[ Top ] [ Modules ]

NAME

  m_xclda

FUNCTION

  LDA or LDA-like XC functionals.

COPYRIGHT

  Copyright (C) 1998-2022 ABINIT group (DCA, XG, GMR, LG, MF, JFD, LK)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_xclda
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27 
28  use m_numeric_tools,      only : invcb
29 
30  implicit none
31 
32  private

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) [[cite:Hedin1971]]

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)

SOURCE

 988 subroutine xchelu(exc,npt,order,rspts,vxc,dvxc)  ! dvxc is optional
 989 
 990 !Arguments ------------------------------------
 991 !scalars
 992  integer,intent(in) :: npt,order
 993 !arrays
 994  real(dp),intent(in) :: rspts(npt)
 995  real(dp),intent(out) :: exc(npt),vxc(npt)
 996  real(dp),intent(out),optional :: dvxc(npt)
 997 
 998 !Local variables-------------------------------
 999 !aa and cc are H-L fitting parameters A and C (C in hartree)
1000 !rs = (3/(4 Pi))**(1/3) * rho(r)**(-1/3).
1001 !scalars
1002  integer :: ipt
1003  real(dp),parameter :: aa=21_dp,c1_21=one/21_dp,c4_9=4.0_dp/9.0_dp,cc=0.0225_dp
1004  real(dp) :: dfac,efac,rs,rsm1,vfac,xx
1005  character(len=500) :: message
1006 
1007 ! *************************************************************************
1008 
1009 !Checks the values of order
1010  if(order<0 .or. order>2)then
1011    write(message, '(a,a,a,i0)' )&
1012 &   'With Hedin-Lundqvist xc functional, the only',ch10,&
1013 &   'allowed values for order are 0, 1 or 2, while it is found to be',order
1014    ABI_BUG(message)
1015  end if
1016 
1017 !Compute vfac=(3/(2*Pi))^(2/3)
1018  vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
1019 !Compute efac=(3/4)*vfac
1020  efac=0.75_dp*vfac
1021 !Compute dfac=(4*Pi/9)*vfac
1022  dfac=(4.0_dp*pi/9.0_dp)*vfac
1023 !separate cases with respect to order
1024  if (order==2) then
1025 !  Loop over grid points
1026    do ipt=1,npt
1027      rs=rspts(ipt)
1028      rsm1=one/rs
1029 !    compute energy density exc (hartree)
1030      xx=rs*c1_21
1031      exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+&
1032 &     half*xx-xx*xx-third) - efac*rsm1
1033 !    compute xc potential d(rho*exc)/d(rho) (hartree)
1034      vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1
1035 !    compute d(vxc)/d(rho) (hartree*bohr^3)
1036      dvxc(ipt)=-(rs**2)*((c4_9*pi)*cc*rs/(one+xx) + dfac)
1037    end do
1038  else
1039 !  Loop over grid points
1040    do ipt=1,npt
1041      rs=rspts(ipt)
1042      rsm1=one/rs
1043 !    compute energy density exc (hartree)
1044      xx=rs*c1_21
1045      exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+&
1046 &     half*xx-xx*xx-third) - efac*rsm1
1047 !    compute xc potential d(rho*exc)/d(rho) (hartree)
1048      vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1
1049    end do
1050  end if
1051 !
1052 end subroutine xchelu

ABINIT/xclb [ Functions ]

[ Top ] [ Functions ]

NAME

 xclb

FUNCTION

 Computes the GGA like part (vx_lb) of the Leeuwen-Baerends
 exchange-correlation potential (vxc_lb) and adds it to the
 lda exchange-correlation potential (vxc_lda) which
 must be provided as input,
            vxci  <--  vxc_lb =: vxc_lda + vx_lb

 R van Leeuwen and EJ Baerends, Phys Rev A 49, 2421 (1994) [[cite:VanLeeuwen1994]]

 With respect to spin, the van Leeuwen-Baerends
 potential is "exchange-like" : separate contributions from
 spin up and spin down.

INPUTS

  npts= number of points to be computed
  nspden=1 for unpolarized, 2 for spin-polarized
  grho2_updn(npts,2*nspden-1)=square of the gradient of the spin-up,
     and, if nspden==2, spin-down, and total density (Hartree/Bohr**2)
  rho_updn(npts,nspden)=spin-up and spin-down density (Hartree/bohr**3)

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output:
  vxci(npts,nspden)=input xc potential to which Leeuwen-Baerends correction
   is added at output.

SOURCE

1175 subroutine xclb(grho2_updn,npts,nspden,rho_updn,vxci)
1176 
1177 !Arguments ------------------------------------
1178 !scalars
1179  integer,intent(in) :: npts,nspden
1180 !arrays
1181  real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden)
1182  real(dp),intent(inout) :: vxci(npts,nspden)
1183 
1184 !Local variables-------------------------------
1185 !scalars
1186  integer :: ipts,ispden
1187  real(dp),parameter :: beta=0.05_dp
1188  real(dp) :: density,density_gradient,density_t13,s_g_sq,scaled_gradient
1189  real(dp) :: scaling_factor,vx_lb
1190 
1191 ! *************************************************************************
1192 
1193 !DEBUG
1194 !write(std_out,*) ' %xclb: enter'
1195 !ENDDEBUG
1196 
1197 !scale the spin densities for evaluating spin up or down exchange
1198  scaling_factor=one
1199  if(nspden == 2) scaling_factor=two
1200 
1201  do ispden=1,nspden
1202 
1203    do ipts=1,npts
1204 
1205      density= scaling_factor * rho_updn(ipts,ispden)
1206      density_gradient= scaling_factor * sqrt(grho2_updn(ipts,ispden))
1207 
1208      density_t13= density**third
1209      scaled_gradient= density_gradient/max(density*density_t13,1.e-12_dp)
1210 
1211      s_g_sq= scaled_gradient*scaled_gradient
1212 
1213      vx_lb= -beta*density_t13 * s_g_sq/ &
1214 &     (one+3.d0*beta* scaled_gradient*log(scaled_gradient+sqrt(one+s_g_sq*s_g_sq)))
1215 
1216      vxci(ipts,ispden)=vxci(ipts,ispden)+vx_lb
1217    end do
1218 
1219  end do
1220 
1221 end subroutine xclb

ABINIT/xcpzca [ Functions ]

[ Top ] [ Functions ]

NAME

 xcpzca

FUNCTION

 Returns exc, vxc, and d(vxc)/d($\rho$) from input rho.

 NOTE
 Perdew-Zunger parameterization of Ceperly-Alder electron gas energy data.
 J. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981) [[cite:Perdew1981]]
 D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980) [[cite:Ceperley1980]]

INPUTS

  npt=number of real space points on which density is provided
  order=gives the maximal derivative of Exc computed.
  rhor(npt)=electron number density (bohr^-3)
  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)

SOURCE

 73 subroutine xcpzca(exc,npt,order,rhor,rspts,vxc,&  !Mandatory arguments
 74 &                dvxc)                            !Optional arguments
 75 
 76 !Arguments ------------------------------------
 77 !scalars
 78  integer,intent(in) :: npt,order
 79 !arrays
 80  real(dp),intent(in) :: rhor(npt),rspts(npt)
 81  real(dp),intent(out) :: exc(npt),vxc(npt)
 82  real(dp),intent(out),optional :: dvxc(npt)
 83 
 84 !Local variables-------------------------------
 85 !Perdew-Zunger parameters a, b, b1, b2, c, d, gamma
 86 !scalars
 87  integer :: ipt
 88  real(dp),parameter :: aa=0.0311_dp,b1=1.0529_dp,b2=0.3334_dp,bb=-0.048_dp
 89  real(dp),parameter :: c4_3=4.0_dp/3.0_dp,c7_6=7.0_dp/6.0_dp,cc=0.0020_dp
 90  real(dp),parameter :: dd=-0.0116_dp,ga=-0.1423_dp
 91  real(dp) :: den,den3,dfac,efac,logrs,rs,rsm1,t1,t2,vfac
 92  character(len=500) :: message
 93 
 94 ! *************************************************************************
 95 
 96 !Compute vfac=(3/(2*Pi))^(2/3)
 97  vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
 98 !Compute efac=(3/4)*vfac
 99  efac=0.75_dp*vfac
100 !Compute dfac=(4*Pi/9)*vfac
101  dfac=(4.0_dp*pi/9.0_dp)*vfac
102 
103 !Checks the values of order
104  if(order<0 .or. order>2)then
105    write(message, '(a,a,a,i0)' )&
106 &   'With Perdew-Zunger Ceperley-Alder xc functional, the only',ch10,&
107 &   'allowed values for order are 0, 1 or 2, while it is found to be',order
108    ABI_BUG(message)
109  end if
110 
111 !Checks the compatibility between the order and the presence of the optional arguments
112  if(order <= 1 .and. present(dvxc))then
113    write(message, '(a,a,a,i0)' )&
114 &   'The order chosen does not need the presence',ch10,&
115 &   'of the vector dvxc, that is needed only with order=2 , while we have',order
116    ABI_BUG(message)
117  end if
118 
119 !separate cases with respect to order
120  if(order==2) then
121 !  Loop over grid points
122    do ipt=1,npt
123      rs=rspts(ipt)
124      rsm1=1.0_dp/rs
125 !    Consider two regimes: rs<1 or rs>=1
126      if (rs<1._dp) then
127        logrs=log(rs)
128 !      compute energy density exc (hartree)
129        exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1
130 !      compute potential vxc=d(rho*exc)/d(rho) (hartree)
131        vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+&
132 &       (bb-aa*third)-vfac*rsm1
133 !      compute d(vxc)/d(rho) (hartree*bohr^3)
134        dvxc(ipt)=-(3._dp*aa+(cc+dd+dd)*rs+2._dp*cc*rs*logrs)&
135 &       /(9._dp*rhor(ipt))-dfac*rs**2
136      else if (rs<1000._dp) then
137        t1=b1*sqrt(rs)
138        t2=b2*rs
139        den=1._dp/(1._dp+t1+t2)
140        exc(ipt)=ga*den-efac*rsm1
141        vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
142        den3=den**3
143        dvxc(ipt)=(ga*den3/(36._dp*rhor(ipt)))*(5._dp*t1+8._dp*t2+&
144 &       7._dp*t1**2+16._dp*t2**2+21._dp*t1*t2)-dfac*rs**2
145      else
146        t1=b1*sqrt(rs)
147        t2=b2*rs
148        den=1._dp/(1._dp+t1+t2)
149        exc(ipt)=ga*den-efac*rsm1
150        vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
151        dvxc(ipt)=0._dp
152      end if
153    end do
154  else
155 !  Loop over grid points
156    do ipt=1,npt
157      rs=rspts(ipt)
158      rsm1=1.0_dp/rs
159 !    Consider two regimes: rs<1 or rs>=1
160      if (rs<1._dp) then
161        logrs=log(rs)
162 !      compute energy density exc (hartree)
163        exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1
164 !      compute potential vxc=d(rho*exc)/d(rho) (hartree)
165        vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+&
166 &       (bb-aa*third)-vfac*rsm1
167 !      compute d(vxc)/d(rho) (hartree*bohr^3)
168      else
169        t1=b1*sqrt(rs)
170        t2=b2*rs
171        den=1._dp/(1._dp+t1+t2)
172        exc(ipt)=ga*den-efac*rsm1
173        vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
174      end if
175    end do
176  end if
177 !
178 end subroutine xcpzca

ABINIT/xcspol [ Functions ]

[ Top ] [ Functions ]

NAME

 xcspol

FUNCTION

 Spin-polarized exchange and correlation, parameterized by Mike Teter of Corning Incorporated.

INPUTS

  nspden=number of spin-density components
  npts= number of points to be computed
  order=its absolute value gives the maximal derivative of Exc to be computed.
  rspts(npts)=Seitz electron radius (bohr)
  zeta(npts)=$(\rho\uparrow-\rho\downarrow)/(\rho\uparrow+\rho\downarrow)$=degree of polarization
  (ignored if nspden=1, in which case zeta should be 0)

OUTPUT

  if(abs(order)>1) dvxc(npts,1+nspden)=              (Hartree*bohr^3)
   if(nspden=1 .and. order==2): dvxc(:,1)=dvxc/d$\rho$ , dvxc(:,2) empty
   if(nspden=1 .and. order==-2): also compute dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$
   if(nspden=2): dvxc(:,1)=dvxc($\uparrow$)/d$\rho(\uparrow)$,
       dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,3)=dvxc($\downarrow$)/d$\rho(\downarrow)$

  exc(npts)=exchange-correlation energy density (hartree)
  vxc(npts,nspden)=xc potent. (d($\rho$*exc)/d($\rho\uparrow$)) and d/d($\rho\downarrow$) (ha)
  (only overall potential d($\rho$*exc)/d($\rho$) returned in vxc(1) for nspden=1)
  ndvxc= size of dvxc(npts,ndvxc)

 Normalization: Exc=$\int(exc(r)*\rho(r) d^3 r)$ for $\rho$(r)=electron density.

TODO

 To be added later
  d2vxc=derivative $d^2 (Vxc)/d(rho)^2$ (hartree*bohr^6)

NOTES

 This form is based on Mike Teter s rational polynomial
 exc=-(a0+a1*rs+a2*rs**2+a3*rs**3)/(b1*rs+b2*rs**2+b3*rs**3+b4*rs**4)
 where the parameters are fit to reproduce
 (in this case) the Perdew-Wang parameterization of the correlation
 energy given in Phys. Rev. B 45, 13244-13249 (1992) [[cite:Perdew1992]].

 Each parameter is interpolated between zeta=0 and 1 by
 a_i(zeta)=a_i(0)+(a_i(1)-a_i(0))*f_x(zeta) and
 f_x(zeta)=[(1+zeta)$^{4/3}$+(1-zeta)$^{4/3}$-2]/(2*(2$^{1/3}$-1)).

 Beware : in this expression, zeta is actually replaced by zeta*alpha_zeta,
 where alpha_zeta is very close to 1, but slightly lower.
 This is to remove the singularity in the derivatives when abs(zeta) is 1
 Below,  a_i(1)-a_i(0) is called "da" for delta a, same for b s.

 rs = $(3/(4\pi))^{1/3} * \rho(r)^{-1/3}$
 zeta = $(\rho\uparrow-\rho\downarrow)/(\rho\uparrow+\rho\downarrow)$
 b1 must be 1 and a0 must be $(3/4)(3/(2\pi))^{2/3}$.

SOURCE

236 subroutine xcspol(exc,npts,nspden,order,rspts,vxc,zeta,ndvxc,& !Mandatory arguments
237 &                 dvxc)                            !Optional arguments
238 
239 !Arguments ------------------------------------
240 !scalars
241  integer,intent(in) :: ndvxc,npts,nspden,order
242 !arrays
243  real(dp),intent(in) :: rspts(npts),zeta(npts)
244  real(dp),intent(out) :: exc(npts),vxc(npts,nspden)
245  real(dp),intent(out),optional :: dvxc(npts,ndvxc)
246 
247 !Local variables-------------------------------
248 !The generation of density from rs needs rsfac and rsfac^(-3) :
249 !rsfac=(3/(4 Pi))^(1/3) ; rsfacm3=4pi/3
250 !Mike Teter s parameters of 8 April 1993.
251 !New parameters which accomodate spin polarization (fit to P-W)
252 !Paramagnetic limit:a0p,...b4p
253 !(a0=(3/4)(3/(2 Pi))^(2/3)) (note that b1=1 is fixed)
254 !Differences, ferromagnetic - paramagnetic (delta params):da1,da2,da3,db1,db2,db3,db4
255 !scalars
256  integer :: ipts
257  real(dp),parameter :: a0p=.4581652932831429_dp,a1p=2.217058676663745_dp
258  real(dp),parameter :: a2p=0.7405551735357053_dp,a3p=0.01968227878617998_dp
259  real(dp),parameter :: alpha_zeta=one-1.0d-6,b1p=one,b2p=4.504130959426697_dp
260  real(dp),parameter :: b3p=1.110667363742916_dp,b4p=0.02359291751427506_dp
261  real(dp),parameter :: da0=.119086804055547_dp,da1=0.6157402568883345_dp
262  real(dp),parameter :: da2=0.1574201515892867_dp,da3=0.003532336663397157_dp
263  real(dp),parameter :: db1=zero,db2=0.2673612973836267_dp
264  real(dp),parameter :: db3=0.2052004607777787_dp,db4=0.004200005045691381_dp
265  real(dp),parameter :: ft=4._dp/3._dp,rsfac=0.6203504908994000_dp
266  real(dp),parameter :: rsfacm3=rsfac**(-3)
267  real(dp) :: a0,a1,a2,a3,b1,b2,b3,b4,d1,d1m1,d2d1drs2,d2d1drsdf,d2excdf2
268  real(dp) :: d2excdrs2,d2excdrsdf,d2excdz2,d2fxcdz2,d2n1drs2,d2n1drsdf,dd1df
269  real(dp) :: dd1drs,dexcdf,dexcdrs,dexcdz,dfxcdz,dn1df,dn1drs,dvxcdrs
270  real(dp) :: dvxcpdrho,dvxcpdz,excipt,fact,fxc,n1
271  real(dp) :: rhom1,rs,vxcp,zet,zetm,zetm_third
272  real(dp) :: zetp,zetp_third
273  character(len=500) :: message
274 !no_abirules
275 !Set a minimum rho below which terms are 0
276  real(dp),parameter :: rhotol=1.d-28
277 !real(dp) :: delta,rho,rho_dn,rho_dnm,rho_dnp,rho_up,rho_upm,rho_upp,zeta_mean
278 
279 ! *************************************************************************
280 
281 !Checks the compatibility between the presence of dvxc and ndvxc
282  if(ndvxc /=0 .neqv. present(dvxc))then
283    message = 'If ndvxc/=0 there must be the optional argument dvxc'
284    ABI_BUG(message)
285  end if
286 
287 !Checks the compatibility between the inputs and the presence of the optional arguments
288  if(abs(order) <= 1 .and. ndvxc /= 0)then
289    write(message, '(4a,i0)' )ch10,&
290 &   'The order chosen does not need the presence',ch10,&
291 &   'of the vector dvxc, that is needed only with |order|>1 , while we have',order
292    ABI_BUG(message)
293  end if
294 
295  if(nspden == 1 .and. ndvxc /=0 .and. ndvxc /= 2)then
296    write(message,'(a,i0)')' Once nspden=1 we must have ndvxc=2, while we have',ndvxc
297    ABI_BUG(message)
298  end if
299 
300  if(nspden == 2 .and. ndvxc /=0 .and. ndvxc /= 3)then
301    write(message, '(a,i0)' )' Once nspden=2 we must have ndvxc=3, while we have',ndvxc
302    ABI_BUG(message)
303  end if
304 
305 
306 !Although fact is parameter value, some compilers are not able to evaluate
307 !it at compile time.
308  fact=one/(two**(four*third)-two)
309 
310 !DEBUG
311 !Finite-difference debugging, do not take away
312 !debug=1
313 !zeta_mean=0.1_dp
314 !delta=0.0001
315 !if(debug==1)then
316 !do ipts=1,npts,5
317 !rho=ipts*0.01_dp
318 !rho_up=rho*(one+zeta_mean)*half
319 !rho_dn=rho*(one-zeta_mean)*half
320 !rho_upp=rho_up+delta
321 !rho_upm=rho_up-delta
322 !rho_dnp=rho_dn+delta
323 !rho_dnm=rho_dn-delta
324 !First possibility : vary rho up , and then rho down
325 !zeta(ipts  )=(rho_up -rho_dn )/(rho_up +rho_dn )
326 !zeta(ipts+1)=(rho_upp-rho_dn )/(rho_upp+rho_dn )
327 !zeta(ipts+2)=(rho_upm-rho_dn )/(rho_upm+rho_dn )
328 !zeta(ipts+3)=(rho_up -rho_dnp)/(rho_up +rho_dnp)
329 !zeta(ipts+4)=(rho_up -rho_dnm)/(rho_up +rho_dnm)
330 !rspts(ipts  )=rsfac*(rho_up +rho_dn )**(-third)
331 !rspts(ipts+1)=rsfac*(rho_upp+rho_dn )**(-third)
332 !rspts(ipts+2)=rsfac*(rho_upm+rho_dn )**(-third)
333 !rspts(ipts+3)=rsfac*(rho_up +rho_dnp)**(-third)
334 !rspts(ipts+4)=rsfac*(rho_up +rho_dnm)**(-third)
335 !DEBUGBUG : another possibility : vary rho and zeta
336 !zeta(ipts+1)=zeta(ipts  )
337 !zeta(ipts+2)=zeta(ipts  )
338 !zeta(ipts+3)=zeta(ipts  )+delta
339 !zeta(ipts+4)=zeta(ipts  )-delta
340 !rspts(ipts+1)=rsfac*(rho+delta)**(-third)
341 !rspts(ipts+2)=rsfac*(rho-delta )**(-third)
342 !rspts(ipts+3)=rspts(ipts  )
343 !rspts(ipts+4)=rspts(ipts  )
344 !ENDDEBUGBUG
345 !end do
346 !end if
347 !nspden=2
348 !order=2
349 !ENDDEBUG
350 
351  if (nspden==1) then
352 !  separate cases with respect to order
353    if(order==-2) then
354 !    No spin-polarization so skip steps related to zeta not 0
355      do ipts=1,npts
356 
357        rs=rspts(ipts)
358        n1=a0p+rs*(a1p+rs*(a2p+rs*a3p))
359        d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p)))
360        d1m1=one/d1
361 
362 !      Exchange-correlation energy
363        excipt=-n1*d1m1
364        exc(ipts)=excipt
365 
366 !      Exchange-correlation potential
367        dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p))
368        dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p)))
369 
370 !      dexcdrs is d(exc)/d(rs)
371        dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
372        vxc(ipts,1)=excipt-third*rs*dexcdrs
373 
374 !      If the exchange-correlation kernel is needed
375 
376        d2n1drs2=2._dp*a2p+rs*(6._dp*a3p)
377        d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p))
378 !      d2excdrs2 is d2(exc)/d(rs)2
379        d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1
380        dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2)
381 !      And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs)
382        dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs
383 
384 !      dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc)
385        dn1df=da0+rs*(da1+rs*(da2+rs*da3))
386        dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4)))
387        dexcdf=-(dn1df+excipt*dd1df)*d1m1
388 !      d2(fxc)/d(zeta)2
389        d2fxcdz2=ft*third*(alpha_zeta**2)*2._dp*fact
390 !      d2(exc)/d(zeta)2
391        d2excdz2=d2fxcdz2*dexcdf
392        rhom1=rsfacm3*rs**3
393        dvxc(ipts,2)= dvxc(ipts,1) - d2excdz2*rhom1
394      end do
395    else if(order**2>1) then
396 !    No spin-polarization so skip steps related to zeta not 0
397      do ipts=1,npts
398 
399        rs=rspts(ipts)
400        n1=a0p+rs*(a1p+rs*(a2p+rs*a3p))
401        d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p)))
402        d1m1=one/d1
403 
404 !      Exchange-correlation energy
405        excipt=-n1*d1m1
406        exc(ipts)=excipt
407 
408 !      Exchange-correlation potential
409        dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p))
410        dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p)))
411 
412 !      dexcdrs is d(exc)/d(rs)
413        dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
414        vxc(ipts,1)=excipt-third*rs*dexcdrs
415 
416 !      If the exchange-correlation kernel is needed
417        d2n1drs2=2._dp*a2p+rs*(6._dp*a3p)
418        d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p))
419 !      d2excdrs2 is d2(exc)/d(rs)2
420        d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1
421        dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2)
422 !      And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs)
423        dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs
424 
425      end do
426    else
427 !    No spin-polarization so skip steps related to zeta not 0
428      do ipts=1,npts
429 
430        rs=rspts(ipts)
431        n1=a0p+rs*(a1p+rs*(a2p+rs*a3p))
432        d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p)))
433        d1m1=one/d1
434 
435 !      Exchange-correlation energy
436        excipt=-n1*d1m1
437        exc(ipts)=excipt
438 
439 !      Exchange-correlation potential
440        dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p))
441        dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p)))
442 
443 !      dexcdrs is d(exc)/d(rs)
444        dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
445        vxc(ipts,1)=excipt-third*rs*dexcdrs
446      end do
447 
448    end if
449 
450 
451 !  Allows for nspden==1, in the case of testing nspden=1 against nspden=2
452  else if (nspden<=2) then
453 
454 
455 !  DEBUG
456 !  do not take away : allows to compare nspden=1 and nspden=2 coding
457 !  if (nspden==1)then
458 !  zeta(:)=zero
459 !  end if
460 !  ENDDEBUG
461 !  separate cases with respect to order
462    if(abs(order)>1) then
463 !    Allow for spin polarization. This part could be optimized for speed.
464      do ipts=1,npts
465 
466        rs=rspts(ipts)
467        zet=zeta(ipts)
468        zetp=one+zet*alpha_zeta
469        zetm=one-zet*alpha_zeta
470        zetp_third=zetp**third
471        zetm_third=zetm**third
472 !      Exchange energy spin interpolation function f(zeta)
473        fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact
474 
475        a0=a0p+fxc*da0
476        a1=a1p+fxc*da1
477        a2=a2p+fxc*da2
478        a3=a3p+fxc*da3
479        b1=b1p+fxc*db1
480        b2=b2p+fxc*db2
481        b3=b3p+fxc*db3
482        b4=b4p+fxc*db4
483 
484        n1= a0+rs*(a1+rs*(a2+rs*a3))
485        d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
486        d1m1=one/d1
487 
488 !      Exchange-correlation energy
489        excipt=-n1*d1m1
490        exc(ipts)=excipt
491 
492 !      Exchange-correlation potential
493        dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3))
494        dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4)))
495 !      dexcdrs is d(exc)/d(rs)
496        dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
497 
498 !      Only vxcp contributes when paramagnetic
499        vxcp=excipt-third*rs*dexcdrs
500 
501 !      d(fxc)/d(zeta)  (which is 0 at zeta=0)
502        dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact
503 
504 !      dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc)
505        dn1df=da0+rs*(da1+rs*(da2+rs*da3))
506        dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4)))
507 
508 !      dexcdz is d(exc)/d(zeta)
509        dexcdf=-(dn1df+excipt*dd1df)*d1m1
510        dexcdz=dfxcdz*dexcdf
511 
512 !      Compute Vxc for both spin channels
513 
514        vxc(ipts,1)=vxcp - (zet-one)*dexcdz
515        vxc(ipts,2)=vxcp - (zet+one)*dexcdz
516 
517 !      DEBUG Allow to check the variation of rho and zeta
518 !      vxc(ipts,1)=vxcp
519 !      vxc(ipts,2)=dexcdz
520 !      ENDDEBUG
521 !      Compute second derivative with respect to rho
522        d2n1drs2=2._dp*a2+rs*(6._dp*a3)
523        d2d1drs2=2._dp*b2+rs*(6._dp*b3+rs*(12._dp*b4))
524 !      d2excdrs2 is d2(exc)/d(rs)2
525        d2excdrs2=-(d2n1drs2+two*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1
526        dvxcdrs=third*(two*dexcdrs-rs*d2excdrs2)
527 !      And d(vxc)/d(rho) paramagnetic =(-rs/(3*rho))*d(vxcp)/d(rs)
528 !      remember : 1/rho=(4pi/3)*rs**3=rsfacm3*rs**3
529        rhom1=rsfacm3*rs**3
530        dvxcpdrho= -rs*rhom1*third * dvxcdrs
531 
532 !      Compute mixed second derivative with respect to rho and zeta
533        d2n1drsdf=da1+rs*(2._dp*da2+rs*(3._dp*da3))
534        d2d1drsdf=db1+rs*(2._dp*db2+rs*(3._dp*db3+rs*(4._dp*db4)))
535 !      d2excdrsdf is d2(exc)/d(rs)df
536        d2excdrsdf=-(d2n1drsdf+dexcdrs*dd1df+dexcdf*dd1drs+excipt*d2d1drsdf)*d1m1
537 !      d(vxc)/d(zeta) paramagnetic
538        dvxcpdz=dexcdz-third*rs*dfxcdz*d2excdrsdf
539 
540 !      Compute second derivative with respect to zeta
541 !      the second derivative of n1 and d1 wrt f vanishes
542        d2excdf2=-(two*dexcdf*dd1df)*d1m1
543 !      d2(fxc)/d(zeta)2
544        d2fxcdz2=ft*third*(alpha_zeta**2)*(zetp_third**(-2)+zetm_third**(-2))*fact
545 !      d2(exc)/d(zeta)2
546        d2excdz2=d2fxcdz2*dexcdf+dfxcdz**2*d2excdf2
547 
548 !      Compute now the three second derivatives of the Exc energy with respect
549 !      to : wrt twice spin-up ; wrt spin-up and spin-dn ; wrt twice spin-down
550        dvxc(ipts,1)= dvxcpdrho   &
551 &       +two*rhom1*( one-zet)*(dvxcpdz-dexcdz) &
552 &       +d2excdz2*rhom1*(one-zet)**2
553        dvxc(ipts,2)= dvxcpdrho   &
554 &       +two*rhom1*(    -zet)*(dvxcpdz-dexcdz) &
555 &       +d2excdz2*rhom1*(one-zet)*(-one-zet)
556 !      if(nspden==2)then
557        dvxc(ipts,3)= dvxcpdrho   &
558 &       +two*rhom1*(-one-zet)*(dvxcpdz-dexcdz) &
559 &       +d2excdz2*rhom1*(-one-zet)**2
560 !      else
561 !      !    For testing purposes, need the spin-averaged quantity
562 !      dvxc(ipts,1)= ( dvxc(ipts,1) + dvxc(ipts,2) ) * half
563 !      end if
564 
565 !      DEBUG Allow to check the variation of rho and zeta
566 !      dvxc(ipts,1)=dvxcpdrho
567 !      dvxc(ipts,2)=d2excdz2
568 !      dvxc(ipts,3)=dvxcpdz
569 !      ENDDEBUG
570      end do
571    else
572 !    Allow for spin polarization. This part could be optimized for speed.
573      do ipts=1,npts
574 
575        rs=rspts(ipts)
576        zet=zeta(ipts)
577        zetp=one+zet*alpha_zeta
578        zetm=one-zet*alpha_zeta
579        zetp_third=zetp**third
580        zetm_third=zetm**third
581 !      Exchange energy spin interpolation function f(zeta)
582        fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact
583 
584        a0=a0p+fxc*da0
585        a1=a1p+fxc*da1
586        a2=a2p+fxc*da2
587        a3=a3p+fxc*da3
588        b1=b1p+fxc*db1
589        b2=b2p+fxc*db2
590        b3=b3p+fxc*db3
591        b4=b4p+fxc*db4
592 
593        n1= a0+rs*(a1+rs*(a2+rs*a3))
594        d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
595        d1m1=one/d1
596 
597 !      Exchange-correlation energy
598        excipt=-n1*d1m1
599        exc(ipts)=excipt
600 
601 !      Exchange-correlation potential
602        dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3))
603        dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4)))
604 !      dexcdrs is d(exc)/d(rs)
605        dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
606 
607 !      Only vxcp contributes when paramagnetic
608        vxcp=excipt-third*rs*dexcdrs
609 
610 !      d(fxc)/d(zeta)  (which is 0 at zeta=0)
611        dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact
612 
613 !      dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc)
614        dn1df=da0+rs*(da1+rs*(da2+rs*da3))
615        dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4)))
616 
617 !      dexcdz is d(exc)/d(zeta)
618        dexcdf=-(dn1df+excipt*dd1df)*d1m1
619        dexcdz=dfxcdz*dexcdf
620 
621 !      Compute Vxc for both spin channels
622 
623        vxc(ipts,1)=vxcp - (zet-one)*dexcdz
624        vxc(ipts,2)=vxcp - (zet+one)*dexcdz
625 
626 !      DEBUG Allow to check the variation of rho and zeta
627 !      vxc(ipts,1)=vxcp
628 !      vxc(ipts,2)=dexcdz
629 !      ENDDEBUG
630      end do
631    end if
632  else
633 
634 !  Disallowed value for nspden
635    write(message, '(3a,i0)' )&
636 &   ' Argument nspden must be 1 or 2; ',ch10,&
637 &   ' Value provided as argument was ',nspden
638    ABI_BUG(message)
639  end if
640 
641 !DEBUG
642 !Finite-difference debugging, do not take away
643 !if(debug==1)then
644 !write(std_out,*)' delta =',delta
645 !do ipts=1,npts,5
646 !rho=(rspts(ipts)/rsfac)**(-3)
647 !write(std_out,'(a,i5,a,2es16.8)' ) ' Point number',ipts,' with rho,zeta=',rho,zeta(ipts)
648 !write(std_out,'(3es16.8)' )exc(ipts)*rho,vxc(ipts,1),vxc(ipts,2)
649 !write(std_out,'(3es16.8)' )dvxc(ipts,1),dvxc(ipts,3),dvxc(ipts,2)
650 !write(std_out,'(3es16.8)' )exc(ipts)*rho,&
651 !&      ( exc(ipts+1)*(rho+delta) - exc(ipts+2)*(rho-delta) )/2._dp/delta,&
652 !&      ( exc(ipts+3)*(rho+delta) - exc(ipts+4)*(rho-delta) )/2._dp/delta
653 !write(std_out,'(4es16.8)' )&
654 !&    ( vxc(ipts+1,1) - vxc(ipts+2,1) )/2._dp/delta,&
655 !&    ( vxc(ipts+3,2) - vxc(ipts+4,2) )/2._dp/delta,&
656 !&    ( vxc(ipts+3,1) - vxc(ipts+4,1) )/2._dp/delta,&
657 !&    ( vxc(ipts+1,2) - vxc(ipts+2,2) )/2._dp/delta
658 !end do
659 !stop
660 !end if
661 !ENDDEBUG
662 
663 !DEBUG
664 !if(order==-2)then
665 !write(std_out,*)' xcspol : ipts,npts ',ipts,npts
666 !write(std_out,*)dvxcdrs,d2excdz2,d2fxcdz2,dexcdf
667 !write(std_out,*)rhom1
668 !write(std_out,*)dvxc(1000,1),dvxc(1000,2)
669 !stop
670 !end if
671 !ENDDEBUG
672 
673 end subroutine xcspol

ABINIT/xctetr [ Functions ]

[ Top ] [ Functions ]

NAME

 xctetr

FUNCTION

 Returns exc, vxc, and d(vxc)/d($\rho$) from input $\rho$.
 Also returns $d^2(Vxc)/d(\rho)^2$ as needed for third-order DFPT

INPUTS

  npt=number of real space points on which density is provided
  order=gives the maximal derivative of Exc computed.
  rhor(npt)=electron number density (bohr^-3)
  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)
  if(order>2) d2vxc(npt)=derivative d$^2$(Vxc)/d$(\rho)^2$ (hartree*bohr^6)

NOTES

 Teter exchange and correlation (xc)--Mike Teter s fit
 to Ceperly-Alder electron gas energy data.  Data from
 D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980) [[cite:Ceperley1980]]
 and private communication from authors.
 This form is based on Mike Teter s rational polynomial
 exc=-(a0+a1*rs+a2*rs**2+a3*rs**3)/(b1*rs+b2*rs**2+b3*rs**3+b4*rs**4)
 where the parameters of the fit are fit to reproduce
 Ceperley-Alder data and the high density limit (rs->0)
 of the electron gas (pure exchange).
 rs = $(3/(4\pi))^{1/3} * \rho(r)^{-1/3}$.
 b1 must be 1 and a0 must be $(3/4)(3/(2\pi))^{2/3}$.
 Fit is by Mike Teter, Corning Incorporated.
 Note that d(vxc)/d($\rho$) gets a little wild at small rho.
 d$^2$(Vxc)/d$(\rho)^2$ is probably wilder.

 Some notation:  (XG 990224, sign convention should be changed, see xcspol.f)
  $Exc = N1/D1$ with $N1=-(a0+a1*rs+...)$ given above and
              $D1= (b1*rs+b2*rs^2+...)$ also given above.
  $Vxc = N2/D1^2$ with $N2=d(N1)/d(rs)$.
  $d(Vxc)/d(rs)=(N3-D3*(2*N2/D1))/D1^2 with N3=d(N2)/d(rs)$ and
              $D3=d(D1)/d(rs)$.
  $d(Vxc)/d(\rho) = (-rs/(3*\rho))* d(Vxc)/d(rs)$.
  $d^2(Vxc)/d(rs)^2 = (N4-2*(2*N3*D3+N2*D4-3*N2*D3^2/D1)/D1)/D1^2$
   with $N4=d(N3)/d(rs), D4=d(D3)/d(rs)$.
  $d^2(Vxc)/d(\rho)^2= rs/(3*\rho)^2)*(4*d(Vxc)/d(rs)+rs*d^2(Vxc)/d(rs)^2)$.

SOURCE

726 subroutine xctetr(exc,npt,order,rhor,rspts,vxc,& !Mandatory arguments
727 &                 d2vxc,dvxc)                    !Optional arguments
728 
729 !Arguments ------------------------------------
730 !scalars
731  integer,intent(in) :: npt,order
732 !arrays
733  real(dp),intent(in) :: rhor(npt),rspts(npt)
734  real(dp),intent(out) :: exc(npt),vxc(npt)
735  real(dp),intent(out),optional :: d2vxc(npt),dvxc(npt)
736 
737 !Local variables-------------------------------
738 !rsfac=(3/(4 Pi))^(1/3)
739 !Mike Teter s parameters: (keep 8 digits after decimal)
740 !(a0=(3/4)(3/(2 Pi))^(2/3)
741 !scalars
742  integer :: ipt
743  real(dp),parameter :: a0=.4581652932831429_dp,a1=2.40875407_dp,a2=.88642404_dp
744  real(dp),parameter :: a3=.02600342_dp,b1=1.0_dp,b2=4.91962865_dp
745  real(dp),parameter :: b3=1.34799453_dp,b4=.03120453_dp,c1=4._dp*a0*b1/3.0_dp
746  real(dp),parameter :: c2=5.0_dp*a0*b2/3.0_dp+a1*b1
747  real(dp),parameter :: c3=2.0_dp*a0*b3+4.0_dp*a1*b2/3.0_dp+2.0_dp*a2*b1/3.0_dp
748  real(dp),parameter :: c4=7.0_dp*a0*b4/3.0_dp+5.0_dp*a1*b3/3.0_dp+a2*b2+a3*b1/3.0_dp
749  real(dp),parameter :: c5=2.0_dp*a1*b4+4.0_dp*a2*b3/3.0_dp+2.0_dp*a3*b2/3.0_dp
750  real(dp),parameter :: c6=5.0_dp*a2*b4/3.0_dp+a3*b3,c7=4.0_dp*a3*b4/3.0_dp
751  real(dp),parameter :: rsfac=0.6203504908994000_dp
752  real(dp) :: d1,d1m1,d2vxcr,d3,d4,dvxcdr,n1,n2,n3,n4,rhom1,rs
753  character(len=500) :: message
754 
755 ! *************************************************************************
756 !
757 !Checks the values of order
758  if(order<0 .or. order>3)then
759    write(message, '(a,a,a,i6)' )&
760 &   'With Teter 91 Ceperley-Alder xc functional, the only',ch10,&
761 &   'allowed values for order are 0, 1, 2 or 3, while it is found to be',order
762    ABI_BUG(message)
763  end if
764 
765 !Checks the compatibility between the order and the presence of the optional arguments
766  if(order /=3 .and. present(d2vxc))then
767    write(message, '(a,a,a,i6)' )&
768 &   'The order chosen does not need the presence',ch10,&
769 &   'of the vector d2vxc, that is needed only with order=3, while we have',order
770    ABI_BUG(message)
771  end if
772 
773  if(order <= 1 .and. present(dvxc))then
774    write(message, '(a,a,a,i6)' )&
775 &   'The order chosen does not need the presence',ch10,&
776 &   'of the vector dvxc, that is needed with order > 1, while we have',order
777    ABI_BUG(message)
778  end if
779 
780 !separated cases with respect to order
781 
782  if (order<=1) then
783 !  Loop over grid points
784    do ipt=1,npt
785      rs=rspts(ipt)
786      n1=-(a0+rs*(a1+rs*(a2+rs*a3)))
787      d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
788      d1m1=1.0_dp/d1
789      n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7))))))
790 !
791 !    Exchange-correlation energy
792      exc(ipt)=n1*d1m1
793 !
794 !    Exchange-correlation potential
795      vxc(ipt)=n2*d1m1**2
796    end do
797  else if (order>2) then
798 !  Loop over grid points
799    do ipt=1,npt
800      rs=rspts(ipt)
801      n1=-(a0+rs*(a1+rs*(a2+rs*a3)))
802      d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
803      d1m1=1.0_dp/d1
804      n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7))))))
805 !
806 !    Exchange-correlation energy
807      exc(ipt)=n1*d1m1
808 !
809 !    Exchange-correlation potential
810      vxc(ipt)=n2*d1m1**2
811 !    Assemble derivative of vxc wrt rs
812      n3=-(c1+rs*(2._dp*c2+rs*(3._dp*c3+rs*(4._dp*c4+rs*(5._dp*c5+&
813 &     rs*(6._dp*c6+rs*(7._dp*c7)))))))
814      d3=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4)))
815      dvxcdr=(n3-d3*(2._dp*n2*d1m1))*d1m1**2
816      rhom1=1.0_dp/rhor(ipt)
817 !
818 !    derivative of vxc wrt rho
819      dvxc(ipt)=-dvxcdr*rs*third*rhom1
820 !
821 
822 !    Assemble derivative d^2(Vxc)/d(rs)^2
823      n4=-(2.0_dp*c2+rs*(6.0_dp*c3+rs*(12.0_dp*c4+rs*(20.0_dp*c5+&
824 &     rs*(30.0_dp*c6+rs*(42.0_dp*c7))))))
825      d4=2.0_dp*b2+rs*(6.0_dp*b3+rs*(12.0_dp*b4))
826      d2vxcr=(n4-2.0_dp*(2.0_dp*n3*d3+n2*d4-3.0_dp*n2*d3**2*d1m1)*d1m1)*d1m1**2
827 
828 !    Derivative d^2(Vxc)/d(rho)^2
829      d2vxc(ipt)=(rs*third*rhom1)*(4.0_dp*dvxcdr+rs*d2vxcr)*third*rhom1
830 
831    end do
832  else if (order>1) then
833 !  Loop over grid points
834    do ipt=1,npt
835      rs=rspts(ipt)
836      n1=-(a0+rs*(a1+rs*(a2+rs*a3)))
837      d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
838      d1m1=1.0_dp/d1
839      n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7))))))
840 !
841 !    Exchange-correlation energy
842      exc(ipt)=n1*d1m1
843 !
844 !    Exchange-correlation potential
845      vxc(ipt)=n2*d1m1**2
846 !    Assemble derivative of vxc wrt rs
847      n3=-(c1+rs*(2._dp*c2+rs*(3._dp*c3+rs*(4._dp*c4+rs*(5._dp*c5+&
848 &     rs*(6._dp*c6+rs*(7._dp*c7)))))))
849      d3=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4)))
850      dvxcdr=(n3-d3*(2._dp*n2*d1m1))*d1m1**2
851      rhom1=1.0_dp/rhor(ipt)
852 !
853 !    derivative of vxc wrt rho
854      dvxc(ipt)=-dvxcdr*rs*third*rhom1
855 !
856    end do
857  end if
858 end subroutine xctetr

ABINIT/xctfw [ Functions ]

[ Top ] [ Functions ]

NAME

 xctfw

FUNCTION

 Add gradient part of the Thomas-Fermi-Weizsacker functional
 Perrot F., Phys. Rev. A20, 586-594 (1979) [[cite:Perrot1979]]

INPUTS

  ndvxcdgr= size of dvxcdgr(npts,ndvxcdgr)
  npts= number of points to be computed
  nspden=number if spin density component (necessarily 1 here)
  grho2_updn(npts,2*nspden-1)=square of the gradient of the spin-up,
     and, if nspden==2, spin-down, and total density (Hartree/Bohr**2),
     only used if gradient corrected functional (option=2,-2,-4 and 4 or beyond)
  rho_updn(npts,nspden)=spin-up and spin-down density (Hartree/bohr**3)
  temp= electronic temperature
  usefxc=1 if free energy fxc is used

SIDE EFFECTS

  The following arrays are modified (gradient correction added):
  dvxcdgr(npts,3)=partial derivative of the XC energy divided by the norm of the gradient
  exci(npts)=exchange-correlation energy density
  fxci(npts)=free energy energy density
  vxci(npts,nspden)=exchange-correlation potential

SOURCE

1252 subroutine xctfw(temp,exci,fxci,usefxc,rho_updn,vxci,npts,nspden,dvxcdgr,ndvxcdgr,grho2_updn)
1253 
1254 !Arguments ------------------------------------
1255 !scalars
1256  integer,intent(in) :: ndvxcdgr,npts,nspden,usefxc
1257  real(dp),intent(in) :: temp
1258 !arrays
1259  real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden)
1260  real(dp),intent(inout) :: dvxcdgr(npts,ndvxcdgr),fxci(npts*usefxc),exci(npts),vxci(npts,nspden)
1261 
1262 !Local variables-------------------------------
1263 !scalars
1264  integer :: iperrot,ipts
1265  logical :: has_fxc,has_dvxcdgr
1266  real(dp) :: etfw,rho,rho_inv,rhomot,yperrot0,vtfw
1267  real(dp) :: yperrot,uperrot,dyperrotdn,duperrotdyperrot
1268  real(dp) :: hperrot,dhperrotdyperrot,dhperrotdn,dhperrotduperrot
1269 !arrays
1270  real(dp) :: wpy(0:7), wpu(0:7)
1271  real(dp),allocatable :: rho_updnm1_3(:,:)
1272 
1273 ! *************************************************************************
1274 
1275  has_fxc=(usefxc/=0)
1276  has_dvxcdgr=(ndvxcdgr/=0)
1277 
1278  yperrot0=1.666081101_dp
1279 
1280  wpy(0)=0.5_dp; wpy(1)=-0.1999176316_dp
1281  wpy(2)=0.09765615709_dp; wpy(3)=-0.06237609924_dp
1282  wpy(4)=0.05801466322_dp; wpy(5)=-0.04449287774_dp
1283  wpy(6)=0.01903211697_dp; wpy(7)=-0.003284096926_dp
1284 
1285  wpu(0)=one/6._dp; wpu(1)=0.311590799_dp
1286  wpu(2)=3.295662439_dp; wpu(3)=-29.22038326_dp
1287  wpu(4)=116.1084531_dp; wpu(5)=-250.4543147_dp
1288  wpu(6)=281.433688_dp; wpu(7)=-128.8784806_dp
1289 
1290  ABI_MALLOC(rho_updnm1_3,(npts,2))
1291 
1292  call invcb(rho_updn(:,1),rho_updnm1_3(:,1),npts)
1293 
1294  do ipts=1,npts
1295 
1296    rho   =rho_updn(ipts,1)
1297    rhomot=rho_updnm1_3(ipts,1)
1298    rho_inv=rhomot*rhomot*rhomot
1299 
1300    yperrot=pi*pi/sqrt2/temp**1.5*two*rho
1301    uperrot=yperrot**(2./3.)
1302 
1303    dyperrotdn=pi*pi/sqrt2/temp**1.5*2.0_dp
1304 
1305    hperrot=zero
1306    dhperrotdyperrot=zero
1307    dhperrotduperrot=zero
1308    if(yperrot<=yperrot0)then
1309      do iperrot=0,7
1310        hperrot=hperrot+wpy(iperrot)*yperrot**iperrot
1311        dhperrotdyperrot=dhperrotdyperrot+iperrot*wpy(iperrot)*yperrot**(iperrot-1)
1312      end do
1313      hperrot=one/12.0_dp*hperrot
1314      dhperrotdyperrot=one/12.0_dp*dhperrotdyperrot
1315      dhperrotdn=dhperrotdyperrot*dyperrotdn
1316    else
1317      do iperrot=0,7
1318        hperrot=hperrot+wpu(iperrot)/uperrot**(2*iperrot)
1319        dhperrotduperrot=dhperrotduperrot-2.*iperrot*wpu(iperrot)/uperrot**(2*iperrot+1)
1320      end do
1321      hperrot=one/12.0_dp*hperrot
1322      dhperrotduperrot=one/12.0_dp*dhperrotduperrot
1323      duperrotdyperrot=two/3._dp/yperrot**(1./3.)
1324      dhperrotdn=dhperrotduperrot*duperrotdyperrot*dyperrotdn
1325    end if
1326 
1327    etfw=hperrot*grho2_updn(ipts,1)*rho_inv*rho_inv
1328    vtfw=-etfw + rho/hperrot*dhperrotdn*etfw
1329 
1330    if(yperrot<=yperrot0)then
1331      exci(ipts)   = exci(ipts) + etfw + 1.5_dp*yperrot*dhperrotdyperrot*grho2_updn(ipts,1)*rho_inv*rho_inv
1332    else
1333      exci(ipts)   = exci(ipts) + etfw + uperrot*dhperrotduperrot*grho2_updn(ipts,1)*rho_inv*rho_inv
1334    end if
1335    vxci(ipts,1) = vxci(ipts,1)  + vtfw
1336    if (has_fxc) fxci(ipts)   = fxci(ipts) + etfw
1337    if (has_dvxcdgr) dvxcdgr(ipts,1)= dvxcdgr(ipts,1)+two*hperrot*rho_inv
1338 
1339  end do
1340 
1341  ABI_FREE(rho_updnm1_3)
1342 
1343 end subroutine xctfw

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}$.

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)

SOURCE

884 subroutine xcwign(exc,npt,order,rspts,vxc,& !Mandatory arguments
885 &                dvxc)                           !Optional arguments
886 
887 !Arguments ------------------------------------
888 !scalars
889  integer,intent(in) :: npt,order
890 !arrays
891  real(dp),intent(in) :: rspts(npt)
892  real(dp),intent(out) :: exc(npt),vxc(npt)
893  real(dp),intent(out),optional :: dvxc(npt)
894 
895 !Local variables-------------------------------
896 !c1 and c2 are the Wigner parameters in hartree and bohr resp.
897 !scalars
898  integer :: ipt
899  real(dp),parameter :: c1=0.44_dp,c2=7.8_dp,c4_3=4.0_dp/3.0_dp
900  real(dp),parameter :: c8_27=8.0_dp/27.0_dp
901  real(dp) :: dfac,efac,rs,rsc2m1,rsm1,vfac,vxcnum
902  character(len=500) :: message
903 
904 ! *************************************************************************
905 
906 !Checks the values of order
907  if(order<0 .or. order>2)then
908    write(message, '(a,a,a,i0)' )&
909 &   'With Wigner xc functional, the only',ch10,&
910 &   'allowed values for order are 0, 1 or 2, while it is found to be',order
911    ABI_BUG(message)
912  end if
913 
914 !Checks the compatibility between the order and the presence of the optional arguments
915  if(order <= 1 .and. present(dvxc))then
916    write(message, '(a,a,a,i3)' )&
917 &   'The order chosen does not need the presence',ch10,&
918 &   'of the vector dvxc, that is needed only with order=2 , while we have',order
919    ABI_BUG(message)
920  end if
921 
922 !Compute vfac=(3/(2*Pi))^(2/3)
923  vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
924 !Compute efac=(3/4)*vfac
925  efac=0.75_dp*vfac
926 !Compute dfac=(4*Pi/9)*vfac
927  dfac=(4.0_dp*pi/9.0_dp)*vfac
928 
929 !separate cases with respect to order
930  if (order==2) then
931 
932 !  Loop over grid points
933    do ipt=1,npt
934      rs=rspts(ipt)
935      rsm1=1.0_dp/rs
936      rsc2m1=1.0_dp/(rs+c2)
937 !    compute energy density (hartree)
938      exc(ipt)=-c1*rsc2m1-efac*rsm1
939      vxcnum=-(c4_3*rs+c2)*c1
940 !    compute potential (hartree)
941      vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1
942 !    compute d(vxc)/d(rho) (hartree*bohr^3)
943      dvxc(ipt)=-(c8_27*pi)*(c1*rs**4)*(rs+rs+c2)*rsc2m1**3-dfac*rs**2
944    end do
945  else
946 
947 !  Loop over grid points
948    do ipt=1,npt
949      rs=rspts(ipt)
950      rsm1=1.0_dp/rs
951      rsc2m1=1.0_dp/(rs+c2)
952 !    compute energy density (hartree)
953      exc(ipt)=-c1*rsc2m1-efac*rsm1
954      vxcnum=-(c4_3*rs+c2)*c1
955 !    compute potential (hartree)
956      vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1
957    end do
958 
959  end if
960 !
961 end subroutine xcwign

ABINIT/xcxalp [ Functions ]

[ Top ] [ Functions ]

NAME

 xcxalp

FUNCTION

 Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input $\rho$.
 "X$\alpha$" method is used in this subroutine:
 a single fixed value is chosen for "alpha", set below.
 Expression is exc=-alpha*efac/rs (hartree), efac below.
 rs = $(3/(4\pi))^{1/3}* \rho (r)^{-1/3}$.

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)

SOURCE

1078 subroutine xcxalp(exc,npt,order,rspts,vxc, dvxc)  ! dvxc is optional
1079 
1080 !Arguments ------------------------------------
1081 !scalars
1082  integer,intent(in) :: npt,order
1083 !arrays
1084  real(dp),intent(in) :: rspts(npt)
1085  real(dp),intent(out) :: exc(npt),vxc(npt)
1086  real(dp),intent(out),optional :: dvxc(npt)
1087 
1088 !Local variables-------------------------------
1089 !Set value of alpha in "X-alpha" method
1090 !scalars
1091  integer :: ipt
1092  real(dp),parameter :: alpha=1.0_dp
1093  real(dp) :: dfac,efac,rs,rsm1,vfac
1094  character(len=500) :: message
1095 
1096 ! *************************************************************************
1097 
1098 !Checks the values of order
1099  if(order<0 .or. order>2)then
1100    write(message, '(a,a,a,i3)' )&
1101 &   'With X-alpha xc functional, the only',ch10,&
1102 &   'allowed values for order are 0, 1 or 2, while it is found to be',order
1103    ABI_BUG(message)
1104  end if
1105 
1106 !Compute vfac=(3/(2*Pi))^(2/3)
1107  vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
1108 !Compute efac=(3/4)*vfac
1109  efac=0.75_dp*vfac
1110 !Compute dfac=(4*Pi/9)*vfac
1111  dfac=(4.0_dp*pi/9.0_dp)*vfac
1112 
1113 !separate cases with respect to order
1114  if(order==2) then
1115 !  Loop over grid points
1116    do ipt=1,npt
1117      rs=rspts(ipt)
1118      rsm1=1.0_dp/rs
1119 !    compute energy density (hartree)
1120      exc(ipt)=-alpha*efac*rsm1
1121 !    compute potential (hartree)
1122      vxc(ipt)=-alpha*vfac*rsm1
1123 !    compute d(vxc)/d(rho) (hartree*bohr^3)
1124      dvxc(ipt)=-alpha*dfac*rs**2
1125    end do
1126  else
1127 !  Loop over grid points
1128    do ipt=1,npt
1129      rs=rspts(ipt)
1130      rsm1=1.0_dp/rs
1131 !    compute energy density (hartree)
1132      exc(ipt)=-alpha*efac*rsm1
1133 !    compute potential (hartree)
1134      vxc(ipt)=-alpha*vfac*rsm1
1135    end do
1136  end if
1137 !
1138 end subroutine xcxalp