Coverage for gpaw/xc/gga.py: 91%
345 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-12 00:18 +0000
1from math import pi
3import numpy as np
5from gpaw.xc.lda import (calculate_paw_correction, stress_integral,
6 stress_lda_term)
7from gpaw.utilities.blas import axpy
8from gpaw.fd_operators import Gradient
9from gpaw.sphere.lebedev import Y_nL, weight_n
10from gpaw.xc.pawcorrection import rnablaY_nLv
11from gpaw.xc.functional import XCFunctional
14def stress_gga_term(gd, sigma_xg, gradn_svg, dedsigma_xg):
15 P = 0
16 nspins = gradn_svg.shape[0]
17 for sigma_g, dedsigma_g in zip(sigma_xg, dedsigma_xg):
18 P -= 2 * stress_integral(gd, sigma_g, dedsigma_g)
19 stress_vv = P * np.eye(3)
20 for v1 in range(3):
21 for v2 in range(3):
22 stress_vv[v1, v2] -= 2 * stress_integral(
23 gd, gradn_svg[0, v1] * gradn_svg[0, v2], dedsigma_xg[0]
24 )
25 if nspins == 2:
26 stress_vv[v1, v2] -= 2 * stress_integral(
27 gd, gradn_svg[0, v1] * gradn_svg[1, v2], dedsigma_xg[1]
28 )
29 stress_vv[v1, v2] -= 2 * stress_integral(
30 gd, gradn_svg[1, v1] * gradn_svg[1, v2], dedsigma_xg[2]
31 )
32 return stress_vv
35class GGARadialExpansion:
36 def __init__(self, rcalc, *args):
37 self.rcalc = rcalc
38 self.args = args
40 def __call__(self, rgd, D_sLq, n_qg, nc0_sg):
41 n_sLg = np.dot(D_sLq, n_qg)
42 n_sLg[:, 0] += nc0_sg
44 dndr_sLg = np.empty_like(n_sLg)
45 for n_Lg, dndr_Lg in zip(n_sLg, dndr_sLg):
46 for n_g, dndr_g in zip(n_Lg, dndr_Lg):
47 rgd.derivative(n_g, dndr_g)
49 nspins, Lmax, nq = D_sLq.shape
50 dEdD_sqL = np.zeros((nspins, nq, Lmax))
52 E = 0.0
53 for n, Y_L in enumerate(Y_nL[:, :Lmax]):
54 w = weight_n[n]
55 rnablaY_Lv = rnablaY_nLv[n, :Lmax]
56 e_g, dedn_sg, b_vsg, dedsigma_xg = \
57 self.rcalc(rgd, n_sLg, Y_L, dndr_sLg, rnablaY_Lv, n,
58 *self.args)
59 dEdD_sqL += np.dot(rgd.dv_g * dedn_sg,
60 n_qg.T)[:, :, np.newaxis] * (w * Y_L)
61 dedsigma_xg *= rgd.dr_g
62 B_vsg = dedsigma_xg[::2] * b_vsg
63 if nspins == 2:
64 B_vsg += 0.5 * dedsigma_xg[1] * b_vsg[:, ::-1]
65 B_vsq = np.dot(B_vsg, n_qg.T)
66 dEdD_sqL += 8 * pi * w * np.inner(rnablaY_Lv, B_vsq.T).T
67 E += w * rgd.integrate(e_g)
69 return E, dEdD_sqL
72# First part of gga_calculate_radial - initializes some quantities.
73def radial_gga_vars(rgd, n_sLg, Y_L, dndr_sLg, rnablaY_Lv):
74 nspins = len(n_sLg)
76 n_sg = np.dot(Y_L, n_sLg)
78 a_sg = np.dot(Y_L, dndr_sLg)
79 b_vsg = np.dot(rnablaY_Lv.T, n_sLg)
81 sigma_xg = rgd.empty(2 * nspins - 1)
82 sigma_xg[::2] = (b_vsg ** 2).sum(0)
83 if nspins == 2:
84 sigma_xg[1] = (b_vsg[:, 0] * b_vsg[:, 1]).sum(0)
85 sigma_xg[:, 1:] /= rgd.r_g[1:] ** 2
86 sigma_xg[:, 0] = sigma_xg[:, 1]
87 sigma_xg[::2] += a_sg ** 2
88 if nspins == 2:
89 sigma_xg[1] += a_sg[0] * a_sg[1]
91 e_g = rgd.empty()
92 dedn_sg = rgd.zeros(nspins)
93 dedsigma_xg = rgd.zeros(2 * nspins - 1)
94 return e_g, n_sg, dedn_sg, sigma_xg, dedsigma_xg, a_sg, b_vsg
97def add_radial_gradient_correction(rgd, sigma_xg, dedsigma_xg, a_sg):
98 nspins = len(a_sg)
99 vv_sg = sigma_xg[:nspins] # reuse array
100 for s in range(nspins):
101 rgd.derivative2(-2 * rgd.dv_g * dedsigma_xg[2 * s] * a_sg[s],
102 vv_sg[s])
103 if nspins == 2:
104 v_g = sigma_xg[2]
105 rgd.derivative2(rgd.dv_g * dedsigma_xg[1] * a_sg[1], v_g)
106 vv_sg[0] -= v_g
107 rgd.derivative2(rgd.dv_g * dedsigma_xg[1] * a_sg[0], v_g)
108 vv_sg[1] -= v_g
110 vv_sg[:, 1:] /= rgd.dv_g[1:]
111 vv_sg[:, 0] = vv_sg[:, 1]
112 return vv_sg
115class GGARadialCalculator:
116 def __init__(self, kernel):
117 self.kernel = kernel
119 def __call__(self, rgd, n_sLg, Y_L, dndr_sLg, rnablaY_Lv, n):
120 (e_g, n_sg, dedn_sg, sigma_xg, dedsigma_xg, a_sg,
121 b_vsg) = radial_gga_vars(rgd, n_sLg, Y_L, dndr_sLg, rnablaY_Lv)
122 self.kernel.calculate(e_g, n_sg, dedn_sg, sigma_xg, dedsigma_xg)
123 vv_sg = add_radial_gradient_correction(rgd, sigma_xg,
124 dedsigma_xg, a_sg)
125 return e_g, dedn_sg + vv_sg, b_vsg, dedsigma_xg
128def calculate_sigma(gd, grad_v, n_sg):
129 r"""Calculate sigma(r) and grad n(r).
130 _ __ _ 2 __ _
131 Returns sigma(r) = |\/ n(r)| and \/ n (r).
133 With multiple spins, sigma has the three elements
135 _ __ _ 2
136 sigma (r) = |\/ n (r)| ,
137 0 up
139 _ __ _ __ _
140 sigma (r) = \/ n (r) . \/ n (r) ,
141 1 up dn
143 _ __ _ 2
144 sigma (r) = |\/ n (r)| .
145 2 dn
146 """
147 nspins = len(n_sg)
148 gradn_svg = gd.empty((nspins, 3))
149 sigma_xg = gd.zeros(nspins * 2 - 1)
150 for v in range(3):
151 for s in range(nspins):
152 grad_v[v](n_sg[s], gradn_svg[s, v])
153 axpy(1.0, gradn_svg[s, v]**2, sigma_xg[2 * s])
154 if nspins == 2:
155 axpy(1.0, gradn_svg[0, v] * gradn_svg[1, v], sigma_xg[1])
156 return sigma_xg, gradn_svg
159def add_gradient_correction(grad_v, gradn_svg, sigma_xg, dedsigma_xg, v_sg):
160 r"""Add gradient correction to potential.
162 ::
164 __ / de(r) __ \
165 v (r) += -2 \/ . | --------- \/ n(r) |
166 xc \ dsigma(r) /
168 Adds arbitrary data to sigma_xg. Be sure to pass a copy if you need
169 sigma_xg after calling this function.
170 """
171 nspins = len(v_sg)
172 # vv_g is a calculation buffer. Its contents will be corrupted.
173 vv_g = sigma_xg[0]
174 for v in range(3):
175 for s in range(nspins):
176 grad_v[v](dedsigma_xg[2 * s] * gradn_svg[s, v], vv_g)
177 vv_g *= 2.0
178 v_sg[s] -= vv_g
179 # axpy(-2.0, vv_g, v_sg[s])
180 if nspins == 2:
181 grad_v[v](dedsigma_xg[1] * gradn_svg[s, v], vv_g)
182 v_sg[1 - s] -= vv_g
183 # axpy(-1.0, vv_g, v_sg[1 - s])
184 # TODO: can the number of gradient evaluations be reduced?
187def gga_vars(gd, grad_v, n_sg):
188 nspins = len(n_sg)
189 sigma_xg, gradn_svg = calculate_sigma(gd, grad_v, n_sg)
190 dedsigma_xg = gd.empty(nspins * 2 - 1)
191 return sigma_xg, dedsigma_xg, gradn_svg
194def get_gradient_ops(gd, nn, xp):
195 return [Gradient(gd, v, n=nn, xp=xp).apply for v in range(3)]
198class GGA(XCFunctional):
199 def __init__(self, kernel, stencil=2):
200 XCFunctional.__init__(self, kernel.name, kernel.type)
201 self.kernel = kernel
202 self.stencil_range = stencil
204 def set_grid_descriptor(self, gd):
205 XCFunctional.set_grid_descriptor(self, gd)
206 self.grad_v = get_gradient_ops(gd, self.stencil_range, self.xp)
208 def todict(self):
209 d = super().todict()
210 d['stencil'] = self.stencil_range
211 return d
213 def get_description(self):
214 return ('{} with {} nearest neighbor stencil'
215 .format(self.name, self.stencil_range))
217 def calculate_impl(self, gd, n_sg, v_sg, e_g):
218 sigma_xg, dedsigma_xg, gradn_svg = gga_vars(gd, self.grad_v, n_sg)
219 self.kernel.calculate(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg)
220 add_gradient_correction(self.grad_v, gradn_svg, sigma_xg,
221 dedsigma_xg, v_sg)
223 def calculate_paw_correction(self, setup, D_sp, dEdD_sp=None,
224 addcoredensity=True, a=None):
225 rcalc = GGARadialCalculator(self.kernel)
226 expansion = GGARadialExpansion(rcalc)
227 return calculate_paw_correction(expansion,
228 setup, D_sp, dEdD_sp,
229 addcoredensity, a)
231 def stress_tensor_contribution(self, n_sg, skip_sum=False):
232 sigma_xg, gradn_svg = calculate_sigma(self.gd, self.grad_v, n_sg)
233 nspins = len(n_sg)
234 dedsigma_xg = self.gd.empty(nspins * 2 - 1)
235 v_sg = self.gd.zeros(nspins)
236 e_g = self.gd.empty()
237 self.kernel.calculate(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg)
239 stress_vv = stress_lda_term(self.gd, e_g, n_sg, v_sg)
240 stress_vv[:] += stress_gga_term(self.gd, sigma_xg, gradn_svg,
241 dedsigma_xg)
243 if not skip_sum:
244 self.gd.comm.sum(stress_vv)
245 return stress_vv
247 def calculate_spherical(self, rgd, n_sg, v_sg, e_g=None):
248 dndr_sg = np.empty_like(n_sg)
249 for n_g, dndr_g in zip(n_sg, dndr_sg):
250 rgd.derivative(n_g, dndr_g)
251 if e_g is None:
252 e_g = rgd.empty()
254 rcalc = GGARadialCalculator(self.kernel)
256 e_g[:], dedn_sg = rcalc(rgd, n_sg[:, np.newaxis],
257 [1.0],
258 dndr_sg[:, np.newaxis],
259 np.zeros((1, 3)), n=None)[:2]
260 v_sg[:] = dedn_sg
261 return rgd.integrate(e_g)
264class PurePythonGGAKernel:
265 def __init__(self, name):
266 self.type = 'GGA'
267 self.name, self.kappa, self.mu, self.beta = pbe_constants(name)
269 def calculate(self, e_g, n_sg, dedn_sg,
270 sigma_xg, dedsigma_xg,
271 tau_sg=None, dedtau_sg=None):
272 e_g[:] = 0.
273 dedsigma_xg[:] = 0.
275 # spin-paired:
276 if len(n_sg) == 1:
277 n = n_sg[0]
278 n[n < 1e-20] = 1e-40
280 # exchange
281 res = gga_x(self.name, 0, n, sigma_xg[0], self.kappa, self.mu)
282 ex, rs, dexdrs, dexda2 = res
283 # correlation
284 res = gga_c(self.name, 0, n, sigma_xg[0], 0, self.beta)
285 ec, rs_, decdrs, decda2, decdzeta = res
287 e_g[:] += n * (ex + ec)
288 dedn_sg[:] += ex + ec - rs * (dexdrs + decdrs) / 3.
289 dedsigma_xg[:] += n * (dexda2 + decda2)
291 # spin-polarized:
292 else:
293 na = 2. * n_sg[0]
294 na[na < 1e-20] = 1e-40
296 nb = 2. * n_sg[1]
297 nb[nb < 1e-20] = 1e-40
299 n = 0.5 * (na + nb)
300 zeta = 0.5 * (na - nb) / n
302 # exchange
303 exa, rsa, dexadrs, dexada2 = gga_x(
304 self.name, 1, na, 4.0 * sigma_xg[0], self.kappa, self.mu)
305 exb, rsb, dexbdrs, dexbda2 = gga_x(
306 self.name, 1, nb, 4.0 * sigma_xg[2], self.kappa, self.mu)
307 a2 = sigma_xg[0] + 2.0 * sigma_xg[1] + sigma_xg[2]
308 # correlation
309 ec, rs, decdrs, decda2, decdzeta = gga_c(
310 self.name, 1, n, a2, zeta, self.beta)
312 e_g[:] += 0.5 * (na * exa + nb * exb) + n * ec
313 dedn_sg[0][:] += (exa + ec - (rsa * dexadrs + rs * decdrs) / 3.0
314 - (zeta - 1.0) * decdzeta)
315 dedn_sg[1][:] += (exb + ec - (rsb * dexbdrs + rs * decdrs) / 3.0
316 - (zeta + 1.0) * decdzeta)
317 dedsigma_xg[0][:] += 2.0 * na * dexada2 + n * decda2
318 dedsigma_xg[1][:] += 2.0 * n * decda2
319 dedsigma_xg[2][:] += 2.0 * nb * dexbda2 + n * decda2
322def pbe_constants(name):
323 if name == 'pyPBE':
324 name = 'PBE'
325 kappa = 0.804
326 mu = 0.2195149727645171
327 beta = 0.06672455060314922
328 elif name in ['pyPBEsol', 'pyzvPBEsol']:
329 name = name[2:]
330 kappa = 0.804
331 mu = 10. / 81.
332 beta = 0.046
333 elif name == 'pyRPBE':
334 name = 'RPBE'
335 kappa = 0.804
336 mu = 0.2195149727645171
337 beta = 0.06672455060314922
338 else:
339 raise NotImplementedError(name)
341 return name, kappa, mu, beta
344# a2 = |grad n|^2
345def gga_x(name, spin, n, a2, kappa, mu, dedmu_g=None):
346 # if dedmu is given, calculate also d(e_x)/d(mu)
347 assert spin in [0, 1]
349 C0I, C1, C2, C3, CC1, CC2, IF2, GAMMA = gga_constants()
350 rs = (C0I / n)**(1 / 3.)
352 # lda part
353 ex = C1 / rs
354 dexdrs = -ex / rs
356 # gga part
357 c = (C2 * rs / n)**2.
358 s2 = a2 * c
360 if name in ['PBE', 'PBEsol', 'zvPBEsol', 'QNA']:
361 x = 1.0 + mu * s2 / kappa
362 Fx = 1.0 + kappa - kappa / x
363 dFxds2 = mu / (x**2.)
364 elif name == 'RPBE':
365 arg = np.maximum(-mu * s2 / kappa, -5.e2)
366 x = np.exp(arg)
367 Fx = 1.0 + kappa * (1.0 - x)
368 dFxds2 = mu * x
369 else:
370 raise NotImplementedError
372 ds2drs = 8.0 * c * a2 / rs
373 dexdrs = dexdrs * Fx + ex * dFxds2 * ds2drs
374 dexda2 = ex * dFxds2 * c
375 if dedmu_g is not None:
376 dedmu_g[:] = ex * s2 / (1 + mu * s2 / kappa)**2
377 ex *= Fx
379 return ex, rs, dexdrs, dexda2
382def gga_c(name, spin, n, a2, zeta, BETA, decdbeta_g=None):
383 # If decdbeta is specified, calculate also d(e_c)/d(beta)
384 assert spin in [0, 1]
385 from gpaw.xc.lda import G
387 C0I, C1, C2, C3, CC1, CC2, IF2, GAMMA = gga_constants()
388 rs = (C0I / n)**(1 / 3.)
390 if name == 'zvPBEsol':
391 zv_a = 1.8
392 zv_o = 9. / 2.
393 zv_x = 1. / 6.
395 # lda part
396 ec, decdrs_0 = G(rs**0.5, 0.031091, 0.21370, 7.5957,
397 3.5876, 1.6382, 0.49294)
399 if spin == 0:
400 decdrs = decdrs_0
401 decdzeta = 0. # dummy
402 else:
403 e1, decdrs_1 = G(rs**0.5, 0.015545, 0.20548, 14.1189,
404 6.1977, 3.3662, 0.62517)
405 alpha, dalphadrs = G(rs**0.5, 0.016887, 0.11125, 10.357,
406 3.6231, 0.88026, 0.49671)
407 alpha *= -1.
408 dalphadrs *= -1.
409 zp = 1.0 + zeta
410 zm = 1.0 - zeta
411 xp = zp**(1 / 3.)
412 xm = zm**(1 / 3.)
413 f = CC1 * (zp * xp + zm * xm - 2.0)
414 f1 = CC2 * (xp - xm)
415 zeta3 = zeta * zeta * zeta
416 zeta4 = zeta * zeta * zeta * zeta
417 x = 1.0 - zeta4
418 decdrs = (decdrs_0 * (1.0 - f * zeta4) +
419 decdrs_1 * f * zeta4 +
420 dalphadrs * f * x * IF2)
421 decdzeta = (4.0 * zeta3 * f * (e1 - ec - alpha * IF2) +
422 f1 * (zeta4 * e1 - zeta4 * ec + x * alpha * IF2))
423 ec += alpha * IF2 * f * x + (e1 - ec) * f * zeta4
425 # gga part
426 n2 = n * n
427 if spin == 1:
428 phi = 0.5 * (xp * xp + xm * xm)
429 phi2 = phi * phi
430 phi3 = phi * phi2
431 t2 = C3 * a2 * rs / (n2 * phi2)
432 y = -ec / (GAMMA * phi3)
433 if name == 'zvPBEsol':
434 u3 = t2**(3. / 2.) * phi3 * (rs / 3.)**(-3. * zv_x)
435 zvarg = -zv_a * u3 * abs(zeta)**zv_o
436 zvf = np.exp(zvarg)
437 else:
438 t2 = C3 * a2 * rs / n2
439 y = -ec / GAMMA
441 x = np.exp(y)
443 A = np.zeros_like(x)
444 indices = np.nonzero(y)
445 if isinstance(BETA, np.ndarray):
446 A[indices] = (BETA[indices] / (GAMMA * (x[indices] - 1.0)))
447 else:
448 A[indices] = (BETA / (GAMMA * (x[indices] - 1.0)))
450 At2 = A * t2
451 nom = 1.0 + At2
452 denom = nom + At2 * At2
453 X = 1.0 + BETA * t2 * nom / (denom * GAMMA)
454 H = GAMMA * np.log(X)
455 tmp = (GAMMA * BETA / (denom * (BETA * t2 * nom + GAMMA * denom)))
456 tmp2 = A * A * x / BETA
457 dAdrs = tmp2 * decdrs
458 if spin == 1:
459 H *= phi3
460 tmp *= phi3
461 dAdrs /= phi3
462 if name == 'zvPBEsol':
463 H_ = H.copy()
464 H *= zvf
465 dHdt2 = (1.0 + 2.0 * At2) * tmp
466 dHdA = -At2 * t2 * t2 * (2.0 + At2) * tmp
467 decdrs += dHdt2 * 7.0 * t2 / rs + dHdA * dAdrs
468 decda2 = dHdt2 * C3 * rs / n2
469 if spin == 1:
470 dphidzeta = np.zeros_like(x)
471 ind1 = np.nonzero(xp)
472 ind2 = np.nonzero(xm)
473 dphidzeta[ind1] += 1.0 / (3.0 * xp[ind1])
474 dphidzeta[ind2] -= 1.0 / (3.0 * xm[ind2])
475 dAdzeta = tmp2 * (decdzeta - 3.0 * ec * dphidzeta / phi) / phi3
476 decdzeta += ((3.0 * H / phi - dHdt2 * 2.0 * t2 / phi) * dphidzeta
477 + dHdA * dAdzeta)
478 decda2 /= phi2
479 if name == 'zvPBEsol':
480 u3_ = (t2**(3. / 2.) *
481 phi3 * rs**(-3. * zv_x - 1.) / 3.**(-3. * zv_x))
482 zvarg_ = -zv_a * u3_ * abs(zeta)**zv_o
483 dzvfdrs = -3. * zv_x * zvf * zvarg_
484 decdrs *= zvf
485 decdrs += H_ * dzvfdrs
487 dt2da2 = C3 * rs / n2
488 assert np.shape(dt2da2) == np.shape(t2)
489 dzvfda2 = dt2da2 * zvf * zvarg * 3. / (2. * t2)
490 decda2 *= zvf
491 decda2 += H_ * dzvfda2
493 dadz = zv_o * abs(zeta)**(zv_o - 1.)
494 dbdphi = 3. * u3 / phi
495 dPdzeta = u3 * dadz + abs(zeta)**(zv_o) * dbdphi * dphidzeta
496 dzvfdzeta = -zv_a * zvf * dPdzeta
497 decdzeta *= zvf
498 decdzeta += H_ * dzvfdzeta
499 ec += H
501 if decdbeta_g is not None:
502 if spin == 0:
503 phi3 = 1.0
504 Y = GAMMA * phi3
505 decdbeta_g[:] = Y / X * t2 / GAMMA
506 decdbeta_g *= ((1 + 2 * At2) / (1 + At2 + At2**2) -
507 (1 + At2) * (At2 + 2 * At2**2) / (1 + At2 + At2**2)**2)
509 return ec, rs, decdrs, decda2, decdzeta
512def gga_constants():
513 from gpaw.xc.lda import lda_constants
514 C0I, C1, CC1, CC2, IF2 = lda_constants()
515 C2 = 0.26053088059892404
516 C3 = 0.10231023756535741
517 GAMMA = 0.0310906908697
519 return C0I, C1, C2, C3, CC1, CC2, IF2, GAMMA