Coverage for gpaw/xc/mgga.py: 90%
305 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-08 00:17 +0000
1from math import sqrt, pi
3import numpy as np
4from scipy.special import eval_legendre
6from gpaw.xc.gga import (add_gradient_correction, gga_vars,
7 GGARadialExpansion, GGARadialCalculator,
8 get_gradient_ops, stress_gga_term)
9from gpaw.xc.lda import (calculate_paw_correction, stress_integral,
10 stress_lda_term)
11from gpaw.xc.functional import XCFunctional
12from gpaw.sphere.lebedev import weight_n
15class MGGA(XCFunctional):
16 orbital_dependent = True
18 def __init__(self, kernel, stencil=2):
19 """Meta GGA functional."""
20 XCFunctional.__init__(self, kernel.name, kernel.type)
21 self.kernel = kernel
22 self.stencil_range = stencil
23 self.fixed_ke = False
25 def set_grid_descriptor(self, gd):
26 self.grad_v = get_gradient_ops(gd, self.stencil_range, np)
27 XCFunctional.set_grid_descriptor(self, gd)
29 def get_setup_name(self):
30 return 'PBE'
32 # This method exists on GGA class as well. Try to solve this
33 # kind of problem when refactoring MGGAs one day.
34 def get_description(self):
35 return ('{} with {} nearest neighbor stencil'
36 .format(self.name, self.stencil_range))
38 def initialize(self, density, hamiltonian, wfs):
39 self.wfs = wfs
40 self.tauct = density.get_pseudo_core_kinetic_energy_density_lfc()
41 self.tauct_G = None
42 self.dedtaut_sG = None
43 if ((not hasattr(hamiltonian, 'xc_redistributor'))
44 or (hamiltonian.xc_redistributor is None)):
45 self.restrict_and_collect = hamiltonian.restrict_and_collect
46 self.distribute_and_interpolate = \
47 density.distribute_and_interpolate
48 else:
49 def _distribute_and_interpolate(in_xR, out_xR=None):
50 tmp_xR = density.interpolate(in_xR)
51 if hamiltonian.xc_redistributor.enabled:
52 out_xR = hamiltonian.xc_redistributor.distribute(tmp_xR,
53 out_xR)
54 elif out_xR is None:
55 out_xR = tmp_xR
56 else:
57 out_xR[:] = tmp_xR
58 return out_xR
60 def _restrict_and_collect(in_xR, out_xR=None):
61 if hamiltonian.xc_redistributor.enabled:
62 in_xR = hamiltonian.xc_redistributor.collect(in_xR)
63 return hamiltonian.restrict(in_xR, out_xR)
64 self.restrict_and_collect = _restrict_and_collect
65 self.distribute_and_interpolate = _distribute_and_interpolate
67 def set_positions(self, spos_ac):
68 self.tauct.set_positions(spos_ac, self.wfs.atom_partition)
69 if self.tauct_G is None:
70 self.tauct_G = self.wfs.gd.empty()
71 self.tauct_G[:] = 0.0
72 self.tauct.add(self.tauct_G)
74 def calculate_impl(self, gd, n_sg, v_sg, e_g):
75 sigma_xg, dedsigma_xg, gradn_svg = gga_vars(gd, self.grad_v, n_sg)
76 self.process_mgga(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg)
77 add_gradient_correction(self.grad_v, gradn_svg, sigma_xg,
78 dedsigma_xg, v_sg)
80 def fix_kinetic_energy_density(self, taut_sG):
81 self.fixed_ke = True
82 self._taut_gradv_init = False
83 self._fixed_taut_sG = taut_sG.copy()
85 def process_mgga(self, e_g, nt_sg, v_sg, sigma_xg, dedsigma_xg):
86 if self.fixed_ke:
87 taut_sG = self._fixed_taut_sG
88 if not self._taut_gradv_init:
89 self._taut_gradv_init = True
90 # ensure initialization for calculation potential
91 self.wfs.calculate_kinetic_energy_density()
92 else:
93 taut_sG = self.wfs.calculate_kinetic_energy_density()
95 if taut_sG is None:
96 taut_sG = self.wfs.gd.zeros(len(nt_sg))
98 if 0: # taut_sG is None:
99 # Below code disabled because it produces garbage in at least
100 # some cases.
101 #
102 # See https://gitlab.com/gpaw/gpaw/issues/124
103 #
104 # Initialize with von Weizsaecker kinetic energy density:
105 nt0_sg = nt_sg.copy()
106 nt0_sg[nt0_sg < 1e-10] = np.inf
107 taut_sg = sigma_xg[::2] / 8 / nt0_sg
108 nspins = self.wfs.nspins
109 taut_sG = self.wfs.gd.empty(nspins)
110 for taut_G, taut_g in zip(taut_sG, taut_sg):
111 self.restrict_and_collect(taut_g, taut_G)
112 else:
113 taut_sg = np.empty_like(nt_sg)
115 for taut_G, taut_g in zip(taut_sG, taut_sg):
116 taut_G += 1.0 / self.wfs.nspins * self.tauct_G
117 self.distribute_and_interpolate(taut_G, taut_g)
119 # bad = taut_sg < tautW_sg + 1e-11
120 # taut_sg[bad] = tautW_sg[bad]
122 # m = 12.0
123 # taut_sg = (taut_sg**m + (tautW_sg / 2)**m)**(1 / m)
125 dedtaut_sg = np.empty_like(nt_sg)
126 self.kernel.calculate(e_g, nt_sg, v_sg, sigma_xg, dedsigma_xg,
127 taut_sg, dedtaut_sg)
129 self.dedtaut_sG = self.wfs.gd.empty(self.wfs.nspins)
130 self.ekin = 0.0
131 for s in range(self.wfs.nspins):
132 self.restrict_and_collect(dedtaut_sg[s], self.dedtaut_sG[s])
133 self.ekin -= self.wfs.gd.integrate(
134 self.dedtaut_sG[s] * (taut_sG[s] -
135 self.tauct_G / self.wfs.nspins))
137 def stress_mgga_term(self, taut_sg, dedtaut_sg):
138 nspins = len(taut_sg)
139 tau_swG = self.wfs.calculate_kinetic_energy_density_crossterms()
140 tau_cross_g = self.gd.empty()
141 P = 0
142 for taut_g, dedtaut_g in zip(taut_sg, dedtaut_sg):
143 P -= stress_integral(self.gd, taut_g, dedtaut_g)
144 stress_vv = P * np.eye(3)
145 for s in range(nspins):
146 w = 0
147 for v1 in range(3):
148 for v2 in range(v1, 3):
149 self.distribute_and_interpolate(
150 tau_swG[s, w], tau_cross_g)
151 stress_contrib = stress_integral(
152 self.gd, tau_cross_g, dedtaut_sg[s]
153 )
154 stress_vv[v1, v2] -= stress_contrib
155 if v1 != v2:
156 stress_vv[v2, v1] -= stress_contrib
157 w += 1
158 self.dedtaut_sG = self.wfs.gd.empty(self.wfs.nspins)
159 for s in range(self.wfs.nspins):
160 self.restrict_and_collect(dedtaut_sg[s], self.dedtaut_sG[s])
161 return stress_vv
163 def stress_tensor_contribution(self, n_sg, skip_sum=False):
164 sigma_xg, dedsigma_xg, gradn_svg = gga_vars(self.gd, self.grad_v, n_sg)
165 taut_sG = self.wfs.calculate_kinetic_energy_density()
166 if taut_sG is None:
167 taut_sG = self.wfs.gd.zeros(len(n_sg))
168 taut_sg = np.empty_like(n_sg)
169 for taut_G, taut_g in zip(taut_sG, taut_sg):
170 taut_G += 1.0 / self.wfs.nspins * self.tauct_G
171 self.distribute_and_interpolate(taut_G, taut_g)
173 nspins = len(n_sg)
174 dedtaut_sg = np.empty_like(n_sg)
175 v_sg = self.gd.zeros(nspins)
176 e_g = self.gd.empty()
177 self.kernel.calculate(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg,
178 taut_sg, dedtaut_sg)
180 stress_vv = stress_lda_term(self.gd, e_g, n_sg, v_sg)
181 stress_vv[:] += stress_gga_term(self.gd, sigma_xg, gradn_svg,
182 dedsigma_xg)
183 stress_vv[:] += self.stress_mgga_term(taut_sg, dedtaut_sg)
185 if not skip_sum:
186 self.gd.comm.sum(stress_vv)
187 return stress_vv
189 def apply_orbital_dependent_hamiltonian(self, kpt, psit_xG,
190 Htpsit_xG, dH_asp=None):
191 self.wfs.apply_mgga_orbital_dependent_hamiltonian(
192 kpt, psit_xG,
193 Htpsit_xG, dH_asp,
194 self.dedtaut_sG[kpt.s])
196 def calculate_paw_correction(self, setup, D_sp, dEdD_sp=None,
197 addcoredensity=True, a=None):
198 assert not hasattr(self, 'D_sp')
199 self.D_sp = D_sp
200 self.n = 0
201 self.ae = True
202 self.xcc = setup.xc_correction
203 self.dEdD_sp = dEdD_sp
205 if self.xcc is not None and self.xcc.tau_npg is None:
206 self.xcc.tau_npg, self.xcc.taut_npg = self.initialize_kinetic(
207 self.xcc)
209 rcalc = self.create_mgga_radial_calculator()
210 expansion = GGARadialExpansion(rcalc)
211 # The damn thing uses too many 'self' variables to define a clean
212 # integrator object.
213 E = calculate_paw_correction(expansion,
214 setup, D_sp, dEdD_sp,
215 addcoredensity, a)
216 del self.D_sp, self.n, self.ae, self.xcc, self.dEdD_sp
217 return E
219 def create_mgga_radial_calculator(self):
220 class MockKernel:
221 def __init__(self, mgga):
222 self.mgga = mgga
224 def calculate(self, e_g, n_sg, v_sg, sigma_xg, dedsigma_xg):
225 self.mgga.mgga_radial(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg)
227 return GGARadialCalculator(MockKernel(self))
229 def mgga_radial(self, e_g, n_sg, v_sg, sigma_xg, dedsigma_xg):
230 n = self.n
231 nspins = len(n_sg)
232 if self.ae:
233 tau_pg = self.xcc.tau_npg[self.n]
234 tauc_g = self.xcc.tauc_g / (sqrt(4 * pi) * nspins)
235 sign = 1.0
236 else:
237 tau_pg = self.xcc.taut_npg[self.n]
238 tauc_g = self.xcc.tauct_g / (sqrt(4 * pi) * nspins)
239 sign = -1.0
240 tau_sg = np.dot(self.D_sp, tau_pg) + tauc_g
242 if 0: # not self.ae:
243 m = 12
244 for tau_g, n_g, sigma_g in zip(tau_sg, n_sg, sigma_xg[::2]):
245 tauw_g = sigma_g / 8 / n_g
246 tau_g[:] = (tau_g**m + (tauw_g / 2)**m)**(1.0 / m)
247 break
249 dedtau_sg = np.empty_like(tau_sg)
250 self.kernel.calculate(e_g, n_sg, v_sg, sigma_xg, dedsigma_xg,
251 tau_sg, dedtau_sg)
252 if self.dEdD_sp is not None:
253 self.dEdD_sp += (sign * weight_n[self.n] *
254 np.inner(dedtau_sg * self.xcc.rgd.dv_g, tau_pg))
255 assert n == self.n
256 self.n += 1
257 if self.n == len(weight_n):
258 self.n = 0
259 self.ae = False
261 def add_forces(self, F_av):
262 dF_av = self.tauct.dict(derivative=True)
263 self.tauct.derivative(self.dedtaut_sG.sum(0), dF_av)
264 for a, dF_v in dF_av.items():
265 F_av[a] += dF_v[0] / self.wfs.nspins
267 def estimate_memory(self, mem):
268 bytecount = self.wfs.gd.bytecount()
269 mem.subnode('MGGA arrays', (1 + self.wfs.nspins) * bytecount)
271 def initialize_kinetic(self, xcc):
272 nii = xcc.nii
273 nn = len(xcc.rnablaY_nLv)
274 ng = len(xcc.phi_jg[0])
276 tau_npg = np.zeros((nn, nii, ng))
277 taut_npg = np.zeros((nn, nii, ng))
278 create_kinetic(xcc, nn, xcc.phi_jg, tau_npg)
279 create_kinetic(xcc, nn, xcc.phit_jg, taut_npg)
280 return tau_npg, taut_npg
283def create_kinetic(xcc, ny, phi_jg, tau_ypg):
284 r"""Short title here.
286 kinetic expression is::
288 __ __
289 tau_s = 1/2 Sum_{i1,i2} D(s,i1,i2) \/phi_i1 . \/phi_i2 +tauc_s
291 here the orbital dependent part is calculated::
293 __ __
294 \/phi_i1 . \/phi_i2 =
295 __ __
296 \/YL1.\/YL2 phi_j1 phi_j2 +YL1 YL2 dphi_j1 dphi_j2
297 ------ ------
298 dr dr
299 __ __
300 \/YL1.\/YL2 [y] = Sum_c A[L1,c,y] A[L2,c,y] / r**2
302 """
303 nj = len(phi_jg)
304 dphidr_jg = np.zeros(np.shape(phi_jg))
305 for j in range(nj):
306 phi_g = phi_jg[j]
307 xcc.rgd.derivative(phi_g, dphidr_jg[j])
309 # Second term:
310 for y in range(ny):
311 i1 = 0
312 p = 0
313 Y_L = xcc.Y_nL[y]
314 for j1, l1, L1 in xcc.jlL:
315 for j2, l2, L2 in xcc.jlL[i1:]:
316 c = Y_L[L1] * Y_L[L2]
317 temp = c * dphidr_jg[j1] * dphidr_jg[j2]
318 tau_ypg[y, p, :] += temp
319 p += 1
320 i1 += 1
321 # first term
322 for y in range(ny):
323 i1 = 0
324 p = 0
325 rnablaY_Lv = xcc.rnablaY_nLv[y, :xcc.Lmax]
326 Ax_L = rnablaY_Lv[:, 0]
327 Ay_L = rnablaY_Lv[:, 1]
328 Az_L = rnablaY_Lv[:, 2]
329 for j1, l1, L1 in xcc.jlL:
330 for j2, l2, L2 in xcc.jlL[i1:]:
331 temp = (Ax_L[L1] * Ax_L[L2] + Ay_L[L1] * Ay_L[L2] +
332 Az_L[L1] * Az_L[L2])
333 temp *= phi_jg[j1] * phi_jg[j2]
334 temp[1:] /= xcc.rgd.r_g[1:]**2
335 temp[0] = temp[1]
336 tau_ypg[y, p, :] += temp
337 p += 1
338 i1 += 1
339 tau_ypg *= 0.5
342class PurePython2DMGGAKernel:
343 def __init__(self, name, pars=None):
344 self.name = name
345 self.pars = pars
346 self.type = 'MGGA'
347 assert self.pars is not None
349 def calculate(self, e_g, n_sg, dedn_sg,
350 sigma_xg, dedsigma_xg,
351 tau_sg, dedtau_sg):
353 e_g[:] = 0.
354 dedsigma_xg[:] = 0.
355 dedtau_sg[:] = 0.
357 # spin-paired:
358 if len(n_sg) == 1:
359 n = n_sg[0]
360 n[n < 1e-20] = 1e-40
361 sigma = sigma_xg[0]
362 sigma[sigma < 1e-20] = 1e-40
363 tau = tau_sg[0]
364 tau[tau < 1e-20] = 1e-40
366 # exchange
367 e_x = twodexchange(n, sigma, tau, self.pars)
368 e_g[:] += e_x * n
370 # spin-polarized:
371 else:
372 n = n_sg
373 n[n < 1e-20] = 1e-40
374 sigma = sigma_xg
375 sigma[sigma < 1e-20] = 1e-40
376 tau = tau_sg
377 tau[tau < 1e-20] = 1e-40
379 # The spin polarized version is handle using the exact spin scaling
380 # Ex[n1, n2] = (Ex[2*n1] + Ex[2*n2])/2
381 na = 2.0 * n[0]
382 nb = 2.0 * n[1]
384 e2na_x = twodexchange(na, 4. * sigma[0], 2. * tau[0], self.pars)
385 e2nb_x = twodexchange(nb, 4. * sigma[2], 2. * tau[1], self.pars)
386 ea_x = e2na_x * na
387 eb_x = e2nb_x * nb
389 e_g[:] += (ea_x + eb_x) / 2.0
392def twodexchange(n, sigma, tau, pars):
393 # parameters for 2 Legendre polynomials
394 parlen_i = int(pars[0])
395 parlen_j = pars[2 + 2 * parlen_i]
396 assert parlen_i == parlen_j
397 pars_i = pars[1:2 + 2 * parlen_i]
398 pars_j = pars[3 + 2 * parlen_i:]
399 trans_i = pars_i[0]
400 trans_j = pars_j[0]
401 orders_i, coefs_i = np.split(pars_i[1:], 2)
402 orders_j, coefs_j = np.split(pars_j[1:], 2)
403 assert len(coefs_i) == len(orders_i)
404 assert len(coefs_j) == len(orders_j)
405 assert len(orders_i) == len(orders_j)
407 # product Legendre expansion of Fx(s, alpha)
408 e_x_ueg, rs = ueg_x(n)
409 Fx = LegendreFx2(n, rs, sigma, tau,
410 trans_i, orders_i, coefs_i, trans_j, orders_j, coefs_j)
411 return e_x_ueg * Fx
414def LegendreFx2(n, rs, sigma, tau,
415 trans_i, orders_i, coefs_i, trans_j, orders_j, coefs_j):
416 # Legendre polynomial basis expansion in 2D
418 # reduced density gradient in transformation t1(s)
419 C2 = 0.26053088059892404
420 s2 = sigma * (C2 * np.divide(rs, n))**2.
421 x_i = transformation(s2, trans_i)
422 assert x_i.all() >= -1.0 and x_i.all() <= 1.0
424 # kinetic energy density parameter alpha in transformation t2(s)
425 alpha = get_alpha(n, sigma, tau)
426 x_j = transformation(alpha, trans_j)
427 assert x_j.all() >= -1.0 and x_j.all() <= 1.0
429 # product exchange enhancement factor
430 Fx_i = legendre_polynomial(x_i, orders_i, coefs_i)
431 # print(Fx_i);asdf
432 Fx_j = legendre_polynomial(x_j, orders_j, coefs_j)
433 Fx = Fx_i * Fx_j
434 return Fx
437def transformation(x, t):
438 if t > 0:
439 tmp = t + x
440 x = 2.0 * np.divide(x, tmp) - 1.0
441 elif int(t) == -1:
442 tmp1 = (1.0 - x**2.0)**3.0
443 tmp2 = (1.0 + x**3.0 + x**6.0)
444 x = -1.0 * np.divide(tmp1, tmp2)
445 else:
446 raise KeyError('transformation %i unknown!' % t)
447 return x
450def get_alpha(n, sigma, tau):
451 # tau LSDA
452 aux = (3. / 10.) * (3.0 * np.pi * np.pi)**(2. / 3.)
453 tau_lsda = aux * n**(5. / 3.)
455 # von Weisaecker
456 ind = (n != 0.).nonzero()
457 gdms = np.maximum(sigma, 1e-40) # |nabla rho|^2
458 tau_w = np.zeros(np.shape(n))
459 tau_w[ind] = np.maximum(np.divide(gdms[ind], 8.0 * n[ind]), 1e-40)
461 # z and alpha
462 tau_ = np.maximum(tau_w, tau)
463 alpha = np.divide(tau_ - tau_w, tau_lsda)
464 assert alpha.all() >= 0.0
465 return alpha
468def ueg_x(n):
469 C0I = 0.238732414637843
470 C1 = -0.45816529328314287
471 rs = (C0I / n)**(1 / 3.)
472 ex = C1 / rs
473 return ex, rs
476def legendre_polynomial(x, orders, coefs):
477 assert len(orders) == len(coefs) == 1
478 return eval_legendre(orders[0], x) * coefs[0]