Coverage for gpaw/xc/fxc_kernels.py: 43%
84 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1import numpy as np
4def get_fHxc_Gr(xcflags, rs, q, qF, s2_g):
5 xc = xcflags.xc
6 if xc == 'rALDA':
7 return fHxc_ralda(q, qF)
9 elif xcflags.is_apbe:
10 return fHxc_apbe(rs, q, s2_g)
11 else:
12 raise ValueError(f'Unknown xc: {xc}')
15def fHxc_ralda(q, qF):
16 # rALDA (exchange only) kernel
17 # Olsen and Thygesen, Phys. Rev. B 88, 115131 (2013)
18 # ALDA up to 2*qF, -vc for q >2qF (such that fHxc vanishes)
20 rxalda_A = 0.25
21 rxalda_qcut = qF * np.sqrt(1.0 / rxalda_A)
23 # construct fHxc(k,r)
24 return (0.5 + 0.0j) * (
25 (1.0 + np.sign(rxalda_qcut - q[:, np.newaxis])) *
26 (1.0 + (-1.0) * rxalda_A * (q[:, np.newaxis] / qF)**2.0))
29def fHxc_apbe(rs, q, s2_g):
30 # Olsen and Thygesen, Phys. Rev. Lett. 112, 203001 (2014)
31 # Exchange only part of the PBE XC kernel, neglecting the terms
32 # arising from the variation of the density gradient
33 # i.e. second functional derivative
34 # d2/drho^2 -> \partial^2/\partial rho^2 at fixed \nabla \rho
36 fxc_PBE = get_pbe_fxc(
37 pbe_rho=3.0 / (4.0 * np.pi * rs**3.0),
38 pbe_s2_g=s2_g)
39 rxapbe_qcut = np.sqrt(-4.0 * np.pi / fxc_PBE)
41 return (0.5 + 0.0j) * (
42 (1.0 + np.sign(rxapbe_qcut - q[:, np.newaxis])) *
43 (1.0 + fxc_PBE / (4.0 * np.pi) * (q[:, np.newaxis])**2.0))
46def get_fspinHxc_Gr_rALDA(qF, q):
47 rxalda_A = 0.25
48 rxalda_qcut = qF * np.sqrt(1.0 / rxalda_A)
50 fspinHxc_Gr = ((0.5 + 0.0j) *
51 (1.0 + np.sign(rxalda_qcut - q[:, np.newaxis])) *
52 (-1.0) * rxalda_A * (q[:, np.newaxis] / qF)**2.0)
53 return fspinHxc_Gr
56def get_fspinHxc_Gr_rAPBE(rs, q, s2_g):
57 fxc_PBE = get_pbe_fxc(pbe_rho=(3.0 / (4.0 * np.pi * rs**3.0)),
58 pbe_s2_g=s2_g)
59 rxapbe_qcut = np.sqrt(-4.0 * np.pi / fxc_PBE)
60 fspinHxc_Gr = ((0.5 + 0.0j) *
61 (1.0 + np.sign(rxapbe_qcut - q[:, np.newaxis])) *
62 fxc_PBE / (4.0 * np.pi) * q[:, np.newaxis]**2.0)
63 return fspinHxc_Gr
66def get_heg_A(rs):
67 # Returns the A coefficient, where the
68 # q ->0 limiting value of static fxc
69 # of the HEG = -\frac{4\pi A }{q_F^2} = f_{xc}^{ALDA}.
70 # We need correlation energy per electron and first and second derivs
71 # w.r.t. rs
72 # See for instance Moroni, Ceperley and Senatore,
73 # Phys. Rev. Lett. 75, 689 (1995)
74 # (and also Kohn and Sham, Phys. Rev. 140, A1133 (1965) equation 2.7)
76 # Exchange contribution
77 heg_A = 0.25
79 # Correlation contribution
80 A_ec, A_dec, A_d2ec = get_pw_lda(rs)
81 heg_A += (1.0 / 27.0 * rs**2.0 * (9.0 * np.pi / 4.0)**(2.0 / 3.0) *
82 (2 * A_dec - rs * A_d2ec))
84 return heg_A
87def get_heg_B(rs):
88 # Returns the B coefficient, where the
89 # q -> oo limiting behaviour of static fxc
90 # of the HEG is -\frac{4\pi B}{q^2} - \frac{4\pi C}{q_F^2}.
91 # Use the parametrisation of Moroni, Ceperley and Senatore,
92 # Phys. Rev. Lett. 75, 689 (1995)
94 mcs_xs = np.sqrt(rs)
96 mcs_a = (1.0, 2.15, 0.0, 0.435)
97 mcs_b = (3.0, 1.57, 0.0, 0.409)
99 mcs_num = 0
101 for mcs_j, mcs_coeff in enumerate(mcs_a):
102 mcs_num += mcs_coeff * mcs_xs**mcs_j
104 mcs_denom = 0
106 for mcs_j, mcs_coeff in enumerate(mcs_b):
107 mcs_denom += mcs_coeff * mcs_xs**mcs_j
109 heg_B = mcs_num / mcs_denom
110 return heg_B
113def get_heg_C(rs):
114 # Returns the C coefficient, where the
115 # q -> oo limiting behaviour of static fxc
116 # of the HEG is -\frac{4\pi B}{q^2} - \frac{4\pi C}{q_F^2}.
117 # Again see Moroni, Ceperley and Senatore,
118 # Phys. Rev. Lett. 75, 689 (1995)
120 C_ec, C_dec, Cd2ec = get_pw_lda(rs)
122 heg_C = ((-1.0) * np.pi**(2.0 / 3.0) * (1.0 / 18.0)**(1.0 / 3.0) *
123 (rs * C_ec + rs**2.0 * C_dec))
125 return heg_C
128def get_heg_D(rs):
129 # Returns a 'D' coefficient, where the
130 # q->0 omega -> oo limiting behaviour
131 # of the frequency dependent fxc is -\frac{4\pi D}{q_F^2}
132 # see Constantin & Pitarke Phys. Rev. B 75, 245127 (2007) equation 7
134 D_ec, D_dec, D_d2ec = get_pw_lda(rs)
136 # Exchange contribution
137 heg_D = 0.15
138 # Correlation contribution
139 heg_D += ((9.0 * np.pi / 4.0)**(2.0 / 3.0) * rs / 3.0 *
140 (22.0 / 15.0 * D_ec + 26.0 / 15.0 * rs * D_dec))
141 return heg_D
144def get_pw_lda(rs):
145 # Returns LDA correlation energy and its first and second
146 # derivatives with respect to rs according to the parametrisation
147 # of Perdew and Wang, Phys. Rev. B 45, 13244 (1992)
149 pw_A = 0.031091
150 pw_alp = 0.21370
151 pw_beta = (7.5957, 3.5876, 1.6382, 0.49294)
153 pw_logdenom = 2.0 * pw_A * (
154 pw_beta[0] * rs**0.5 + pw_beta[1] * rs**1.0 +
155 pw_beta[2] * rs**1.5 + pw_beta[3] * rs**2.0)
157 pw_dlogdenom = 2.0 * pw_A * (0.5 * pw_beta[0] * rs**(-0.5) +
158 1.0 * pw_beta[1] +
159 1.5 * pw_beta[2] * rs**0.5 +
160 2.0 * pw_beta[3] * rs)
162 pw_d2logdenom = 2.0 * pw_A * (-0.25 * pw_beta[0] * rs**(-1.5) +
163 0.75 * pw_beta[2] * rs**(-0.5) +
164 2.0 * pw_beta[3])
166 pw_logarg = 1.0 + 1.0 / pw_logdenom
167 pw_dlogarg = (-1.0) / (pw_logdenom**2.0) * pw_dlogdenom
168 pw_d2logarg = 2.0 / (pw_logdenom**3.0) * (pw_dlogdenom**2.0)
169 pw_d2logarg += (-1.0) / (pw_logdenom**2.0) * pw_d2logdenom
171 # pw_ec = the correlation energy (per electron)
172 pw_ec = -2.0 * pw_A * (1 + pw_alp * rs) * np.log(pw_logarg)
174 # pw_dec = first derivative
176 pw_dec = -2.0 * pw_A * (1 + pw_alp * rs) * pw_dlogarg / pw_logarg
177 pw_dec += (-2.0 * pw_A * pw_alp) * np.log(pw_logarg)
179 # pw_d2ec = second derivative
181 pw_d2ec = (-2.0) * pw_A * pw_alp * pw_dlogarg / pw_logarg
182 pw_d2ec += (-2.0) * pw_A * ((1 + pw_alp * rs) *
183 (pw_d2logarg / pw_logarg -
184 (pw_dlogarg**2.0) / (pw_logarg**2.0)))
185 pw_d2ec += (-2.0 * pw_A * pw_alp) * pw_dlogarg / pw_logarg
187 return pw_ec, pw_dec, pw_d2ec
190def get_pbe_fxc(pbe_rho, pbe_s2_g):
191 return get_pbe_fxc_and_intermediate_derivatives(pbe_rho, pbe_s2_g)[0]
194def get_pbe_fxc_and_intermediate_derivatives(pbe_rho, pbe_s2_g):
195 # The intermediate derivatives are only used for testing.
196 pbe_kappa = 0.804
197 pbe_mu = 0.2195149727645171
199 pbe_denom_g = 1.0 + pbe_mu * pbe_s2_g / pbe_kappa
201 F_g = 1.0 + pbe_kappa - pbe_kappa / pbe_denom_g
202 Fn_g = -8.0 / 3.0 * pbe_mu * pbe_s2_g / pbe_rho / pbe_denom_g**2.0
203 Fnn_g = (-11.0 / 3.0 / pbe_rho * Fn_g -
204 2.0 / pbe_kappa * Fn_g**2.0 * pbe_denom_g)
206 e_g = -3.0 / 4.0 * (3.0 / np.pi)**(1.0 / 3.0) * pbe_rho**(4.0 / 3.0)
207 v_g = 4.0 / 3.0 * e_g / pbe_rho
208 f_g = 1.0 / 3.0 * v_g / pbe_rho
210 pbe_f_g = f_g * F_g + 2.0 * v_g * Fn_g + e_g * Fnn_g
211 return pbe_f_g, F_g, Fn_g, Fnn_g