Branch data 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 : 209904 : subroutine xcspol(exc,npts,nspden,order,rspts,vxc,zeta,ndvxc,& !Mandatory arguments
240 : 5666 : & 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 [ - + ]: 104952 : 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 [ + + - + ]: 104952 : 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 [ + + + + : 104952 : 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 [ + + + + : 104952 : 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 : 104952 : 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 [ + + ]: 104952 : if (nspden==1) then
355 : : ! separate cases with respect to order
356 [ + + ]: 81503 : 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 [ + + ]: 81425 : else if(order**2>1) then
399 : : ! No spin-polarization so skip steps related to zeta not 0
400 [ + + ]: 13108382 : do ipts=1,npts
401 : :
402 : 13104796 : rs=rspts(ipts)
403 : 13104796 : n1=a0p+rs*(a1p+rs*(a2p+rs*a3p))
404 : 13104796 : d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p)))
405 : 13104796 : d1m1=one/d1
406 : :
407 : : ! Exchange-correlation energy
408 : 13104796 : excipt=-n1*d1m1
409 : 13104796 : exc(ipts)=excipt
410 : :
411 : : ! Exchange-correlation potential
412 : 13104796 : dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p))
413 : 13104796 : dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p)))
414 : :
415 : : ! dexcdrs is d(exc)/d(rs)
416 : 13104796 : dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
417 : 13104796 : vxc(ipts,1)=excipt-third*rs*dexcdrs
418 : :
419 : : ! If the exchange-correlation kernel is needed
420 : 13104796 : d2n1drs2=2._dp*a2p+rs*(6._dp*a3p)
421 : 13104796 : d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p))
422 : : ! d2excdrs2 is d2(exc)/d(rs)2
423 : 13104796 : d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1
424 : 13104796 : dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2)
425 : : ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs)
426 : 13108382 : 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 [ + + ]: 256645279 : do ipts=1,npts
432 : :
433 : 256567440 : rs=rspts(ipts)
434 : 256567440 : n1=a0p+rs*(a1p+rs*(a2p+rs*a3p))
435 : 256567440 : d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p)))
436 : 256567440 : d1m1=one/d1
437 : :
438 : : ! Exchange-correlation energy
439 : 256567440 : excipt=-n1*d1m1
440 : 256567440 : exc(ipts)=excipt
441 : :
442 : : ! Exchange-correlation potential
443 : 256567440 : dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p))
444 : 256567440 : dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p)))
445 : :
446 : : ! dexcdrs is d(exc)/d(rs)
447 : 256567440 : dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1
448 : 256645279 : 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 [ + + ]: 104952 : 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 : : !!***
|
