Line data Source code
1 : !!****m* ABINIT/m_xclda
2 : !! NAME
3 : !! m_xclda
4 : !!
5 : !! FUNCTION
6 : !! LDA or LDA-like XC functionals.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA,XG,GMR,LG,MF,JFD,LK,AB)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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 : use m_special_funcs, only : tildeAx
28 : use m_numeric_tools, only : invcb
29 :
30 : implicit none
31 :
32 : private
33 : !!***
34 :
35 : public :: xcpzca ! Perdew-Zunger parameterization of Ceperly-Alder electron gas energy data.
36 : public :: xcspol ! Spin-polarized exchange and correlation, parameterized by Mike Teter
37 : public :: xctetr ! Teter exchange and correlation --Mike Teter s fit
38 : public :: xcwign ! Wigner exchange and correlation.
39 : public :: xchelu ! Hedin-Lundqvist exchange and correlation
40 : public :: xcxalp ! X$\alpha$ method.
41 : public :: xclb ! GGA like part (vx_lb) of the Leeuwen-Baerends XC potential.
42 : public :: xctfw ! Thomas-Fermi-Weizsacker functional
43 : public :: xcksdt ! corrKSDT finite-temperature xc functional
44 : public :: fec_ksdt ! corrKSDT finite-temperature xc functional (exchange energies)
45 : public :: fxc_ksdt ! corrKSDT finite-temperature xc functional (exchange-correlation energies)
46 : !!***
47 :
48 : contains
49 : !!***
50 :
51 : !!****f* ABINIT/xcpzca
52 : !! NAME
53 : !! xcpzca
54 : !!
55 : !! FUNCTION
56 : !! Returns exc, vxc, and d(vxc)/d($\rho$) from input rho.
57 : !!
58 : !! NOTE
59 : !! Perdew-Zunger parameterization of Ceperly-Alder electron gas energy data.
60 : !! J. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981) [[cite:Perdew1981]]
61 : !! D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980) [[cite:Ceperley1980]]
62 : !!
63 : !! INPUTS
64 : !! npt=number of real space points on which density is provided
65 : !! order=gives the maximal derivative of Exc computed.
66 : !! rhor(npt)=electron number density (bohr^-3)
67 : !! rspts(npt)=corresponding Wigner-Seitz radii, precomputed
68 : !!
69 : !! OUTPUT
70 : !! exc(npt)=exchange-correlation energy density (hartree)
71 : !! vxc(npt)=xc potential (d($\rho$*exc)/d($\rho$)) (hartree)
72 : !! if(order>1) dvxc(npt)=derivative d(vxc)/d($\rho$) (hartree*bohr^3)
73 : !!
74 : !! SOURCE
75 :
76 88406 : subroutine xcpzca(exc,npt,order,rhor,rspts,vxc,& !Mandatory arguments
77 : & dvxc) !Optional arguments
78 :
79 : !Arguments ------------------------------------
80 : !scalars
81 : integer,intent(in) :: npt,order
82 : !arrays
83 : real(dp),intent(in) :: rhor(npt),rspts(npt)
84 : real(dp),intent(out) :: exc(npt),vxc(npt)
85 : real(dp),intent(out),optional :: dvxc(npt)
86 :
87 : !Local variables-------------------------------
88 : !Perdew-Zunger parameters a, b, b1, b2, c, d, gamma
89 : !scalars
90 : integer :: ipt
91 : real(dp),parameter :: aa=0.0311_dp,b1=1.0529_dp,b2=0.3334_dp,bb=-0.048_dp
92 : real(dp),parameter :: c4_3=4.0_dp/3.0_dp,c7_6=7.0_dp/6.0_dp,cc=0.0020_dp
93 : real(dp),parameter :: dd=-0.0116_dp,ga=-0.1423_dp
94 : real(dp) :: den,den3,dfac,efac,logrs,rs,rsm1,t1,t2,vfac
95 : character(len=500) :: message
96 :
97 : ! *************************************************************************
98 :
99 : !Compute vfac=(3/(2*Pi))^(2/3)
100 88406 : vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
101 : !Compute efac=(3/4)*vfac
102 88406 : efac=0.75_dp*vfac
103 : !Compute dfac=(4*Pi/9)*vfac
104 88406 : dfac=(4.0_dp*pi/9.0_dp)*vfac
105 :
106 : !Checks the values of order
107 88406 : if(order<0 .or. order>2)then
108 : write(message, '(a,a,a,i0)' )&
109 0 : & 'With Perdew-Zunger Ceperley-Alder xc functional, the only',ch10,&
110 0 : & 'allowed values for order are 0, 1 or 2, while it is found to be',order
111 0 : ABI_BUG(message)
112 : end if
113 :
114 : !Checks the compatibility between the order and the presence of the optional arguments
115 88406 : if(order <= 1 .and. present(dvxc))then
116 : write(message, '(a,a,a,i0)' )&
117 0 : & 'The order chosen does not need the presence',ch10,&
118 0 : & 'of the vector dvxc, that is needed only with order=2 , while we have',order
119 0 : ABI_BUG(message)
120 : end if
121 :
122 : !separate cases with respect to order
123 88406 : if(order==2) then
124 : ! Loop over grid points
125 204779910 : do ipt=1,npt
126 204728136 : rs=rspts(ipt)
127 204728136 : rsm1=1.0_dp/rs
128 : ! Consider two regimes: rs<1 or rs>=1
129 204779910 : if (rs<1._dp) then
130 4263956 : logrs=log(rs)
131 : ! compute energy density exc (hartree)
132 4263956 : exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1
133 : ! compute potential vxc=d(rho*exc)/d(rho) (hartree)
134 : vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+&
135 4263956 : & (bb-aa*third)-vfac*rsm1
136 : ! compute d(vxc)/d(rho) (hartree*bohr^3)
137 : dvxc(ipt)=-(3._dp*aa+(cc+dd+dd)*rs+2._dp*cc*rs*logrs)&
138 4263956 : & /(9._dp*rhor(ipt))-dfac*rs**2
139 200464180 : else if (rs<1000._dp) then
140 185474176 : t1=b1*sqrt(rs)
141 185474176 : t2=b2*rs
142 185474176 : den=1._dp/(1._dp+t1+t2)
143 185474176 : exc(ipt)=ga*den-efac*rsm1
144 185474176 : vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
145 185474176 : den3=den**3
146 : dvxc(ipt)=(ga*den3/(36._dp*rhor(ipt)))*(5._dp*t1+8._dp*t2+&
147 185474176 : & 7._dp*t1**2+16._dp*t2**2+21._dp*t1*t2)-dfac*rs**2
148 : else
149 14990004 : t1=b1*sqrt(rs)
150 14990004 : t2=b2*rs
151 14990004 : den=1._dp/(1._dp+t1+t2)
152 14990004 : exc(ipt)=ga*den-efac*rsm1
153 14990004 : vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
154 14990004 : dvxc(ipt)=0._dp
155 : end if
156 : end do
157 : else
158 : ! Loop over grid points
159 23542841 : do ipt=1,npt
160 23506209 : rs=rspts(ipt)
161 23506209 : rsm1=1.0_dp/rs
162 : ! Consider two regimes: rs<1 or rs>=1
163 23542841 : if (rs<1._dp) then
164 13249903 : logrs=log(rs)
165 : ! compute energy density exc (hartree)
166 13249903 : exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1
167 : ! compute potential vxc=d(rho*exc)/d(rho) (hartree)
168 : vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+&
169 13249903 : & (bb-aa*third)-vfac*rsm1
170 : ! compute d(vxc)/d(rho) (hartree*bohr^3)
171 : else
172 10256306 : t1=b1*sqrt(rs)
173 10256306 : t2=b2*rs
174 10256306 : den=1._dp/(1._dp+t1+t2)
175 10256306 : exc(ipt)=ga*den-efac*rsm1
176 10256306 : vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
177 : end if
178 : end do
179 : end if
180 : !
181 88406 : end subroutine xcpzca
182 : !!***
183 :
184 : !!****f* ABINIT/xcspol
185 : !! NAME
186 : !! xcspol
187 : !!
188 : !! FUNCTION
189 : !! Spin-polarized exchange and correlation, parameterized by Mike Teter of Corning Incorporated.
190 : !!
191 : !! INPUTS
192 : !! nspden=number of spin-density components
193 : !! npts= number of points to be computed
194 : !! order=its absolute value gives the maximal derivative of Exc to be computed.
195 : !! rspts(npts)=Seitz electron radius (bohr)
196 : !! zeta(npts)=$(\rho\uparrow-\rho\downarrow)/(\rho\uparrow+\rho\downarrow)$=degree of polarization
197 : !! (ignored if nspden=1, in which case zeta should be 0)
198 : !!
199 : !! OUTPUT
200 : !! if(abs(order)>1) dvxc(npts,1+nspden)= (Hartree*bohr^3)
201 : !! if(nspden=1 .and. order==2): dvxc(:,1)=dvxc/d$\rho$ , dvxc(:,2) empty
202 : !! if(nspden=1 .and. order==-2): also compute dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$
203 : !! if(nspden=2): dvxc(:,1)=dvxc($\uparrow$)/d$\rho(\uparrow)$,
204 : !! dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,3)=dvxc($\downarrow$)/d$\rho(\downarrow)$
205 : !!
206 : !! exc(npts)=exchange-correlation energy density (hartree)
207 : !! vxc(npts,nspden)=xc potent. (d($\rho$*exc)/d($\rho\uparrow$)) and d/d($\rho\downarrow$) (ha)
208 : !! (only overall potential d($\rho$*exc)/d($\rho$) returned in vxc(1) for nspden=1)
209 : !! ndvxc= size of dvxc(npts,ndvxc)
210 : !!
211 : !! Normalization: Exc=$\int(exc(r)*\rho(r) d^3 r)$ for $\rho$(r)=electron density.
212 : !!
213 : !! TODO
214 : !! To be added later
215 : !! d2vxc=derivative $d^2 (Vxc)/d(rho)^2$ (hartree*bohr^6)
216 : !!
217 : !! NOTES
218 : !! This form is based on Mike Teter s rational polynomial
219 : !! exc=-(a0+a1*rs+a2*rs**2+a3*rs**3)/(b1*rs+b2*rs**2+b3*rs**3+b4*rs**4)
220 : !! where the parameters are fit to reproduce
221 : !! (in this case) the Perdew-Wang parameterization of the correlation
222 : !! energy given in Phys. Rev. B 45, 13244-13249 (1992) [[cite:Perdew1992]].
223 : !!
224 : !! Each parameter is interpolated between zeta=0 and 1 by
225 : !! a_i(zeta)=a_i(0)+(a_i(1)-a_i(0))*f_x(zeta) and
226 : !! f_x(zeta)=[(1+zeta)$^{4/3}$+(1-zeta)$^{4/3}$-2]/(2*(2$^{1/3}$-1)).
227 : !!
228 : !! Beware : in this expression, zeta is actually replaced by zeta*alpha_zeta,
229 : !! where alpha_zeta is very close to 1, but slightly lower.
230 : !! This is to remove the singularity in the derivatives when abs(zeta) is 1
231 : !! Below, a_i(1)-a_i(0) is called "da" for delta a, same for b s.
232 : !!
233 : !! rs = $(3/(4\pi))^{1/3} * \rho(r)^{-1/3}$
234 : !! zeta = $(\rho\uparrow-\rho\downarrow)/(\rho\uparrow+\rho\downarrow)$
235 : !! b1 must be 1 and a0 must be $(3/4)(3/(2\pi))^{2/3}$.
236 : !!
237 : !! SOURCE
238 :
239 209860 : subroutine xcspol(exc,npts,nspden,order,rspts,vxc,zeta,ndvxc,& !Mandatory arguments
240 5660 : & dvxc) !Optional arguments
241 :
242 : !Arguments ------------------------------------
243 : !scalars
244 : integer,intent(in) :: ndvxc,npts,nspden,order
245 : !arrays
246 : real(dp),intent(in) :: rspts(npts),zeta(npts)
247 : real(dp),intent(out) :: exc(npts),vxc(npts,nspden)
248 : real(dp),intent(out),optional :: dvxc(npts,ndvxc)
249 :
250 : !Local variables-------------------------------
251 : !The generation of density from rs needs rsfac and rsfac^(-3) :
252 : !rsfac=(3/(4 Pi))^(1/3) ; rsfacm3=4pi/3
253 : !Mike Teter s parameters of 8 April 1993.
254 : !New parameters which accomodate spin polarization (fit to P-W)
255 : !Paramagnetic limit:a0p,...b4p
256 : !(a0=(3/4)(3/(2 Pi))^(2/3)) (note that b1=1 is fixed)
257 : !Differences, ferromagnetic - paramagnetic (delta params):da1,da2,da3,db1,db2,db3,db4
258 : !scalars
259 : integer :: ipts
260 : real(dp),parameter :: a0p=.4581652932831429_dp,a1p=2.217058676663745_dp
261 : real(dp),parameter :: a2p=0.7405551735357053_dp,a3p=0.01968227878617998_dp
262 : real(dp),parameter :: alpha_zeta=one-1.0d-6,b1p=one,b2p=4.504130959426697_dp
263 : real(dp),parameter :: b3p=1.110667363742916_dp,b4p=0.02359291751427506_dp
264 : real(dp),parameter :: da0=.119086804055547_dp,da1=0.6157402568883345_dp
265 : real(dp),parameter :: da2=0.1574201515892867_dp,da3=0.003532336663397157_dp
266 : real(dp),parameter :: db1=zero,db2=0.2673612973836267_dp
267 : real(dp),parameter :: db3=0.2052004607777787_dp,db4=0.004200005045691381_dp
268 : real(dp),parameter :: ft=4._dp/3._dp,rsfac=0.6203504908994000_dp
269 : real(dp),parameter :: rsfacm3=rsfac**(-3)
270 : real(dp) :: a0,a1,a2,a3,b1,b2,b3,b4,d1,d1m1,d2d1drs2,d2d1drsdf,d2excdf2
271 : real(dp) :: d2excdrs2,d2excdrsdf,d2excdz2,d2fxcdz2,d2n1drs2,d2n1drsdf,dd1df
272 : real(dp) :: dd1drs,dexcdf,dexcdrs,dexcdz,dfxcdz,dn1df,dn1drs,dvxcdrs
273 : real(dp) :: dvxcpdrho,dvxcpdz,excipt,fact,fxc,n1
274 : real(dp) :: rhom1,rs,vxcp,zet,zetm,zetm_third
275 : real(dp) :: zetp,zetp_third
276 : character(len=500) :: message
277 : !no_abirules
278 : !Set a minimum rho below which terms are 0
279 : real(dp),parameter :: rhotol=1.d-28
280 : !real(dp) :: delta,rho,rho_dn,rho_dnm,rho_dnp,rho_up,rho_upm,rho_upp,zeta_mean
281 :
282 : ! *************************************************************************
283 :
284 : !Checks the compatibility between the presence of dvxc and ndvxc
285 104930 : if(ndvxc /=0 .neqv. present(dvxc))then
286 0 : message = 'If ndvxc/=0 there must be the optional argument dvxc'
287 0 : ABI_BUG(message)
288 : end if
289 :
290 : !Checks the compatibility between the inputs and the presence of the optional arguments
291 104930 : if(abs(order) <= 1 .and. ndvxc /= 0)then
292 0 : write(message, '(4a,i0)' )ch10,&
293 0 : & 'The order chosen does not need the presence',ch10,&
294 0 : & 'of the vector dvxc, that is needed only with |order|>1 , while we have',order
295 0 : ABI_BUG(message)
296 : end if
297 :
298 104930 : if(nspden == 1 .and. ndvxc /=0 .and. ndvxc /= 2)then
299 0 : write(message,'(a,i0)')' Once nspden=1 we must have ndvxc=2, while we have',ndvxc
300 0 : ABI_BUG(message)
301 : end if
302 :
303 104930 : if(nspden == 2 .and. ndvxc /=0 .and. ndvxc /= 3)then
304 0 : write(message, '(a,i0)' )' Once nspden=2 we must have ndvxc=3, while we have',ndvxc
305 0 : ABI_BUG(message)
306 : end if
307 :
308 :
309 : !Although fact is parameter value, some compilers are not able to evaluate
310 : !it at compile time.
311 104930 : fact=one/(two**(four*third)-two)
312 :
313 : !DEBUG
314 : !Finite-difference debugging, do not take away
315 : !debug=1
316 : !zeta_mean=0.1_dp
317 : !delta=0.0001
318 : !if(debug==1)then
319 : !do ipts=1,npts,5
320 : !rho=ipts*0.01_dp
321 : !rho_up=rho*(one+zeta_mean)*half
322 : !rho_dn=rho*(one-zeta_mean)*half
323 : !rho_upp=rho_up+delta
324 : !rho_upm=rho_up-delta
325 : !rho_dnp=rho_dn+delta
326 : !rho_dnm=rho_dn-delta
327 : !First possibility : vary rho up , and then rho down
328 : !zeta(ipts )=(rho_up -rho_dn )/(rho_up +rho_dn )
329 : !zeta(ipts+1)=(rho_upp-rho_dn )/(rho_upp+rho_dn )
330 : !zeta(ipts+2)=(rho_upm-rho_dn )/(rho_upm+rho_dn )
331 : !zeta(ipts+3)=(rho_up -rho_dnp)/(rho_up +rho_dnp)
332 : !zeta(ipts+4)=(rho_up -rho_dnm)/(rho_up +rho_dnm)
333 : !rspts(ipts )=rsfac*(rho_up +rho_dn )**(-third)
334 : !rspts(ipts+1)=rsfac*(rho_upp+rho_dn )**(-third)
335 : !rspts(ipts+2)=rsfac*(rho_upm+rho_dn )**(-third)
336 : !rspts(ipts+3)=rsfac*(rho_up +rho_dnp)**(-third)
337 : !rspts(ipts+4)=rsfac*(rho_up +rho_dnm)**(-third)
338 : !DEBUGBUG : another possibility : vary rho and zeta
339 : !zeta(ipts+1)=zeta(ipts )
340 : !zeta(ipts+2)=zeta(ipts )
341 : !zeta(ipts+3)=zeta(ipts )+delta
342 : !zeta(ipts+4)=zeta(ipts )-delta
343 : !rspts(ipts+1)=rsfac*(rho+delta)**(-third)
344 : !rspts(ipts+2)=rsfac*(rho-delta )**(-third)
345 : !rspts(ipts+3)=rspts(ipts )
346 : !rspts(ipts+4)=rspts(ipts )
347 : !ENDDEBUGBUG
348 : !end do
349 : !end if
350 : !nspden=2
351 : !order=2
352 : !ENDDEBUG
353 :
354 104930 : if (nspden==1) then
355 : ! separate cases with respect to order
356 81481 : if(order==-2) then
357 : ! No spin-polarization so skip steps related to zeta not 0
358 294500 : do ipts=1,npts
359 :
360 294422 : rs=rspts(ipts)
361 294422 : n1=a0p+rs*(a1p+rs*(a2p+rs*a3p))
362 294422 : d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p)))
363 294422 : d1m1=one/d1
364 :
365 : ! Exchange-correlation energy
366 294422 : excipt=-n1*d1m1
367 294422 : exc(ipts)=excipt
368 :
369 : ! Exchange-correlation potential
370 294422 : dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p))
371 294422 : dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p)))
372 :
373 : ! dexcdrs is d(exc)/d(rs)
374 294422 : dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
375 294422 : vxc(ipts,1)=excipt-third*rs*dexcdrs
376 :
377 : ! If the exchange-correlation kernel is needed
378 :
379 294422 : d2n1drs2=2._dp*a2p+rs*(6._dp*a3p)
380 294422 : d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p))
381 : ! d2excdrs2 is d2(exc)/d(rs)2
382 294422 : d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1
383 294422 : dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2)
384 : ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs)
385 294422 : dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs
386 :
387 : ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc)
388 294422 : dn1df=da0+rs*(da1+rs*(da2+rs*da3))
389 294422 : dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4)))
390 294422 : dexcdf=-(dn1df+excipt*dd1df)*d1m1
391 : ! d2(fxc)/d(zeta)2
392 294422 : d2fxcdz2=ft*third*(alpha_zeta**2)*2._dp*fact
393 : ! d2(exc)/d(zeta)2
394 294422 : d2excdz2=d2fxcdz2*dexcdf
395 294422 : rhom1=rsfacm3*rs**3
396 294500 : dvxc(ipts,2)= dvxc(ipts,1) - d2excdz2*rhom1
397 : end do
398 81403 : else if(order**2>1) then
399 : ! No spin-polarization so skip steps related to zeta not 0
400 13086876 : do ipts=1,npts
401 :
402 13083296 : rs=rspts(ipts)
403 13083296 : n1=a0p+rs*(a1p+rs*(a2p+rs*a3p))
404 13083296 : d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p)))
405 13083296 : d1m1=one/d1
406 :
407 : ! Exchange-correlation energy
408 13083296 : excipt=-n1*d1m1
409 13083296 : exc(ipts)=excipt
410 :
411 : ! Exchange-correlation potential
412 13083296 : dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p))
413 13083296 : dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p)))
414 :
415 : ! dexcdrs is d(exc)/d(rs)
416 13083296 : dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
417 13083296 : vxc(ipts,1)=excipt-third*rs*dexcdrs
418 :
419 : ! If the exchange-correlation kernel is needed
420 13083296 : d2n1drs2=2._dp*a2p+rs*(6._dp*a3p)
421 13083296 : d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p))
422 : ! d2excdrs2 is d2(exc)/d(rs)2
423 13083296 : d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1
424 13083296 : dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2)
425 : ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs)
426 13086876 : dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs
427 :
428 : end do
429 : else
430 : ! No spin-polarization so skip steps related to zeta not 0
431 256590013 : do ipts=1,npts
432 :
433 256512190 : rs=rspts(ipts)
434 256512190 : n1=a0p+rs*(a1p+rs*(a2p+rs*a3p))
435 256512190 : d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p)))
436 256512190 : d1m1=one/d1
437 :
438 : ! Exchange-correlation energy
439 256512190 : excipt=-n1*d1m1
440 256512190 : exc(ipts)=excipt
441 :
442 : ! Exchange-correlation potential
443 256512190 : dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p))
444 256512190 : dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p)))
445 :
446 : ! dexcdrs is d(exc)/d(rs)
447 256512190 : dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
448 256590013 : vxc(ipts,1)=excipt-third*rs*dexcdrs
449 : end do
450 :
451 : end if
452 :
453 :
454 : ! Allows for nspden==1, in the case of testing nspden=1 against nspden=2
455 23449 : else if (nspden<=2) then
456 :
457 :
458 : ! DEBUG
459 : ! do not take away : allows to compare nspden=1 and nspden=2 coding
460 : ! if (nspden==1)then
461 : ! zeta(:)=zero
462 : ! end if
463 : ! ENDDEBUG
464 : ! separate cases with respect to order
465 23449 : if(abs(order)>1) then
466 : ! Allow for spin polarization. This part could be optimized for speed.
467 7787370 : do ipts=1,npts
468 :
469 7785368 : rs=rspts(ipts)
470 7785368 : zet=zeta(ipts)
471 7785368 : zetp=one+zet*alpha_zeta
472 7785368 : zetm=one-zet*alpha_zeta
473 7785368 : zetp_third=zetp**third
474 7785368 : zetm_third=zetm**third
475 : ! Exchange energy spin interpolation function f(zeta)
476 7785368 : fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact
477 :
478 7785368 : a0=a0p+fxc*da0
479 7785368 : a1=a1p+fxc*da1
480 7785368 : a2=a2p+fxc*da2
481 7785368 : a3=a3p+fxc*da3
482 7785368 : b1=b1p+fxc*db1
483 7785368 : b2=b2p+fxc*db2
484 7785368 : b3=b3p+fxc*db3
485 7785368 : b4=b4p+fxc*db4
486 :
487 7785368 : n1= a0+rs*(a1+rs*(a2+rs*a3))
488 7785368 : d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
489 7785368 : d1m1=one/d1
490 :
491 : ! Exchange-correlation energy
492 7785368 : excipt=-n1*d1m1
493 7785368 : exc(ipts)=excipt
494 :
495 : ! Exchange-correlation potential
496 7785368 : dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3))
497 7785368 : dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4)))
498 : ! dexcdrs is d(exc)/d(rs)
499 7785368 : dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
500 :
501 : ! Only vxcp contributes when paramagnetic
502 7785368 : vxcp=excipt-third*rs*dexcdrs
503 :
504 : ! d(fxc)/d(zeta) (which is 0 at zeta=0)
505 7785368 : dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact
506 :
507 : ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc)
508 7785368 : dn1df=da0+rs*(da1+rs*(da2+rs*da3))
509 7785368 : dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4)))
510 :
511 : ! dexcdz is d(exc)/d(zeta)
512 7785368 : dexcdf=-(dn1df+excipt*dd1df)*d1m1
513 7785368 : dexcdz=dfxcdz*dexcdf
514 :
515 : ! Compute Vxc for both spin channels
516 :
517 7785368 : vxc(ipts,1)=vxcp - (zet-one)*dexcdz
518 7785368 : vxc(ipts,2)=vxcp - (zet+one)*dexcdz
519 :
520 : ! DEBUG Allow to check the variation of rho and zeta
521 : ! vxc(ipts,1)=vxcp
522 : ! vxc(ipts,2)=dexcdz
523 : ! ENDDEBUG
524 : ! Compute second derivative with respect to rho
525 7785368 : d2n1drs2=2._dp*a2+rs*(6._dp*a3)
526 7785368 : d2d1drs2=2._dp*b2+rs*(6._dp*b3+rs*(12._dp*b4))
527 : ! d2excdrs2 is d2(exc)/d(rs)2
528 7785368 : d2excdrs2=-(d2n1drs2+two*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1
529 7785368 : dvxcdrs=third*(two*dexcdrs-rs*d2excdrs2)
530 : ! And d(vxc)/d(rho) paramagnetic =(-rs/(3*rho))*d(vxcp)/d(rs)
531 : ! remember : 1/rho=(4pi/3)*rs**3=rsfacm3*rs**3
532 7785368 : rhom1=rsfacm3*rs**3
533 7785368 : dvxcpdrho= -rs*rhom1*third * dvxcdrs
534 :
535 : ! Compute mixed second derivative with respect to rho and zeta
536 7785368 : d2n1drsdf=da1+rs*(2._dp*da2+rs*(3._dp*da3))
537 7785368 : d2d1drsdf=db1+rs*(2._dp*db2+rs*(3._dp*db3+rs*(4._dp*db4)))
538 : ! d2excdrsdf is d2(exc)/d(rs)df
539 7785368 : d2excdrsdf=-(d2n1drsdf+dexcdrs*dd1df+dexcdf*dd1drs+excipt*d2d1drsdf)*d1m1
540 : ! d(vxc)/d(zeta) paramagnetic
541 7785368 : dvxcpdz=dexcdz-third*rs*dfxcdz*d2excdrsdf
542 :
543 : ! Compute second derivative with respect to zeta
544 : ! the second derivative of n1 and d1 wrt f vanishes
545 7785368 : d2excdf2=-(two*dexcdf*dd1df)*d1m1
546 : ! d2(fxc)/d(zeta)2
547 7785368 : d2fxcdz2=ft*third*(alpha_zeta**2)*(zetp_third**(-2)+zetm_third**(-2))*fact
548 : ! d2(exc)/d(zeta)2
549 7785368 : d2excdz2=d2fxcdz2*dexcdf+dfxcdz**2*d2excdf2
550 :
551 : ! Compute now the three second derivatives of the Exc energy with respect
552 : ! to : wrt twice spin-up ; wrt spin-up and spin-dn ; wrt twice spin-down
553 : dvxc(ipts,1)= dvxcpdrho &
554 : & +two*rhom1*( one-zet)*(dvxcpdz-dexcdz) &
555 7785368 : & +d2excdz2*rhom1*(one-zet)**2
556 : dvxc(ipts,2)= dvxcpdrho &
557 : & +two*rhom1*( -zet)*(dvxcpdz-dexcdz) &
558 7785368 : & +d2excdz2*rhom1*(one-zet)*(-one-zet)
559 : ! if(nspden==2)then
560 : dvxc(ipts,3)= dvxcpdrho &
561 : & +two*rhom1*(-one-zet)*(dvxcpdz-dexcdz) &
562 7787370 : & +d2excdz2*rhom1*(-one-zet)**2
563 : ! else
564 : ! ! For testing purposes, need the spin-averaged quantity
565 : ! dvxc(ipts,1)= ( dvxc(ipts,1) + dvxc(ipts,2) ) * half
566 : ! end if
567 :
568 : ! DEBUG Allow to check the variation of rho and zeta
569 : ! dvxc(ipts,1)=dvxcpdrho
570 : ! dvxc(ipts,2)=d2excdz2
571 : ! dvxc(ipts,3)=dvxcpdz
572 : ! ENDDEBUG
573 : end do
574 : else
575 : ! Allow for spin polarization. This part could be optimized for speed.
576 71748438 : do ipts=1,npts
577 :
578 71726991 : rs=rspts(ipts)
579 71726991 : zet=zeta(ipts)
580 71726991 : zetp=one+zet*alpha_zeta
581 71726991 : zetm=one-zet*alpha_zeta
582 71726991 : zetp_third=zetp**third
583 71726991 : zetm_third=zetm**third
584 : ! Exchange energy spin interpolation function f(zeta)
585 71726991 : fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact
586 :
587 71726991 : a0=a0p+fxc*da0
588 71726991 : a1=a1p+fxc*da1
589 71726991 : a2=a2p+fxc*da2
590 71726991 : a3=a3p+fxc*da3
591 71726991 : b1=b1p+fxc*db1
592 71726991 : b2=b2p+fxc*db2
593 71726991 : b3=b3p+fxc*db3
594 71726991 : b4=b4p+fxc*db4
595 :
596 71726991 : n1= a0+rs*(a1+rs*(a2+rs*a3))
597 71726991 : d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
598 71726991 : d1m1=one/d1
599 :
600 : ! Exchange-correlation energy
601 71726991 : excipt=-n1*d1m1
602 71726991 : exc(ipts)=excipt
603 :
604 : ! Exchange-correlation potential
605 71726991 : dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3))
606 71726991 : dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4)))
607 : ! dexcdrs is d(exc)/d(rs)
608 71726991 : dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
609 :
610 : ! Only vxcp contributes when paramagnetic
611 71726991 : vxcp=excipt-third*rs*dexcdrs
612 :
613 : ! d(fxc)/d(zeta) (which is 0 at zeta=0)
614 71726991 : dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact
615 :
616 : ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc)
617 71726991 : dn1df=da0+rs*(da1+rs*(da2+rs*da3))
618 71726991 : dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4)))
619 :
620 : ! dexcdz is d(exc)/d(zeta)
621 71726991 : dexcdf=-(dn1df+excipt*dd1df)*d1m1
622 71726991 : dexcdz=dfxcdz*dexcdf
623 :
624 : ! Compute Vxc for both spin channels
625 :
626 71726991 : vxc(ipts,1)=vxcp - (zet-one)*dexcdz
627 71748438 : vxc(ipts,2)=vxcp - (zet+one)*dexcdz
628 :
629 : ! DEBUG Allow to check the variation of rho and zeta
630 : ! vxc(ipts,1)=vxcp
631 : ! vxc(ipts,2)=dexcdz
632 : ! ENDDEBUG
633 : end do
634 : end if
635 : else
636 :
637 : ! Disallowed value for nspden
638 : write(message, '(3a,i0)' )&
639 0 : & ' Argument nspden must be 1 or 2; ',ch10,&
640 0 : & ' Value provided as argument was ',nspden
641 0 : ABI_BUG(message)
642 : end if
643 :
644 : !DEBUG
645 : !Finite-difference debugging, do not take away
646 : !if(debug==1)then
647 : !write(std_out,*)' delta =',delta
648 : !do ipts=1,npts,5
649 : !rho=(rspts(ipts)/rsfac)**(-3)
650 : !write(std_out,'(a,i5,a,2es16.8)' ) ' Point number',ipts,' with rho,zeta=',rho,zeta(ipts)
651 : !write(std_out,'(3es16.8)' )exc(ipts)*rho,vxc(ipts,1),vxc(ipts,2)
652 : !write(std_out,'(3es16.8)' )dvxc(ipts,1),dvxc(ipts,3),dvxc(ipts,2)
653 : !write(std_out,'(3es16.8)' )exc(ipts)*rho,&
654 : !& ( exc(ipts+1)*(rho+delta) - exc(ipts+2)*(rho-delta) )/2._dp/delta,&
655 : !& ( exc(ipts+3)*(rho+delta) - exc(ipts+4)*(rho-delta) )/2._dp/delta
656 : !write(std_out,'(4es16.8)' )&
657 : !& ( vxc(ipts+1,1) - vxc(ipts+2,1) )/2._dp/delta,&
658 : !& ( vxc(ipts+3,2) - vxc(ipts+4,2) )/2._dp/delta,&
659 : !& ( vxc(ipts+3,1) - vxc(ipts+4,1) )/2._dp/delta,&
660 : !& ( vxc(ipts+1,2) - vxc(ipts+2,2) )/2._dp/delta
661 : !end do
662 : !stop
663 : !end if
664 : !ENDDEBUG
665 :
666 : !DEBUG
667 : !if(order==-2)then
668 : !write(std_out,*)' xcspol : ipts,npts ',ipts,npts
669 : !write(std_out,*)dvxcdrs,d2excdz2,d2fxcdz2,dexcdf
670 : !write(std_out,*)rhom1
671 : !write(std_out,*)dvxc(1000,1),dvxc(1000,2)
672 : !stop
673 : !end if
674 : !ENDDEBUG
675 :
676 104930 : end subroutine xcspol
677 : !!***
678 :
679 :
680 : !!****f* ABINIT/xctetr
681 : !! NAME
682 : !! xctetr
683 : !!
684 : !! FUNCTION
685 : !! Returns exc, vxc, and d(vxc)/d($\rho$) from input $\rho$.
686 : !! Also returns $d^2(Vxc)/d(\rho)^2$ as needed for third-order DFPT
687 : !!
688 : !! INPUTS
689 : !! npt=number of real space points on which density is provided
690 : !! order=gives the maximal derivative of Exc computed.
691 : !! rhor(npt)=electron number density (bohr^-3)
692 : !! rspts(npt)=corresponding Wigner-Seitz radii, precomputed
693 : !!
694 : !! OUTPUT
695 : !! exc(npt)=exchange-correlation energy density (hartree)
696 : !! vxc(npt)=xc potential (d(rho*exc)/d(rho)) (hartree)
697 : !! if(order>1) dvxc(npt)=derivative d(vxc)/d($\rho$) (hartree*bohr^3)
698 : !! if(order>2) d2vxc(npt)=derivative d$^2$(Vxc)/d$(\rho)^2$ (hartree*bohr^6)
699 : !!
700 : !! NOTES
701 : !! Teter exchange and correlation (xc)--Mike Teter s fit
702 : !! to Ceperly-Alder electron gas energy data. Data from
703 : !! D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980) [[cite:Ceperley1980]]
704 : !! and private communication from authors.
705 : !! This form is based on Mike Teter s rational polynomial
706 : !! exc=-(a0+a1*rs+a2*rs**2+a3*rs**3)/(b1*rs+b2*rs**2+b3*rs**3+b4*rs**4)
707 : !! where the parameters of the fit are fit to reproduce
708 : !! Ceperley-Alder data and the high density limit (rs->0)
709 : !! of the electron gas (pure exchange).
710 : !! rs = $(3/(4\pi))^{1/3} * \rho(r)^{-1/3}$.
711 : !! b1 must be 1 and a0 must be $(3/4)(3/(2\pi))^{2/3}$.
712 : !! Fit is by Mike Teter, Corning Incorporated.
713 : !! Note that d(vxc)/d($\rho$) gets a little wild at small rho.
714 : !! d$^2$(Vxc)/d$(\rho)^2$ is probably wilder.
715 : !!
716 : !! Some notation: (XG 990224, sign convention should be changed, see xcspol.f)
717 : !! $Exc = N1/D1$ with $N1=-(a0+a1*rs+...)$ given above and
718 : !! $D1= (b1*rs+b2*rs^2+...)$ also given above.
719 : !! $Vxc = N2/D1^2$ with $N2=d(N1)/d(rs)$.
720 : !! $d(Vxc)/d(rs)=(N3-D3*(2*N2/D1))/D1^2 with N3=d(N2)/d(rs)$ and
721 : !! $D3=d(D1)/d(rs)$.
722 : !! $d(Vxc)/d(\rho) = (-rs/(3*\rho))* d(Vxc)/d(rs)$.
723 : !! $d^2(Vxc)/d(rs)^2 = (N4-2*(2*N3*D3+N2*D4-3*N2*D3^2/D1)/D1)/D1^2$
724 : !! with $N4=d(N3)/d(rs), D4=d(D3)/d(rs)$.
725 : !! $d^2(Vxc)/d(\rho)^2= rs/(3*\rho)^2)*(4*d(Vxc)/d(rs)+rs*d^2(Vxc)/d(rs)^2)$.
726 : !!
727 : !! SOURCE
728 :
729 5044 : subroutine xctetr(exc,npt,order,rhor,rspts,vxc,& !Mandatory arguments
730 : & d2vxc,dvxc) !Optional arguments
731 :
732 : !Arguments ------------------------------------
733 : !scalars
734 : integer,intent(in) :: npt,order
735 : !arrays
736 : real(dp),intent(in) :: rhor(npt),rspts(npt)
737 : real(dp),intent(out) :: exc(npt),vxc(npt)
738 : real(dp),intent(out),optional :: d2vxc(npt),dvxc(npt)
739 :
740 : !Local variables-------------------------------
741 : !rsfac=(3/(4 Pi))^(1/3)
742 : !Mike Teter s parameters: (keep 8 digits after decimal)
743 : !(a0=(3/4)(3/(2 Pi))^(2/3)
744 : !scalars
745 : integer :: ipt
746 : real(dp),parameter :: a0=.4581652932831429_dp,a1=2.40875407_dp,a2=.88642404_dp
747 : real(dp),parameter :: a3=.02600342_dp,b1=1.0_dp,b2=4.91962865_dp
748 : real(dp),parameter :: b3=1.34799453_dp,b4=.03120453_dp,c1=4._dp*a0*b1/3.0_dp
749 : real(dp),parameter :: c2=5.0_dp*a0*b2/3.0_dp+a1*b1
750 : 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
751 : 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
752 : 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
753 : real(dp),parameter :: c6=5.0_dp*a2*b4/3.0_dp+a3*b3,c7=4.0_dp*a3*b4/3.0_dp
754 : real(dp),parameter :: rsfac=0.6203504908994000_dp
755 : real(dp) :: d1,d1m1,d2vxcr,d3,d4,dvxcdr,n1,n2,n3,n4,rhom1,rs
756 : character(len=500) :: message
757 :
758 : ! *************************************************************************
759 : !
760 : !Checks the values of order
761 5044 : if(order<0 .or. order>3)then
762 : write(message, '(a,a,a,i6)' )&
763 0 : & 'With Teter 91 Ceperley-Alder xc functional, the only',ch10,&
764 0 : & 'allowed values for order are 0, 1, 2 or 3, while it is found to be',order
765 0 : ABI_BUG(message)
766 : end if
767 :
768 : !Checks the compatibility between the order and the presence of the optional arguments
769 5044 : if(order /=3 .and. present(d2vxc))then
770 : write(message, '(a,a,a,i6)' )&
771 0 : & 'The order chosen does not need the presence',ch10,&
772 0 : & 'of the vector d2vxc, that is needed only with order=3, while we have',order
773 0 : ABI_BUG(message)
774 : end if
775 :
776 5044 : if(order <= 1 .and. present(dvxc))then
777 : write(message, '(a,a,a,i6)' )&
778 0 : & 'The order chosen does not need the presence',ch10,&
779 0 : & 'of the vector dvxc, that is needed with order > 1, while we have',order
780 0 : ABI_BUG(message)
781 : end if
782 :
783 : !separated cases with respect to order
784 :
785 5044 : if (order<=1) then
786 : ! Loop over grid points
787 15013460 : do ipt=1,npt
788 15008579 : rs=rspts(ipt)
789 15008579 : n1=-(a0+rs*(a1+rs*(a2+rs*a3)))
790 15008579 : d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
791 15008579 : d1m1=1.0_dp/d1
792 15008579 : n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7))))))
793 : !
794 : ! Exchange-correlation energy
795 15008579 : exc(ipt)=n1*d1m1
796 : !
797 : ! Exchange-correlation potential
798 15013460 : vxc(ipt)=n2*d1m1**2
799 : end do
800 163 : else if (order>2) then
801 : ! Loop over grid points
802 64066 : do ipt=1,npt
803 64040 : rs=rspts(ipt)
804 64040 : n1=-(a0+rs*(a1+rs*(a2+rs*a3)))
805 64040 : d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
806 64040 : d1m1=1.0_dp/d1
807 64040 : n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7))))))
808 : !
809 : ! Exchange-correlation energy
810 64040 : exc(ipt)=n1*d1m1
811 : !
812 : ! Exchange-correlation potential
813 64040 : vxc(ipt)=n2*d1m1**2
814 : ! Assemble derivative of vxc wrt rs
815 : n3=-(c1+rs*(2._dp*c2+rs*(3._dp*c3+rs*(4._dp*c4+rs*(5._dp*c5+&
816 64040 : & rs*(6._dp*c6+rs*(7._dp*c7)))))))
817 64040 : d3=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4)))
818 64040 : dvxcdr=(n3-d3*(2._dp*n2*d1m1))*d1m1**2
819 64040 : rhom1=1.0_dp/rhor(ipt)
820 : !
821 : ! derivative of vxc wrt rho
822 64040 : dvxc(ipt)=-dvxcdr*rs*third*rhom1
823 : !
824 :
825 : ! Assemble derivative d^2(Vxc)/d(rs)^2
826 : n4=-(2.0_dp*c2+rs*(6.0_dp*c3+rs*(12.0_dp*c4+rs*(20.0_dp*c5+&
827 64040 : & rs*(30.0_dp*c6+rs*(42.0_dp*c7))))))
828 64040 : d4=2.0_dp*b2+rs*(6.0_dp*b3+rs*(12.0_dp*b4))
829 64040 : d2vxcr=(n4-2.0_dp*(2.0_dp*n3*d3+n2*d4-3.0_dp*n2*d3**2*d1m1)*d1m1)*d1m1**2
830 :
831 : ! Derivative d^2(Vxc)/d(rho)^2
832 64066 : d2vxc(ipt)=(rs*third*rhom1)*(4.0_dp*dvxcdr+rs*d2vxcr)*third*rhom1
833 :
834 : end do
835 : else if (order>1) then
836 : ! Loop over grid points
837 241851 : do ipt=1,npt
838 241714 : rs=rspts(ipt)
839 241714 : n1=-(a0+rs*(a1+rs*(a2+rs*a3)))
840 241714 : d1=rs*(b1+rs*(b2+rs*(b3+rs*b4)))
841 241714 : d1m1=1.0_dp/d1
842 241714 : n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7))))))
843 : !
844 : ! Exchange-correlation energy
845 241714 : exc(ipt)=n1*d1m1
846 : !
847 : ! Exchange-correlation potential
848 241714 : vxc(ipt)=n2*d1m1**2
849 : ! Assemble derivative of vxc wrt rs
850 : n3=-(c1+rs*(2._dp*c2+rs*(3._dp*c3+rs*(4._dp*c4+rs*(5._dp*c5+&
851 241714 : & rs*(6._dp*c6+rs*(7._dp*c7)))))))
852 241714 : d3=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4)))
853 241714 : dvxcdr=(n3-d3*(2._dp*n2*d1m1))*d1m1**2
854 241714 : rhom1=1.0_dp/rhor(ipt)
855 : !
856 : ! derivative of vxc wrt rho
857 241851 : dvxc(ipt)=-dvxcdr*rs*third*rhom1
858 : !
859 : end do
860 : end if
861 5044 : end subroutine xctetr
862 : !!***
863 :
864 : !!****f* ABINIT/xcwign
865 : !! NAME
866 : !! xcwign
867 : !!
868 : !! FUNCTION
869 : !! Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input $\rho$.
870 : !! Wigner exchange and correlation (xc)--see e.g. David Pines,
871 : !! Elementary Excitations in Solids, p. 94, NY 1964.
872 : !! Expression is exc=-(0.44)/(rs+7.8)-efac/rs (hartree), efac below.
873 : !! rs = $(3/(4\pi))^{1/3}* \rho (r)^{-1/3}$.
874 : !!
875 : !! INPUTS
876 : !! npt=number of real space points on which density is provided
877 : !! order=gives the maximal derivative of Exc computed.
878 : !! rspts(npt)=corresponding Wigner-Seitz radii, precomputed
879 : !!
880 : !! OUTPUT
881 : !! exc(npt)=exchange-correlation energy density (hartree)
882 : !! vxc(npt)=xc potential (d($\rho$*exc)/d($\rho$)) (hartree)
883 : !! if(order>1) dvxc(npt)=derivative d(vxc)/d($\rho$) (hartree*bohr^3)
884 : !!
885 : !! SOURCE
886 :
887 396 : subroutine xcwign(exc,npt,order,rspts,vxc,& !Mandatory arguments
888 : & dvxc) !Optional arguments
889 :
890 : !Arguments ------------------------------------
891 : !scalars
892 : integer,intent(in) :: npt,order
893 : !arrays
894 : real(dp),intent(in) :: rspts(npt)
895 : real(dp),intent(out) :: exc(npt),vxc(npt)
896 : real(dp),intent(out),optional :: dvxc(npt)
897 :
898 : !Local variables-------------------------------
899 : !c1 and c2 are the Wigner parameters in hartree and bohr resp.
900 : !scalars
901 : integer :: ipt
902 : real(dp),parameter :: c1=0.44_dp,c2=7.8_dp,c4_3=4.0_dp/3.0_dp
903 : real(dp),parameter :: c8_27=8.0_dp/27.0_dp
904 : real(dp) :: dfac,efac,rs,rsc2m1,rsm1,vfac,vxcnum
905 : character(len=500) :: message
906 :
907 : ! *************************************************************************
908 :
909 : !Checks the values of order
910 396 : if(order<0 .or. order>2)then
911 : write(message, '(a,a,a,i0)' )&
912 0 : & 'With Wigner xc functional, the only',ch10,&
913 0 : & 'allowed values for order are 0, 1 or 2, while it is found to be',order
914 0 : ABI_BUG(message)
915 : end if
916 :
917 : !Checks the compatibility between the order and the presence of the optional arguments
918 396 : if(order <= 1 .and. present(dvxc))then
919 : write(message, '(a,a,a,i3)' )&
920 0 : & 'The order chosen does not need the presence',ch10,&
921 0 : & 'of the vector dvxc, that is needed only with order=2 , while we have',order
922 0 : ABI_BUG(message)
923 : end if
924 :
925 : !Compute vfac=(3/(2*Pi))^(2/3)
926 396 : vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
927 : !Compute efac=(3/4)*vfac
928 396 : efac=0.75_dp*vfac
929 : !Compute dfac=(4*Pi/9)*vfac
930 396 : dfac=(4.0_dp*pi/9.0_dp)*vfac
931 :
932 : !separate cases with respect to order
933 396 : if (order==2) then
934 :
935 : ! Loop over grid points
936 0 : do ipt=1,npt
937 0 : rs=rspts(ipt)
938 0 : rsm1=1.0_dp/rs
939 0 : rsc2m1=1.0_dp/(rs+c2)
940 : ! compute energy density (hartree)
941 0 : exc(ipt)=-c1*rsc2m1-efac*rsm1
942 0 : vxcnum=-(c4_3*rs+c2)*c1
943 : ! compute potential (hartree)
944 0 : vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1
945 : ! compute d(vxc)/d(rho) (hartree*bohr^3)
946 0 : dvxc(ipt)=-(c8_27*pi)*(c1*rs**4)*(rs+rs+c2)*rsc2m1**3-dfac*rs**2
947 : end do
948 : else
949 :
950 : ! Loop over grid points
951 1573260 : do ipt=1,npt
952 1572864 : rs=rspts(ipt)
953 1572864 : rsm1=1.0_dp/rs
954 1572864 : rsc2m1=1.0_dp/(rs+c2)
955 : ! compute energy density (hartree)
956 1572864 : exc(ipt)=-c1*rsc2m1-efac*rsm1
957 1572864 : vxcnum=-(c4_3*rs+c2)*c1
958 : ! compute potential (hartree)
959 1573260 : vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1
960 : end do
961 :
962 : end if
963 : !
964 396 : end subroutine xcwign
965 : !!***
966 :
967 :
968 : !!****f* ABINIT/xchelu
969 : !! NAME
970 : !! xchelu
971 : !!
972 : !! FUNCTION
973 : !! Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input rho.
974 : !!
975 : !! NOTES
976 : !! Hedin-Lundqvist exchange and correlation (xc)--
977 : !! L. Hedin and B.I. Lundqvist, J. Phys. C. 4, 2064 (1971) [[cite:Hedin1971]]
978 : !!
979 : !! INPUTS
980 : !! npt=number of real space points on which density is provided
981 : !! order=gives the maximal derivative of Exc computed.
982 : !! rspts(npt)=Wigner-Seitz radii at each point
983 : !!
984 : !! OUTPUT
985 : !! exc(npt)=exchange-correlation energy density (hartree)
986 : !! vxc(npt)=xc potential (d($\rho$*exc)/d($\rho$)) (hartree)
987 : !! if(order>1) dvxc(npt)=derivative d(vxc)/d($\rho$) (hartree*bohr^3)
988 : !!
989 : !! SOURCE
990 :
991 396 : subroutine xchelu(exc,npt,order,rspts,vxc,dvxc) ! dvxc is optional
992 :
993 : !Arguments ------------------------------------
994 : !scalars
995 : integer,intent(in) :: npt,order
996 : !arrays
997 : real(dp),intent(in) :: rspts(npt)
998 : real(dp),intent(out) :: exc(npt),vxc(npt)
999 : real(dp),intent(out),optional :: dvxc(npt)
1000 :
1001 : !Local variables-------------------------------
1002 : !aa and cc are H-L fitting parameters A and C (C in hartree)
1003 : !rs = (3/(4 Pi))**(1/3) * rho(r)**(-1/3).
1004 : !scalars
1005 : integer :: ipt
1006 : real(dp),parameter :: aa=21_dp,c1_21=one/21_dp,c4_9=4.0_dp/9.0_dp,cc=0.0225_dp
1007 : real(dp) :: dfac,efac,rs,rsm1,vfac,xx
1008 : character(len=500) :: message
1009 :
1010 : ! *************************************************************************
1011 :
1012 : !Checks the values of order
1013 396 : if(order<0 .or. order>2)then
1014 : write(message, '(a,a,a,i0)' )&
1015 0 : & 'With Hedin-Lundqvist xc functional, the only',ch10,&
1016 0 : & 'allowed values for order are 0, 1 or 2, while it is found to be',order
1017 0 : ABI_BUG(message)
1018 : end if
1019 :
1020 : !Compute vfac=(3/(2*Pi))^(2/3)
1021 396 : vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
1022 : !Compute efac=(3/4)*vfac
1023 396 : efac=0.75_dp*vfac
1024 : !Compute dfac=(4*Pi/9)*vfac
1025 396 : dfac=(4.0_dp*pi/9.0_dp)*vfac
1026 : !separate cases with respect to order
1027 396 : if (order==2) then
1028 : ! Loop over grid points
1029 0 : do ipt=1,npt
1030 0 : rs=rspts(ipt)
1031 0 : rsm1=one/rs
1032 : ! compute energy density exc (hartree)
1033 0 : xx=rs*c1_21
1034 : exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+&
1035 0 : & half*xx-xx*xx-third) - efac*rsm1
1036 : ! compute xc potential d(rho*exc)/d(rho) (hartree)
1037 0 : vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1
1038 : ! compute d(vxc)/d(rho) (hartree*bohr^3)
1039 0 : dvxc(ipt)=-(rs**2)*((c4_9*pi)*cc*rs/(one+xx) + dfac)
1040 : end do
1041 : else
1042 : ! Loop over grid points
1043 1573260 : do ipt=1,npt
1044 1572864 : rs=rspts(ipt)
1045 1572864 : rsm1=one/rs
1046 : ! compute energy density exc (hartree)
1047 1572864 : xx=rs*c1_21
1048 : exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+&
1049 1572864 : & half*xx-xx*xx-third) - efac*rsm1
1050 : ! compute xc potential d(rho*exc)/d(rho) (hartree)
1051 1573260 : vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1
1052 : end do
1053 : end if
1054 : !
1055 396 : end subroutine xchelu
1056 : !!***
1057 :
1058 : !!****f* ABINIT/xcxalp
1059 : !! NAME
1060 : !! xcxalp
1061 : !!
1062 : !! FUNCTION
1063 : !! Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input $\rho$.
1064 : !! "X$\alpha$" method is used in this subroutine:
1065 : !! a single fixed value is chosen for "alpha", set below.
1066 : !! Expression is exc=-alpha*efac/rs (hartree), efac below.
1067 : !! rs = $(3/(4\pi))^{1/3}* \rho (r)^{-1/3}$.
1068 : !!
1069 : !! INPUTS
1070 : !! npt=number of real space points on which density is provided
1071 : !! order=gives the maximal derivative of Exc computed.
1072 : !! rspts(npt)=Wigner-Seitz radii, at each point
1073 : !!
1074 : !! OUTPUT
1075 : !! exc(npt)=exchange-correlation energy density (hartree)
1076 : !! vxc(npt)=xc potential (d($\rho$*exc)/d($\rho$)) (hartree)
1077 : !! if(order>1) dvxc(npt)=derivative d(vxc)/d($\rho$) (hartree*bohr^3)
1078 : !!
1079 : !! SOURCE
1080 :
1081 198 : subroutine xcxalp(exc,npt,order,rspts,vxc, dvxc) ! dvxc is optional
1082 :
1083 : !Arguments ------------------------------------
1084 : !scalars
1085 : integer,intent(in) :: npt,order
1086 : !arrays
1087 : real(dp),intent(in) :: rspts(npt)
1088 : real(dp),intent(out) :: exc(npt),vxc(npt)
1089 : real(dp),intent(out),optional :: dvxc(npt)
1090 :
1091 : !Local variables-------------------------------
1092 : !Set value of alpha in "X-alpha" method
1093 : !scalars
1094 : integer :: ipt
1095 : real(dp),parameter :: alpha=1.0_dp
1096 : real(dp) :: dfac,efac,rs,rsm1,vfac
1097 : character(len=500) :: message
1098 :
1099 : ! *************************************************************************
1100 :
1101 : !Checks the values of order
1102 198 : if(order<0 .or. order>2)then
1103 : write(message, '(a,a,a,i3)' )&
1104 0 : & 'With X-alpha xc functional, the only',ch10,&
1105 0 : & 'allowed values for order are 0, 1 or 2, while it is found to be',order
1106 0 : ABI_BUG(message)
1107 : end if
1108 :
1109 : !Compute vfac=(3/(2*Pi))^(2/3)
1110 198 : vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
1111 : !Compute efac=(3/4)*vfac
1112 198 : efac=0.75_dp*vfac
1113 : !Compute dfac=(4*Pi/9)*vfac
1114 198 : dfac=(4.0_dp*pi/9.0_dp)*vfac
1115 :
1116 : !separate cases with respect to order
1117 198 : if(order==2) then
1118 : ! Loop over grid points
1119 0 : do ipt=1,npt
1120 0 : rs=rspts(ipt)
1121 0 : rsm1=1.0_dp/rs
1122 : ! compute energy density (hartree)
1123 0 : exc(ipt)=-alpha*efac*rsm1
1124 : ! compute potential (hartree)
1125 0 : vxc(ipt)=-alpha*vfac*rsm1
1126 : ! compute d(vxc)/d(rho) (hartree*bohr^3)
1127 0 : dvxc(ipt)=-alpha*dfac*rs**2
1128 : end do
1129 : else
1130 : ! Loop over grid points
1131 786630 : do ipt=1,npt
1132 786432 : rs=rspts(ipt)
1133 786432 : rsm1=1.0_dp/rs
1134 : ! compute energy density (hartree)
1135 786432 : exc(ipt)=-alpha*efac*rsm1
1136 : ! compute potential (hartree)
1137 786630 : vxc(ipt)=-alpha*vfac*rsm1
1138 : end do
1139 : end if
1140 : !
1141 198 : end subroutine xcxalp
1142 : !!***
1143 :
1144 : !!****f* ABINIT/xclb
1145 : !! NAME
1146 : !! xclb
1147 : !!
1148 : !! FUNCTION
1149 : !! Computes the GGA like part (vx_lb) of the Leeuwen-Baerends
1150 : !! exchange-correlation potential (vxc_lb) and adds it to the
1151 : !! lda exchange-correlation potential (vxc_lda) which
1152 : !! must be provided as input,
1153 : !! vxci <-- vxc_lb =: vxc_lda + vx_lb
1154 : !!
1155 : !! R van Leeuwen and EJ Baerends, Phys Rev A 49, 2421 (1994) [[cite:VanLeeuwen1994]]
1156 : !!
1157 : !! With respect to spin, the van Leeuwen-Baerends
1158 : !! potential is "exchange-like" : separate contributions from
1159 : !! spin up and spin down.
1160 : !!
1161 : !! INPUTS
1162 : !! npts= number of points to be computed
1163 : !! nspden=1 for unpolarized, 2 for spin-polarized
1164 : !! grho2_updn(npts,2*nspden-1)=square of the gradient of the spin-up,
1165 : !! and, if nspden==2, spin-down, and total density (Hartree/Bohr**2)
1166 : !! rho_updn(npts,nspden)=spin-up and spin-down density (Hartree/bohr**3)
1167 : !!
1168 : !! OUTPUT
1169 : !! (see side effects)
1170 : !!
1171 : !! SIDE EFFECTS
1172 : !! Input/Output:
1173 : !! vxci(npts,nspden)=input xc potential to which Leeuwen-Baerends correction
1174 : !! is added at output.
1175 : !!
1176 : !! SOURCE
1177 :
1178 0 : subroutine xclb(grho2_updn,npts,nspden,rho_updn,vxci)
1179 :
1180 : !Arguments ------------------------------------
1181 : !scalars
1182 : integer,intent(in) :: npts,nspden
1183 : !arrays
1184 : real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden)
1185 : real(dp),intent(inout) :: vxci(npts,nspden)
1186 :
1187 : !Local variables-------------------------------
1188 : !scalars
1189 : integer :: ipts,ispden
1190 : real(dp),parameter :: beta=0.05_dp
1191 : real(dp) :: density,density_gradient,density_t13,s_g_sq,scaled_gradient
1192 : real(dp) :: scaling_factor,vx_lb
1193 :
1194 : ! *************************************************************************
1195 :
1196 : !DEBUG
1197 : !write(std_out,*) ' %xclb: enter'
1198 : !ENDDEBUG
1199 :
1200 : !scale the spin densities for evaluating spin up or down exchange
1201 0 : scaling_factor=one
1202 0 : if(nspden == 2) scaling_factor=two
1203 :
1204 0 : do ispden=1,nspden
1205 :
1206 0 : do ipts=1,npts
1207 :
1208 0 : density= scaling_factor * rho_updn(ipts,ispden)
1209 0 : density_gradient= scaling_factor * sqrt(grho2_updn(ipts,ispden))
1210 :
1211 0 : density_t13= density**third
1212 0 : scaled_gradient= density_gradient/max(density*density_t13,1.e-12_dp)
1213 :
1214 0 : s_g_sq= scaled_gradient*scaled_gradient
1215 :
1216 : vx_lb= -beta*density_t13 * s_g_sq/ &
1217 0 : & (one+3.d0*beta* scaled_gradient*log(scaled_gradient+sqrt(one+s_g_sq*s_g_sq)))
1218 :
1219 0 : vxci(ipts,ispden)=vxci(ipts,ispden)+vx_lb
1220 : end do
1221 :
1222 : end do
1223 :
1224 0 : end subroutine xclb
1225 : !!***
1226 :
1227 : !!****f* ABINIT/xctfw
1228 : !! NAME
1229 : !! xctfw
1230 : !!
1231 : !! FUNCTION
1232 : !! Add gradient part of the Thomas-Fermi-Weizsacker functional
1233 : !! Perrot F., Phys. Rev. A20, 586-594 (1979) [[cite:Perrot1979]]
1234 : !!
1235 : !! INPUTS
1236 : !! ndvxcdgr= size of dvxcdgr(npts,ndvxcdgr)
1237 : !! npts= number of points to be computed
1238 : !! nspden=number if spin density component (necessarily 1 here)
1239 : !! grho2_updn(npts,2*nspden-1)=square of the gradient of the spin-up,
1240 : !! and, if nspden==2, spin-down, and total density (Hartree/Bohr**2),
1241 : !! only used if gradient corrected functional (option=2,-2,-4 and 4 or beyond)
1242 : !! rho_updn(npts,nspden)=spin-up and spin-down density (Hartree/bohr**3)
1243 : !! temp= electronic temperature
1244 : !!
1245 : !! SIDE EFFECTS
1246 : !! The following arrays are modified (gradient correction added):
1247 : !! dvxcdgr(npts,3)=partial derivative of the XC energy divided by the norm of the gradient
1248 : !! fxci(npts)=free energy energy density
1249 : !! tsxci(npts)=entropy energy density
1250 : !! vxci(npts,nspden)=exchange-correlation potential
1251 : !!
1252 : !! SOURCE
1253 :
1254 189 : subroutine xctfw(temp,fxci,tsxci,rho_updn,vxci,npts,nspden,dvxcdgr,ndvxcdgr,grho2_updn)
1255 :
1256 : !Arguments ------------------------------------
1257 : !scalars
1258 : integer,intent(in) :: ndvxcdgr,npts,nspden
1259 : real(dp),intent(in) :: temp
1260 : !arrays
1261 : real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden)
1262 : real(dp),intent(inout) :: dvxcdgr(npts,ndvxcdgr),fxci(npts),tsxci(npts),vxci(npts,nspden)
1263 :
1264 : !Local variables-------------------------------
1265 : !scalars
1266 : integer :: iperrot,ipts
1267 : logical :: has_dvxcdgr
1268 : real(dp) :: etfw,rho,rho_inv,rhomot,yperrot0,vtfw
1269 : real(dp) :: yperrot,uperrot,dyperrotdn,duperrotdyperrot
1270 : real(dp) :: hperrot,dhperrotdyperrot,dhperrotdn,dhperrotduperrot
1271 : !arrays
1272 : real(dp) :: wpy(0:7), wpu(0:7)
1273 189 : real(dp),allocatable :: rho_updnm1_3(:,:),exci(:)
1274 :
1275 : ! *************************************************************************
1276 :
1277 : !We would rather work with exc than with tsxc
1278 567 : ABI_MALLOC(exci,(npts))
1279 729378 : exci=fxci+tsxci
1280 189 : has_dvxcdgr=(ndvxcdgr/=0)
1281 :
1282 189 : yperrot0=1.666081101_dp
1283 :
1284 189 : wpy(0)=0.5_dp; wpy(1)=-0.1999176316_dp
1285 189 : wpy(2)=0.09765615709_dp; wpy(3)=-0.06237609924_dp
1286 189 : wpy(4)=0.05801466322_dp; wpy(5)=-0.04449287774_dp
1287 189 : wpy(6)=0.01903211697_dp; wpy(7)=-0.003284096926_dp
1288 :
1289 189 : wpu(0)=one/6._dp; wpu(1)=0.311590799_dp
1290 189 : wpu(2)=3.295662439_dp; wpu(3)=-29.22038326_dp
1291 189 : wpu(4)=116.1084531_dp; wpu(5)=-250.4543147_dp
1292 189 : wpu(6)=281.433688_dp; wpu(7)=-128.8784806_dp
1293 :
1294 567 : ABI_MALLOC(rho_updnm1_3,(npts,2))
1295 :
1296 189 : call invcb(rho_updn(:,1),rho_updnm1_3(:,1),npts)
1297 :
1298 729189 : do ipts=1,npts
1299 729000 : rho =rho_updn(ipts,1)
1300 729000 : rhomot=rho_updnm1_3(ipts,1)
1301 729000 : rho_inv=rhomot*rhomot*rhomot
1302 :
1303 729000 : yperrot=pi*pi/sqrt2/temp**1.5*two*rho
1304 729000 : uperrot=yperrot**(2./3.)
1305 :
1306 729000 : dyperrotdn=pi*pi/sqrt2/temp**1.5*2.0_dp
1307 :
1308 729000 : hperrot=zero
1309 729000 : dhperrotdyperrot=zero
1310 729000 : dhperrotduperrot=zero
1311 729000 : if(yperrot<=yperrot0)then
1312 0 : do iperrot=0,7
1313 0 : hperrot=hperrot+wpy(iperrot)*yperrot**iperrot
1314 0 : dhperrotdyperrot=dhperrotdyperrot+iperrot*wpy(iperrot)*yperrot**(iperrot-1)
1315 : end do
1316 0 : hperrot=one/12.0_dp*hperrot
1317 0 : dhperrotdyperrot=one/12.0_dp*dhperrotdyperrot
1318 0 : dhperrotdn=dhperrotdyperrot*dyperrotdn
1319 : else
1320 6561000 : do iperrot=0,7
1321 5832000 : hperrot=hperrot+wpu(iperrot)/uperrot**(2*iperrot)
1322 6561000 : dhperrotduperrot=dhperrotduperrot-2.*iperrot*wpu(iperrot)/uperrot**(2*iperrot+1)
1323 : end do
1324 729000 : hperrot=one/12.0_dp*hperrot
1325 729000 : dhperrotduperrot=one/12.0_dp*dhperrotduperrot
1326 729000 : duperrotdyperrot=two/3._dp/yperrot**(1./3.)
1327 729000 : dhperrotdn=dhperrotduperrot*duperrotdyperrot*dyperrotdn
1328 : end if
1329 :
1330 729000 : etfw=hperrot*grho2_updn(ipts,1)*rho_inv*rho_inv
1331 729000 : vtfw=-etfw + rho/hperrot*dhperrotdn*etfw
1332 :
1333 729000 : if(yperrot<=yperrot0)then
1334 0 : exci(ipts) = exci(ipts) + etfw + 1.5_dp*yperrot*dhperrotdyperrot*grho2_updn(ipts,1)*rho_inv*rho_inv
1335 : else
1336 729000 : exci(ipts) = exci(ipts) + etfw + uperrot*dhperrotduperrot*grho2_updn(ipts,1)*rho_inv*rho_inv
1337 : end if
1338 729000 : vxci(ipts,1) = vxci(ipts,1) + vtfw
1339 729000 : fxci(ipts) = fxci(ipts) + etfw
1340 729189 : if (has_dvxcdgr) dvxcdgr(ipts,1)= dvxcdgr(ipts,1)+two*hperrot*rho_inv
1341 : end do
1342 :
1343 729189 : tsxci=exci-fxci
1344 189 : ABI_FREE(rho_updnm1_3)
1345 189 : ABI_FREE(exci)
1346 :
1347 189 : end subroutine xctfw
1348 : !!***
1349 :
1350 : !!****f* ABINIT/xcksdt
1351 : !! NAME
1352 : !! xcksdt
1353 : !!
1354 : !! FUNCTION
1355 : !! Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input rho.
1356 : !!
1357 : !! NOTES
1358 : !! Karasiev-Sjostrom-Dufty-Trickey (KSDT) TLDA xc-functional
1359 : !! V.V. Karasiev, T. Sjostrom, J. Dufty, and S.B. Trickey, PRL 112, 076403 (2014) [[cite:Karasiev2014]]
1360 : !! Copyright (C) 2014 Orbital-free DFT group at University of Florida (V.V. Karasiev)
1361 : !!
1362 : !! INPUTS
1363 : !! npt=number of real space points on which density is provided
1364 : !! order=gives the maximal derivative of Exc computed.
1365 : !! rhor=value of electronic density at each point
1366 : !! rspts(npt)=Wigner-Seitz radii at each point
1367 : !! el_temp=electronic temperature (hartree)
1368 : !!
1369 : !! OUTPUT
1370 : !! exc(npt)=exchange-correlation free energy density (hartree)
1371 : !! tsxc(npt)=exchange-correlation entropy energy density (hartree)
1372 : !! vxc(npt)=xc potential (d($\rho$*exc)/d($\rho$)) (hartree)
1373 : !! if(order>1) dvxc(npt)=derivative d(vxc)/d($\rho$) (hartree*bohr^3)
1374 : !!
1375 : !! SOURCE
1376 56 : subroutine xcksdt(exc,tsxc,npt,order,rhor,rspts,el_temp,vxc,&
1377 : & dvxc) ! optional arguments
1378 : !Arguments ------------------------------------
1379 : !scalars
1380 : integer,intent(in) :: npt,order
1381 : real(dp),intent(in) :: el_temp
1382 : !arrays
1383 : real(dp),intent(in) :: rhor(npt),rspts(npt)
1384 : real(dp),intent(out) :: exc(npt),tsxc(npt),vxc(npt)
1385 : real(dp),intent(out),optional :: dvxc(npt)
1386 :
1387 : !Local variables ------------------------------
1388 : !scalars
1389 : integer :: ipt
1390 : real(dp) :: tfac,rs,rho,tempf,tred,fxc,einxc,tsxctmp
1391 : real(dp) :: drho
1392 : integer :: i
1393 : character(len=500) :: message
1394 : !arrays
1395 : real(dp) :: vxctmp(5)
1396 :
1397 : ! *************************************************************************
1398 :
1399 56 : tfac=(3._dp*pi**2)**(2._dp/3._dp)/2._dp
1400 : !Checks the values of order
1401 56 : if(order<0.or.order>2)then
1402 : write(message,'(a,a,a,i0)')&
1403 0 : & 'With Karasiev-Sjostrom-Dufty-Trickey xc functional, the only',ch10,&
1404 0 : & 'allowed values for order are 0, 1 or 2, while it is found to be',order
1405 0 : ABI_BUG(message)
1406 : end if
1407 :
1408 : !Checks the compatibility between the order and the presence of the optional arguments
1409 56 : if(order<=1.and.present(dvxc))then
1410 : write(message,'(a,a,a,i0)')&
1411 0 : & 'The order chosen does not need the presence',ch10,&
1412 0 : & 'of the vector dvxc, that is needed only with order=2, while we have',order
1413 0 : ABI_BUG(message)
1414 : end if
1415 :
1416 : !calculate exc=fxc, vxc, and tsxc (orders 1 and 2)
1417 : !Loop over grid points
1418 28814 : do ipt=1,npt
1419 28758 : rs=rspts(ipt)
1420 28758 : rho=rhor(ipt) !0.75_dp/pi/(rs**3)
1421 28758 : tempf=tfac*rho**(2._dp/3._dp) !(3._dp*pi**2*rho)**(2._dp/3._dp)/2._dp
1422 28758 : tred=el_temp/tempf
1423 28758 : call fxc_ksdt(fxc,vxc(ipt),einxc,tsxctmp,rs,tred,0)
1424 28758 : exc(ipt)=fxc
1425 28758 : tsxc(ipt)=tsxctmp
1426 28814 : if( (exc(ipt)/=exc(ipt)).or.(vxc(ipt)/=vxc(ipt)) ) then
1427 0 : exc(ipt)=0._dp
1428 0 : vxc(ipt)=0._dp
1429 : write(message, '(a,2d12.5)' )&
1430 0 : & 'fxc or vxc = NaN: rs,tred=',rs,tred
1431 0 : ABI_BUG(message)
1432 : endif
1433 : end do
1434 : !for order==2, use numerical derivative
1435 56 : if(order==2) then
1436 : ! Loop over grid points
1437 7007 : do ipt=1,npt
1438 7000 : drho=0.01_dp*rhor(ipt)
1439 42000 : do i=1,5
1440 35000 : rho=rhor(ipt)+drho*dble(i-3)
1441 35000 : rs=(0.75_dp/pi/rho)**(1._dp/3._dp) ! density
1442 35000 : tempf=tfac*rho**(2._dp/3._dp) !(3._dp*pi**2*rho)**(2._dp/3._dp)/2._dp
1443 35000 : tred=el_temp/tempf
1444 42000 : call fxc_ksdt(fxc,vxctmp(i),einxc,tsxctmp,rs,tred,0)
1445 : enddo
1446 7000 : dvxc(ipt)=vxctmp(1)-8._dp*vxctmp(2)+8._dp*vxctmp(4)-vxctmp(5)
1447 7000 : dvxc(ipt)=dvxc(ipt)/(12._dp*drho)
1448 7007 : if( dvxc(ipt)/=dvxc(ipt) ) then
1449 0 : dvxc(ipt)=0._dp
1450 : write(message, '(a,2d12.5)' )&
1451 0 : & 'dvxc = NaN: rs,tred=',rs,tred
1452 0 : ABI_BUG(message)
1453 7000 : elseif(dvxc(ipt)>huge(1._dp)) then
1454 0 : dvxc(ipt)=0._dp
1455 : write(message, '(a,2d12.5)' )&
1456 0 : & 'dvxc = Inf: rs,tred=',rs,tred
1457 0 : ABI_BUG(message)
1458 : endif
1459 : enddo
1460 : endif
1461 56 : end subroutine xcksdt
1462 : !!***
1463 :
1464 : !!****f* ABINIT/fxc_ksdt
1465 : !! NAME
1466 : !! fxc_ksdt
1467 : !!
1468 : !! FUNCTION
1469 : !! LDA XC free-energy parameterization from Monte Carlo data (unpol/pol)
1470 : !!
1471 : !! NOTES
1472 : !! Karasiev-Sjostrom-Dufty-Trickey (KSDT) TLDA xc-functional
1473 : !! V.V. Karasiev, T. Sjostrom, J. Dufty, and S.B. Trickey, PRL 112, 076403 (2014) [[cite:Karasiev2014]]
1474 : !! Copyright (C) 2014 Orbital-free DFT group at University of Florida (V.V. Karasiev)
1475 : !!
1476 : !! INPUTS
1477 : !! rs=Wigner-Seitz radius (bohr)
1478 : !! t=reduced temperature
1479 : !! iz=spin polarization (0 - spin-unpolarized, 1 - fully polarized)
1480 : !!
1481 : !! OUTPUT
1482 : !! fxc=exchange-correlation free energy per particle (hartree)
1483 : !! vxc=xc potential (d($\rho$*exc)/d($\rho$)) (hartree)
1484 : !! exc=exchange-correlation internal energy per particle (hartree)
1485 : !! tsxc=exchange-correlation entropy energy per particle (hartree)
1486 : !!
1487 : !! SOURCE
1488 158274 : subroutine fxc_ksdt(fxc,vxc,exc,tsxc,rs,t,iz)
1489 : !Arguments ------------------------------------
1490 : !scalars
1491 : integer,intent(in) :: iz
1492 : real(dp),intent(out) :: fxc,vxc,exc,tsxc
1493 : real(dp),intent(in) :: rs,t
1494 : !Local variables ------------------------------
1495 : !scalars
1496 : real(dp),parameter :: onethird=1._dp/3._dp
1497 : real(dp),parameter :: threehalf=1.5_dp
1498 : real(dp),parameter :: lambda=(4._dp/9._dp/pi)**onethird
1499 : real(dp),parameter :: a0=1._dp/(pi*lambda)
1500 : real(dp) :: aa,daa,bb,dbb,cc,dcc,dd,ddd,ee,dee,tempF,sxc
1501 : real(dp) :: dtdn,tanht,dtanht,tanhsqrt,dtanhsqrt,f1
1502 : real(dp) :: num,dnum,den,dden,dnumdrs,ddendrs,n,drsdn,omega
1503 : !arrays
1504 : real(dp) :: a(6),b(0:1,4),c(0:1,3),d(0:1,5),e(0:1,5)
1505 :
1506 : ! *************************************************************************
1507 :
1508 : data a/0.75_dp,3.04363_dp,-0.092270_dp,1.70350_dp,8.31051_dp,5.1105_dp/
1509 : !
1510 : data b(0,:)/0.342554_dp,9.141315_dp,0.448483_dp,18.553096_dp/
1511 : data c(0,:)/0.875130_dp,-0.256320_dp,0.953988_dp/
1512 : data d(0,:)/0.725917_dp,2.237347_dp,0.280748_dp,4.185911_dp,0.692183_dp/
1513 : data e(0,:)/0.255415_dp,0.931933_dp,0.115398_dp,17.234117_dp,0.451437_dp/
1514 : !
1515 : data b(1,:)/0.329001_dp,111.598308_dp,0.537053_dp,105.086663_dp/
1516 : data c(1,:)/0.848930_dp,0.167952_dp,0.088820_dp/
1517 : data d(1,:)/0.551330_dp,180.213159_dp,134.486231_dp,103.861695_dp,17.750710_dp/
1518 : data e(1,:)/0.153124_dp,19.543945_dp,43.400337_dp,120.255145_dp,15.662836_dp/
1519 : !
1520 158274 : if(iz==0) then
1521 : omega=1._dp
1522 0 : elseif(iz==1) then
1523 0 : omega=2._dp**onethird
1524 : endif
1525 : !
1526 158274 : if(t==0._dp) then
1527 : !fxc
1528 0 : f1=-1._dp/rs
1529 0 : num=omega*a0*a(1)+b(iz,1)*sqrt(rs)+c(iz,1)*e(iz,1)*rs
1530 0 : den=1._dp+d(iz,1)*sqrt(rs)+e(iz,1)*rs
1531 0 : fxc=f1*num/den
1532 : !
1533 0 : dnumdrs=b(iz,1)/sqrt(rs)/2._dp+c(iz,1)*e(iz,1)
1534 0 : ddendrs=d(iz,1)/sqrt(rs)/2._dp+e(iz,1)
1535 :
1536 0 : n=3._dp/(4._dp*pi*rs**3) ! density
1537 0 : drsdn=-onethird*rs/n ! (drs/dn)
1538 : !Fxc=n*fxc=n*f1*num/den=n*A*B*C
1539 : !dFxc/dn=fxc+n*(dA/dn)*B*C+n*a*(dB/dn)*C+n*a*B*(dC/dn)
1540 : vxc=fxc+(onethird*f1)*num/den & ! fxc+n*(dA/dn)*B*C
1541 : & + n*f1*(dnumdrs*drsdn)/den & ! n*a*(dB/dn)*C
1542 0 : & - n*f1*num*(ddendrs*drsdn)/den**2 ! n*a*B*(dC/dn)
1543 0 : exc=fxc
1544 0 : tsxc=zero
1545 : else
1546 158274 : tanht=tanh(1._dp/t)
1547 158274 : tanhsqrt=tanh(1._dp/sqrt(t))
1548 158274 : dtanht=(tanht**2-1._dp)/t**2 !d/dt tanh(1/t)
1549 158274 : dtanhsqrt=(tanhsqrt**2-1._dp)/t**threehalf/2._dp !d/dt tanh(1/sqrt(t))
1550 : !
1551 : ! a(t)
1552 158274 : num=a(1)+a(2)*t**2+a(3)*t**3+a(4)*t**4
1553 158274 : den=1._dp+a(5)*t**2+a(6)*t**4
1554 : !
1555 158274 : dnum=a(2)*2._dp*t+a(3)*3._dp*t**2+a(4)*4._dp*t**3
1556 158274 : dden=a(5)*2._dp*t+a(6)*4._dp*t**3
1557 : !
1558 158274 : aa=a0*tanht*num/den
1559 158274 : daa=a0*(dtanht*num/den+tanht*dnum/den-tanht*num*dden/den**2)
1560 : !
1561 : ! b(t)
1562 158274 : num=b(iz,1)+b(iz,2)*t**2+b(iz,3)*t**4
1563 158274 : den=1._dp+b(iz,4)*t**2+omega*sqrt(3._dp)*b(iz,3)/sqrt(2._dp*lambda**2)*t**4
1564 : !
1565 158274 : dnum=b(iz,2)*2._dp*t+b(iz,3)*4._dp*t**3
1566 158274 : dden=b(iz,4)*2._dp*t+omega*sqrt(3._dp)*b(iz,3)/sqrt(2._dp*lambda**2)*4._dp*t**3
1567 : !
1568 158274 : bb=tanhsqrt*num/den
1569 158274 : dbb=dtanhsqrt*num/den+tanhsqrt*dnum/den-tanhsqrt*num*dden/den**2
1570 : !
1571 : ! d(t)
1572 158274 : num=d(iz,1)+d(iz,2)*t**2+d(iz,3)*t**4
1573 158274 : den=1._dp+d(iz,4)*t**2+d(iz,5)*t**4
1574 : !
1575 158274 : dnum=d(iz,2)*2._dp*t+d(iz,3)*4._dp*t**3
1576 158274 : dden=d(iz,4)*2._dp*t+d(iz,5)*4._dp*t**3
1577 : !
1578 158274 : dd=tanhsqrt*num/den
1579 158274 : ddd=dtanhsqrt*num/den+tanhsqrt*dnum/den-tanhsqrt*num*dden/den**2
1580 : !
1581 : ! e(t)
1582 158274 : num=e(iz,1)+e(iz,2)*t**2+e(iz,3)*t**4
1583 158274 : den=1._dp+e(iz,4)*t**2+e(iz,5)*t**4
1584 : !
1585 158274 : dnum=e(iz,2)*2._dp*t+e(iz,3)*4._dp*t**3
1586 158274 : dden=e(iz,4)*2._dp*t+e(iz,5)*4._dp*t**3
1587 : !
1588 158274 : ee=tanht*num/den
1589 158274 : dee=dtanht*num/den+tanht*dnum/den-tanht*num*dden/den**2
1590 : !
1591 : ! c(t)
1592 158274 : num=c(iz,1)+c(iz,2)*exp(-c(iz,3)/t)
1593 158274 : dnum=c(iz,2)*c(iz,3)*exp(-c(iz,3)/t)/t**2
1594 158274 : cc=num*ee
1595 158274 : dcc=dnum*ee+num*dee
1596 : !
1597 : !fxc
1598 158274 : f1=-1._dp/rs
1599 158274 : num=omega*aa+bb*sqrt(rs)+cc*rs
1600 158274 : den=1._dp+dd*sqrt(rs)+ee*rs
1601 158274 : fxc=f1*num/den
1602 : !
1603 158274 : dnum=omega*daa+dbb*sqrt(rs)+dcc*rs
1604 158274 : dnumdrs=bb/sqrt(rs)/2._dp+cc
1605 158274 : dden=ddd*sqrt(rs)+dee*rs
1606 158274 : ddendrs=dd/sqrt(rs)/2._dp+ee
1607 :
1608 158274 : n=3._dp/(4._dp*pi*rs**3) ! density
1609 158274 : tempF = (3._dp*pi**2*n)**(2._dp/3._dp)/2._dp*omega**2
1610 158274 : dtdn = -2._dp/3._dp*t/n ! (dt/dn)
1611 158274 : drsdn=-onethird*rs/n ! (drs/dn)
1612 : !Fxc=n*fxc=n*f1*num/den=n*A*B*C
1613 : !dFxc/dn=fxc+n*(dA/dn)*B*C+n*a*(dB/dn)*C+n*a*B*(dC/dn)
1614 : vxc=fxc+(onethird*f1)*num/den & ! fxc+n*(dA/dn)*B*C
1615 : & + n*f1*(dnum*dtdn+dnumdrs*drsdn)/den & ! n*a*(dB/dn)*C
1616 158274 : & - n*f1*num*(dden*dtdn+ddendrs*drsdn)/den**2 ! n*a*B*(dC/dn)
1617 : sxc=f1*(dnum)/den & ! A*(dB/dn)*C
1618 158274 : & -f1*num*(dden)/den**2 ! A*B*(dC/dn)
1619 158274 : sxc=-sxc/tempF !sxc=-(t/T)*dfxc/dt=-(1/tempF)*dfxc/dt
1620 158274 : exc=fxc+t*tempF*sxc !exc=fxc+T*sxc=fxc+(t*tempF)*sxc
1621 158274 : tsxc=t*tempF*sxc
1622 : endif
1623 158274 : end subroutine fxc_ksdt
1624 : !!***
1625 :
1626 : !!****f* ABINIT/fex_ksdt
1627 : !! NAME
1628 : !! fex_ksdt
1629 : !!
1630 : !! FUNCTION
1631 : !! Returns exchange energy per electron from KSDT xc functional
1632 : !!
1633 : !! NOTES
1634 : !! Karasiev-Sjostrom-Dufty-Trickey (KSDT) TLDA xc-functional
1635 : !! V.V. Karasiev, T. Sjostrom, J. Dufty, and S.B. Trickey, PRL 112, 076403 (2014) [[cite:Karasiev2014]]
1636 : !! Copyright (C) 2014 Orbital-free DFT group at University of Florida (V.V. Karasiev)
1637 : !!
1638 : !! INPUTS
1639 : !! rs=Wigner-Seitz radius (Bohr)
1640 : !! degauss=setup temperature (Rydberg)
1641 : !!
1642 : !! OUTPUT
1643 : !! fx=exchange free energy per particle (hartree)
1644 : !! einx=exchange internal energy per particle (hartree)
1645 : !! tsx=exchange entropy energy per particle (hartree)
1646 : !! vx=exchange potential
1647 : !!
1648 : !! SOURCE
1649 64758 : subroutine fex_ksdt(rs,fx,einx,tsx,vx,degauss)
1650 : !Arguments ------------------------------------
1651 : !scalars
1652 : real(dp),intent(in) :: rs,degauss
1653 : real(dp),intent(out) :: fx,einx,tsx,vx
1654 : !Local variables ------------------------------
1655 : !scalars
1656 : real(dp),parameter :: twothird=two/three
1657 : real(dp) :: ex0,vx0,tF,t,dtdn,rho,sx,Ax,dAx,d2Ax,f_slater,alpha_slater
1658 :
1659 : ! *************************************************************************
1660 :
1661 64758 : rho=3._dp/(4._dp*pi*rs**3)
1662 64758 : tF=(3._dp*pi**2*rho)**twothird/2._dp !tF=Fermi temperature for spin-unpol case
1663 64758 : t = degauss/2.0_dp/tF
1664 : !ef=1.841584276_dp/rs/rs
1665 : !tred=degauss/2.0_dp/ef
1666 64758 : dtdn = -twothird*t/rho ! (dt/dn)
1667 :
1668 64758 : f_slater=-0.687247939924714d0
1669 64758 : alpha_slater=twothird
1670 64758 : ex0=f_slater*alpha_slater/rs
1671 64758 : vx0=four/three*f_slater*alpha_slater/rs
1672 :
1673 64758 : call tildeAx(t,Ax,dAx,d2Ax)
1674 64758 : fx = ex0*Ax ! exchange free-energy per electron
1675 64758 : vx = vx0*Ax + rho*ex0*dAx*dtdn ! d(n*ex0*Ax)/dn = d(n*ex0)/dn + n*ex0*(dAx/dtred)*(dtred/dn)
1676 64758 : sx = -ex0*dAx/tF ! entropy per electron
1677 64758 : einx = fx + t*tF*sx ! internal energy per electron
1678 64758 : tsx = t*tF*sx ! T*entropy per electron
1679 64758 : end subroutine fex_ksdt
1680 : !!***
1681 :
1682 : !!****f* ABINIT/fec_ksdt
1683 : !! NAME
1684 : !! fec_ksdt
1685 : !!
1686 : !! FUNCTION
1687 : !! Returns correlation energy per electron from KSDT xc functional
1688 : !!
1689 : !! NOTES
1690 : !! Karasiev-Sjostrom-Dufty-Trickey (KSDT) TLDA xc-functional
1691 : !! V.V. Karasiev, T. Sjostrom, J. Dufty, and S.B. Trickey, PRL 112, 076403 (2014) [[cite:Karasiev2014]]
1692 : !! Copyright (C) 2014 Orbital-free DFT group at University of Florida (V.V. Karasiev)
1693 : !!
1694 : !! INPUTS
1695 : !! rs=Wigner-Seitz radius (Bohr)
1696 : !! degauss=setup temperature (Rydberg)
1697 : !!
1698 : !! OUTPUT
1699 : !! fc=correlation free energy per particle (hartree)
1700 : !! einc=correlation internal energy per particle (hartree)
1701 : !! tsc=correlation entropy energy per particle (hartree)
1702 : !! vc=correlation potential
1703 : !!
1704 : !! SOURCE
1705 129516 : subroutine fec_ksdt(rs,fc,einc,tsc,vc,degauss)
1706 : !Arguments ------------------------------------
1707 : !scalars
1708 : real(dp),intent(in) :: rs,degauss
1709 : real(dp),intent(out) :: fc,einc,tsc,vc
1710 : !Local variables ------------------------------
1711 : !scalars
1712 : real(dp),parameter :: twothird=two/three
1713 : real(dp) :: fxc,vxc,einxc,tsxc,fx,einx,tsx,vx,t,ef
1714 :
1715 : ! *************************************************************************
1716 :
1717 64758 : ef=1.841584276_dp/rs/rs*1.d0**(2.d0/3.d0)
1718 64758 : t=degauss/2.0_dp/ef
1719 64758 : call fxc_ksdt(fxc,vxc,einxc,tsxc,rs,t,0)
1720 64758 : call fex_ksdt(rs,fx,einx,tsx,vx,degauss)
1721 64758 : fc=fxc-fx
1722 64758 : einc=einxc-einx
1723 64758 : tsc=tsxc-tsx
1724 64758 : vc=vxc-vx
1725 64758 : end subroutine fec_ksdt
1726 : !!***
1727 :
1728 : end module m_xclda
1729 : !!***
|
