Coverage for gpaw/response/matrix_elements.py: 100%
204 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-12 00:18 +0000
1from __future__ import annotations
2from abc import ABC, abstractmethod
4import numpy as np
6from gpaw.sphere.integrate import spherical_truncation_function_collection
7from gpaw.kpt_descriptor import KPointDescriptor
9from gpaw.response import timer
10from gpaw.response.kspair import KohnShamKPointPair
11from gpaw.response.pair import phase_shifted_fft_indices
12from gpaw.response.site_paw import calculate_site_matrix_element_correction
13from gpaw.response.localft import calculate_LSDA_Wxc, add_LSDA_trans_fxc
14from gpaw.response.site_data import AtomicSiteData
17class MatrixElement(ABC):
18 """Data class for transitions distributed Kohn-Sham matrix elements."""
20 def __init__(self, tblocks):
21 self.tblocks = tblocks
22 self.array = self.zeros()
23 assert self.array.shape[0] == tblocks.blocksize
25 @abstractmethod
26 def zeros(self):
27 """Generate matrix element array with zeros."""
29 @property
30 def local_array_view(self):
31 return self.array[:self.tblocks.nlocal]
33 def get_global_array(self):
34 """Get the global (all gathered) matrix element."""
35 return self.tblocks.all_gather(self.array)
38class MatrixElementCalculator(ABC):
39 r"""Abstract base class for matrix element calculators.
41 In the PAW method, Kohn-Sham matrix elements,
42 ˰
43 A_(nks,n'k's') = <ψ_nks|A|ψ_n'k's'>
45 can be evaluated in the space of pseudo waves using the pseudo operator
46 __ __
47 ˷ ˰ \ \ ˷ ˰ ˷ ˰ ˷ ˷
48 A = A + / / |p_ai>[<φ_ai|A|φ_ai'> - <φ_ai|A|φ_ai'>]<p_ai'|
49 ‾‾ ‾‾
50 a i,i'
52 to which effect,
53 ˷ ˷ ˷
54 A_(nks,n'k's') = <ψ_nks|A|ψ_n'k's'>
56 This is an abstract base class for calculating such matrix elements for a
57 number of band and spin transitions t=(n,s)->(n',s') for a given k-point
58 pair k and k + q:
60 A_kt = A_(nks,n'k+qs')
61 """
63 def __init__(self, gs, context):
64 self.gs = gs
65 self.context = context
67 @timer('Calculate matrix element')
68 def __call__(self, kptpair: KohnShamKPointPair, *args) -> MatrixElement:
69 r"""Calculate the matrix element for all transitions t.
71 The calculation is split into a pseudo contribution and a PAW
72 correction:
73 ˷
74 A_kt = A_kt + ΔA_kt,
76 see [PRB 103, 245110 (2021)] for additional details and references.
77 """
78 matrix_element = self.create_matrix_element(kptpair.tblocks, *args)
79 self.add_pseudo_contribution(kptpair, matrix_element)
80 self.add_paw_correction(kptpair, matrix_element)
81 return matrix_element
83 @abstractmethod
84 def create_matrix_element(self, tblocks, *args) -> MatrixElement:
85 """Return a new MatrixElement instance."""
87 def add_pseudo_contribution(self, kptpair, matrix_element):
88 """Add the pseudo matrix element to an output array.
90 The pseudo matrix element is evaluated on the coarse real-space grid
91 and integrated:
93 ˷ ˷ ˰ ˷
94 A_kt = <ψ_nks|A|ψ_n'k+qs'>
96 / ˷ ˰ ˷
97 = | dr ψ_nks^*(r) A ψ_n'k+qs'(r)
98 /
99 """
100 self._add_pseudo_contribution(*self.extract_pseudo_waves(kptpair),
101 matrix_element=matrix_element)
103 def extract_pseudo_waves(self, kptpair):
104 """Extract the pseudo wave functions for each k-point pair transition.
105 """
106 ikpt1 = kptpair.ikpt1
107 ikpt2 = kptpair.ikpt2
109 # Map the k-points from the irreducible part of the BZ to the BZ
110 # k-point K (up to a reciprocal lattice vector)
111 k1_c = self.gs.ibz2bz[kptpair.K1].map_kpoint()
112 k2_c = self.gs.ibz2bz[kptpair.K2].map_kpoint()
114 # Fourier transform the periodic part of the pseudo waves to the coarse
115 # real-space grid and map them to the BZ k-point K (up to the same
116 # reciprocal lattice vector as above)
117 ut1_hR = self.get_periodic_pseudo_waves(kptpair.K1, ikpt1)
118 ut2_hR = self.get_periodic_pseudo_waves(kptpair.K2, ikpt2)
120 # Fold out the pseudo waves to the transition index
121 ut1_mytR = ut1_hR[ikpt1.h_myt]
122 ut2_mytR = ut2_hR[ikpt2.h_myt]
124 return k1_c, k2_c, ut1_mytR, ut2_mytR
126 def add_paw_correction(self, kptpair, matrix_element):
127 r"""Add the matrix element PAW correction to an output array.
129 The PAW correction is calculated using the projector overlaps of the
130 pseudo waves:
131 __ __
132 \ \ ˷ ˷ ˷ ˷
133 ΔA_kt = / / <ψ_nks|p_ai> ΔA_aii' <p_ai'|ψ_n'k+qs'>
134 ‾‾ ‾‾
135 a i,i'
137 where the PAW correction tensor is calculated on a radial grid inside
138 each augmentation sphere of position R_a, using the atom-centered
139 partial waves φ_ai(r):
140 ˰ ˷ ˰ ˷
141 ΔA_aii' = <φ_ai|A|φ_ai'> - <φ_ai|A|φ_ai'>
143 / ˰
144 = | dr [φ_ai^*(r-R_a) A φ_ai'(r-R_a)
145 / ˷ ˰ ˷
146 - φ_ai^*(r-R_a) A φ_ai'(r-R_a)]
147 """
148 ikpt1 = kptpair.ikpt1
149 ikpt2 = kptpair.ikpt2
151 # Map the projections from the irreducible part of the BZ to the BZ
152 # k-point K
153 P1h = self.gs.ibz2bz[kptpair.K1].map_projections(ikpt1.Ph)
154 P2h = self.gs.ibz2bz[kptpair.K2].map_projections(ikpt2.Ph)
156 # Fold out the projectors to the transition index
157 P1_amyti = ikpt1.projectors_in_transition_index(P1h)
158 P2_amyti = ikpt2.projectors_in_transition_index(P2h)
159 assert P1_amyti.atom_partition.comm.size == \
160 P2_amyti.atom_partition.comm.size == 1, \
161 'We need access to the projections of all atoms'
163 self._add_paw_correction(P1_amyti, P2_amyti, matrix_element)
165 @abstractmethod
166 def _add_pseudo_contribution(self, k1_c, k2_c, ut1_mytR, ut2_mytR,
167 matrix_element: MatrixElement):
168 """Add pseudo contribution based on the pseudo waves in real space."""
170 @abstractmethod
171 def _add_paw_correction(self, P1_amyti, P2_amyti,
172 matrix_element: MatrixElement):
173 """Add paw correction based on the projector overlaps."""
175 def get_periodic_pseudo_waves(self, K, ikpt):
176 """FFT the Kohn-Sham orbitals to real space and map them from the
177 irreducible k-point to the k-point in question."""
178 ut_hR = self.gs.gd.empty(ikpt.nh, self.gs.dtype)
179 for h, psit_G in enumerate(ikpt.psit_hG):
180 ut_hR[h] = self.gs.ibz2bz[K].map_pseudo_wave(
181 self.gs.global_pd.ifft(psit_G, ikpt.ik))
183 return ut_hR
186class PlaneWaveMatrixElement(MatrixElement):
187 def __init__(self, tblocks, qpd):
188 self.qpd = qpd
189 super().__init__(tblocks)
191 def zeros(self):
192 return self.qpd.zeros(self.tblocks.blocksize)
195class PlaneWaveMatrixElementCalculator(MatrixElementCalculator):
196 r"""Abstract base class for calculating plane-wave matrix elements.
198 Calculates the following plane-wave matrix element for a given local
199 functional of the electron (spin-)density f[n](r) = f(n(r)) and k-point
200 pair (k, k + q):
202 f_kt(G+q) = f_(nks,n'k+qs')(G+q) = <ψ_nks| e^-i(G+q)r f(r) |ψ_n'k+qs'>
204 /
205 = | dr e^-i(G+q)r ψ_nks^*(r) ψ_n'k+qs'(r) f(r)
206 /
207 """
209 def __init__(self, gs, context,
210 rshelmax: int = -1,
211 rshewmin: float | None = None):
212 """Construct the PlaneWaveMatrixElementCalculator.
214 Parameters
215 ----------
216 rshelmax : int
217 The maximum index l (l < 6) to use in the expansion of f(r) into
218 real spherical harmonics for the PAW correction.
219 rshewmin : float or None
220 If None, the f(r) will be fully expanded up to the chosen lmax.
221 Given as a float (0 < rshewmin < 1), rshewmin indicates what
222 coefficients to use in the expansion. If any (l,m) coefficient
223 contributes with less than a fraction of rshewmin on average, it
224 will not be included.
225 """
226 super().__init__(gs, context)
228 # Expand local functional in real spherical harmonics around each atom
229 rshe_a = []
230 for a, micro_setup in enumerate(self.gs.micro_setups):
231 rshe, info_string = micro_setup.expand_function(
232 self.add_f, lmax=rshelmax, wmin=rshewmin)
233 self.print_rshe_info(a, info_string)
234 rshe_a.append(rshe)
235 self.rshe_a = rshe_a
237 # PAW correction tensor for a given q_c
238 self._currentq_c = None
239 self._F_aGii = None
241 @abstractmethod
242 def add_f(gd, n_sx, f_x):
243 """Add the local functional f(n(r)) to the f_x output array."""
245 def print_rshe_info(self, a, info_string):
246 """Print information about the functional expansion at atom a."""
247 info_string = f'RSHE of atom {a}:\n' + info_string
248 self.context.print(info_string.replace('\n', '\n ') + '\n')
250 def initialize_paw_corrections(self, qpd):
251 """Initialize the PAW corrections ahead of the actual calculation."""
252 self.get_paw_corrections(qpd)
254 def get_paw_corrections(self, qpd):
255 """Get PAW corrections corresponding to a specific q-vector."""
256 if self._currentq_c is None \
257 or not np.allclose(qpd.q_c, self._currentq_c):
258 with self.context.timer('Initialize PAW corrections'):
259 self._F_aGii = self.gs.matrix_element_paw_corrections(
260 qpd, self.rshe_a)
261 self._currentq_c = qpd.q_c
262 return self._F_aGii
264 @staticmethod
265 def create_matrix_element(tblocks, qpd):
266 return PlaneWaveMatrixElement(tblocks, qpd)
268 @timer('Calculate pseudo matrix element')
269 def _add_pseudo_contribution(self, k1_c, k2_c, ut1_mytR, ut2_mytR,
270 matrix_element: PlaneWaveMatrixElement):
271 r"""Add the pseudo matrix element to the output array.
273 The pseudo matrix element is evaluated on the coarse real-space grid
274 and FFT'ed to reciprocal space,
276 ˷ / ˷ ˷
277 f_kt(G+q) = | dr e^-i(G+q)r ψ_nks^*(r) ψ_n'k+qs'(r) f(r)
278 /V0
279 ˷ ˷
280 = FFT_G[e^-iqr ψ_nks^*(r) ψ_n'k+qs'(r) f(r)]
282 where the Kohn-Sham orbitals are normalized to the unit cell and the
283 functional f[n](r+R)=f(r) is lattice periodic.
284 """
285 qpd = matrix_element.qpd
286 # G: reciprocal space
287 f_mytG = matrix_element.local_array_view
288 # R: real space
289 ft_mytR = self._evaluate_pseudo_matrix_element(ut1_mytR, ut2_mytR)
291 # Get the FFT indices corresponding to the pair density Fourier
292 # transform ˷ ˷
293 # FFT_G[e^(-i[k+q-k']r) u_nks^*(r) u_n'k's'(r)]
294 # This includes a (k,k')-dependent phase, since k2_c only is required
295 # to equal k1_c + qpd.q_c modulo a reciprocal lattice vector.
296 Q_G = phase_shifted_fft_indices(k1_c, k2_c, qpd)
298 # Add the desired plane-wave components of the FFT'ed pseudo matrix
299 # element to the output array
300 for f_G, ft_R in zip(f_mytG, ft_mytR):
301 f_G[:] += qpd.fft(ft_R, 0, Q_G) * self.gs.gd.dv
303 @timer('Evaluate pseudo matrix element')
304 def _evaluate_pseudo_matrix_element(self, ut1_mytR, ut2_mytR):
305 """Evaluate the pseudo matrix element in real-space."""
306 # Evaluate the pseudo pair density ˷ ˷
307 nt_mytR = ut1_mytR.conj() * ut2_mytR # u_nks^*(r) u_n'k's'(r)
309 # Evaluate the local functional f(n(r)) on the coarse real-space grid
310 # NB: Here we assume that f(r) is sufficiently smooth to be represented
311 # on a regular grid (unlike the wave functions).
312 n_sR, gd = self.gs.get_all_electron_density(gridrefinement=1)
313 f_R = gd.zeros()
314 self.add_f(gd, n_sR, f_R)
316 return nt_mytR * f_R[np.newaxis]
318 @timer('Calculate the matrix-element PAW corrections')
319 def _add_paw_correction(self, P1_amyti, P2_amyti,
320 matrix_element: PlaneWaveMatrixElement):
321 r"""Add the matrix-element PAW correction to the output array.
323 The correction is calculated from
324 __ __
325 \ \ ˷ ˷ ˷ ˷
326 Δf_kt(G+q) = / / <ψ_nks|p_ai><p_ai'|ψ_n'k+qs'> F_aii'(G+q)
327 ‾‾ ‾‾
328 a i,i'
330 where the matrix-element PAW correction tensor is given by:
332 /
333 F_aii'(G+q) = | dr e^-i(G+q)r [φ_ai^*(r-R_a) φ_ai'(r-R_a)
334 / ˷ ˷
335 - φ_ai^*(r-R_a) φ_ai'(r-R_a)] f[n](r)
336 """
337 f_mytG = matrix_element.local_array_view
338 F_aGii = self.get_paw_corrections(matrix_element.qpd)
339 for a, F_Gii in enumerate(F_aGii):
340 # Make outer product of the projector overlaps
341 P1ccP2_mytii = P1_amyti[a].conj()[..., np.newaxis] \
342 * P2_amyti[a][:, np.newaxis]
343 # Sum over partial wave indices and add correction to the output
344 f_mytG[:] += np.einsum('tij, Gij -> tG', P1ccP2_mytii, F_Gii)
347class NewPairDensityCalculator(PlaneWaveMatrixElementCalculator):
348 r"""Class for calculating pair densities
350 n_kt(G+q) = n_(nks,n'k+qs')(G+q) = <ψ_nks| e^-i(G+q)r |ψ_n'k+qs'>
351 """
353 def __init__(self, gs, context):
354 super().__init__(gs, context,
355 # Expanding f(r) = 1 in real spherical harmonics only
356 # involves l = 0
357 rshelmax=0)
359 def add_f(self, gd, n_sx, f_x):
360 f_x[:] += 1.
362 def print_rshe_info(self, *args):
363 # The expansion in spherical harmonics is trivial (l = 0), so there is
364 # no need to print anything
365 pass
368class TransversePairPotentialCalculator(PlaneWaveMatrixElementCalculator):
369 r"""Calculator for the transverse magnetic pair potential.
371 The transverse magnetic pair potential is a plane-wave matrix element
372 where the local functional is the transverse LDA kernel:
374 W^⟂_kt(G+q) = W^⟂_(nks,n'k+qs')(G+q)
376 = <ψ_nks| e^-i(G+q)r f_LDA^-+(r) |ψ_n'k+qs'>
377 """
379 def add_f(self, gd, n_sx, f_x):
380 return add_LSDA_trans_fxc(gd, n_sx, f_x, fxc='ALDA')
383class SiteMatrixElement(MatrixElement):
384 def __init__(self, tblocks, q_c, sites):
385 self.q_c = q_c
386 self.sites = sites
387 super().__init__(tblocks)
389 def zeros(self):
390 return np.zeros(
391 (self.tblocks.blocksize, len(self.sites), self.sites.npartitions),
392 dtype=complex)
395class SiteMatrixElementCalculator(MatrixElementCalculator):
396 r"""Class for calculating site matrix elements.
398 The site matrix elements are defined as the expectation value of any local
399 functional of the electron density f[n](r) = f(n(r)), evaluated on a given
400 site a for every site partitioning p. The sites are defined in terms of
401 smooth truncation functions Θ(r∊Ω_ap), interpolating smoothly between unity
402 for positions inside the spherical site volume and zero outside it:
404 f^ap_kt = f^ap_(nks,n'k+qs') = <ψ_nks|Θ(r∊Ω_ap)f(r)|ψ_n'k+qs'>
406 /
407 = | dr Θ(r∊Ω_ap) f(r) ψ_nks^*(r) ψ_n'k+qs'(r)
408 /
410 For details, see [publication in preparation].
411 """
413 def __init__(self, gs, context, sites,
414 rshelmax: int = -1,
415 rshewmin: float | None = None):
416 """Construct the SiteMatrixElementCalculator.
418 Parameters
419 ----------
420 rshelmax : int
421 The maximum index l (l < 6) to use in the expansion of f(r) into
422 real spherical harmonics for the PAW correction.
423 rshewmin : float or None
424 If None, the f(r) will be fully expanded up to the chosen lmax.
425 Given as a float (0 < rshewmin < 1), rshewmin indicates what
426 coefficients to use in the expansion. If any (l,m) coefficient
427 contributes with less than a fraction of rshewmin on average, it
428 will not be included.
429 """
430 super().__init__(gs, context)
431 self.sites = sites
432 self.site_data = AtomicSiteData(self.gs, sites)
434 # Expand local functional in real spherical harmonics around each site
435 rshe_a = []
436 for a, micro_setup in enumerate(self.site_data.micro_setup_a):
437 rshe, info_string = micro_setup.expand_function(
438 self.add_f, lmax=rshelmax, wmin=rshewmin)
439 self.print_rshe_info(a, info_string)
440 rshe_a.append(rshe)
441 self.rshe_a = rshe_a
443 # PAW correction tensor
444 self._F_apii = None
446 @abstractmethod
447 def add_f(self, gd, n_sx, f_x):
448 """Add the local functional f(n(r)) to the f_x output array."""
450 def print_rshe_info(self, a, info_string):
451 """Print information about the expansion at site a."""
452 A = self.sites.A_a[a] # Atomic index of site a
453 info_string = f'RSHE of site {a} (atom {A}):\n' + info_string
454 self.context.print(info_string.replace('\n', '\n ') + '\n')
456 def get_paw_correction_tensor(self):
457 if self._F_apii is None:
458 self._F_apii = self.calculate_paw_correction_tensor()
459 return self._F_apii
461 def calculate_paw_correction_tensor(self):
462 """Calculate the site matrix element correction tensor F_ii'^ap."""
463 F_apii = []
464 for rshe, A, rc_p, lambd_p in zip(
465 self.rshe_a, self.sites.A_a,
466 self.sites.rc_ap, self.site_data.lambd_ap):
467 # Calculate the PAW correction
468 pawdata = self.gs.pawdatasets.by_atom[A]
469 F_apii.append(calculate_site_matrix_element_correction(
470 pawdata, rshe, rc_p, self.site_data.drcut, lambd_p))
472 return F_apii
474 def create_matrix_element(self, tblocks, q_c):
475 return SiteMatrixElement(tblocks, q_c, self.sites)
477 @timer('Calculate pseudo site matrix element')
478 def _add_pseudo_contribution(self, k1_c, k2_c, ut1_mytR, ut2_mytR,
479 matrix_element: SiteMatrixElement):
480 """Add the pseudo site matrix element to the output array.
482 The pseudo matrix element is evaluated on the coarse real-space grid
483 and integrated together with the smooth truncation function,
485 ˷ / ˷ ˷
486 f^ap_kt = | dr Θ(r∊Ω_ap) f(r) ψ_nks^*(r) ψ_n'k+qs'(r)
487 /
489 where the Kohn-Sham orbitals are normalized to the unit cell.
490 """
491 # Construct pseudo waves with Bloch phases
492 r_Rc = np.transpose(self.gs.ibz2bz.r_cR, # scaled grid coordinates
493 (1, 2, 3, 0))
494 psit1_mytR = np.exp(2j * np.pi * r_Rc @ k1_c)[np.newaxis] * ut1_mytR
495 psit2_mytR = np.exp(2j * np.pi * r_Rc @ k2_c)[np.newaxis] * ut2_mytR
496 # Calculate real-space pair densities ñ_kt(r)
497 nt_mytR = psit1_mytR.conj() * psit2_mytR
498 # Evaluate the local functional f(n(r)) on the coarse real-space grid
499 # NB: Here we assume that f(r) is sufficiently smooth to be represented
500 # on a regular grid (unlike the wave functions).
501 n_sR, gd = self.gs.get_all_electron_density(gridrefinement=1)
502 f_R = gd.zeros()
503 self.add_f(gd, n_sR, f_R)
505 # Set up spherical truncation function collection on the coarse
506 # real-space grid with a KPointDescriptor including only the q-point.
507 qd = KPointDescriptor([matrix_element.q_c])
508 stfc = spherical_truncation_function_collection(
509 self.gs.gd, self.site_data.spos_ac,
510 self.sites.rc_ap, self.site_data.drcut, self.site_data.lambd_ap,
511 kd=qd, dtype=complex)
513 # Integrate Θ(r∊Ω_ap) f(r) ñ_kt(r)
514 ntlocal = nt_mytR.shape[0]
515 ft_amytp = {a: np.empty((ntlocal, self.sites.npartitions),
516 dtype=complex)
517 for a in range(len(self.sites))}
518 stfc.integrate(nt_mytR * f_R[np.newaxis], ft_amytp, q=0)
520 # Add integral to output array
521 f_mytap = matrix_element.local_array_view
522 for a in range(len(self.sites)):
523 f_mytap[:, a] += ft_amytp[a]
525 @timer('Calculate site matrix element PAW correction')
526 def _add_paw_correction(self, P1_Amyti, P2_Amyti,
527 matrix_element: SiteMatrixElement):
528 r"""Add the site matrix element PAW correction to the output array.
530 For every site a, we only need a PAW correction for that site itself,
531 __
532 \ ˷ ˷ ˷ ˷
533 Δf^ap_kt = / <ψ_nks|p_ai> F_apii' <p_ai'|ψ_n'k+qs'>
534 ‾‾
535 i,i'
537 where F_apii' is the site matrix element correction tensor.
538 """
539 f_mytap = matrix_element.local_array_view
540 F_apii = self.get_paw_correction_tensor()
541 for a, (A, F_pii) in enumerate(zip(
542 self.sites.A_a, F_apii)):
543 # Make outer product of the projector overlaps
544 P1ccP2_mytii = P1_Amyti[A].conj()[..., np.newaxis] \
545 * P2_Amyti[A][:, np.newaxis]
546 # Sum over partial wave indices and add correction to the output
547 f_mytap[:, a] += np.einsum('tij, pij -> tp', P1ccP2_mytii, F_pii)
550class SitePairDensityCalculator(SiteMatrixElementCalculator):
551 """Class for calculating site pair densities.
553 The site pair density corresponds to a site matrix element with f(r) = 1:
555 n^ap_(nks,n'k+qs') = <ψ_nks|Θ(r∊Ω_ap)|ψ_n'k+qs'>
556 """
558 def __init__(self, gs, context, sites):
559 super().__init__(gs, context, sites,
560 # Expanding f(r) = 1 in real spherical harmonics only
561 # involves l = 0
562 rshelmax=0)
564 def add_f(self, gd, n_sx, f_x):
565 f_x[:] += 1.
567 def print_rshe_info(self, *args):
568 # The expansion in spherical harmonics is trivial (l = 0), so there is
569 # no need to print anything
570 pass
573class SiteZeemanPairEnergyCalculator(SiteMatrixElementCalculator):
574 """Class for calculating site Zeeman pair energies.
576 The site Zeeman pair energy is defined as the site matrix element with
577 f(r) = - W_xc^z(r):
579 E^(Z,ap)_(nks,n'k+qs') = - <ψ_nks|Θ(r∊Ω_ap)W_xc^z(r)|ψ_n'k+qs'>
580 """
581 def add_f(self, gd, n_sx, f_x):
582 f_x[:] += - calculate_LSDA_Wxc(gd, n_sx)
585class SiteSpinPairEnergyCalculator(SiteMatrixElementCalculator):
586 """Class for calculating site spin pair energies.
588 The site spin pair energy is defined as the site matrix element with
589 f(r) = B^xc(r) = - |W_xc^z(r)|:
591 d^(xc,ap)_(nks,n'k+qs') = - <ψ_nks|Θ(r∊Ω_ap)|W_xc^z(r)||ψ_n'k+qs'>
592 """
593 def add_f(self, gd, n_sx, f_x):
594 f_x[:] += - np.abs(calculate_LSDA_Wxc(gd, n_sx))