Coverage for gpaw/response/df.py: 91%
307 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1from __future__ import annotations
2from abc import ABC, abstractmethod
3from dataclasses import dataclass
4import sys
6import numpy as np
7from ase.units import Hartree
9import gpaw.mpi as mpi
11from gpaw.response.pw_parallelization import Blocks1D
12from gpaw.response.coulomb_kernels import CoulombKernel
13from gpaw.response.dyson import DysonEquation
14from gpaw.response.density_kernels import DensityXCKernel
15from gpaw.response.chi0 import Chi0Calculator, get_frequency_descriptor
16from gpaw.response.chi0_data import Chi0Data
17from gpaw.response.pair import get_gs_and_context
19from typing import TYPE_CHECKING
21if TYPE_CHECKING:
22 from gpaw.response.groundstate import CellDescriptor
23 from gpaw.response.frequencies import FrequencyDescriptor
24 from gpaw.response.qpd import SingleQPWDescriptor
27"""
28On the notation in this module.
30When calculating properties such as the dielectric function, EELS spectrum and
31polarizability there are many inherent subtleties relating to (ir)reducible
32representations and inclusion of local-field effects. For the reciprocal space
33representation of the Coulomb potential, we use the following notation
35v or v(q): The bare Coulomb interaction, 4π/|G+q|²
37V or V(q): The specified Coulomb interaction. Will usually be either the bare
38 interaction or a truncated version hereof.
39ˍ ˍ
40V or V(q): The modified Coulomb interaction. Equal to V(q) for finite
41 reciprocal wave vectors G > 0, but modified to exclude long-range
42 interactions, that is, equal to 0 for G = 0.
43"""
46@dataclass
47class Chi0DysonEquations:
48 chi0: Chi0Data
49 coulomb: CoulombKernel
50 xc_kernel: DensityXCKernel | None
51 cd: CellDescriptor
53 def __post_init__(self):
54 if self.coulomb.truncation is None:
55 self.bare_coulomb = self.coulomb
56 else:
57 self.bare_coulomb = self.coulomb.new(truncation=None)
58 # When inverting the Dyson equation, we distribute frequencies globally
59 blockdist = self.chi0.body.blockdist.new_distributor(nblocks='max')
60 self.wblocks = Blocks1D(blockdist.blockcomm, len(self.chi0.wd))
62 @staticmethod
63 def _normalize(direction):
64 if isinstance(direction, str):
65 d_v = {'x': [1, 0, 0],
66 'y': [0, 1, 0],
67 'z': [0, 0, 1]}[direction]
68 else:
69 d_v = direction
70 d_v = np.asarray(d_v) / np.linalg.norm(d_v)
71 return d_v
73 def get_chi0_wGG(self, direction='x'):
74 chi0 = self.chi0
75 chi0_wGG = chi0.body.get_distributed_frequencies_array().copy()
76 if chi0.qpd.optical_limit:
77 # Project head and wings along the input direction
78 d_v = self._normalize(direction)
79 W_w = self.wblocks.myslice
80 chi0_wGG[:, 0] = np.dot(d_v, chi0.chi0_WxvG[W_w, 0])
81 chi0_wGG[:, :, 0] = np.dot(d_v, chi0.chi0_WxvG[W_w, 1])
82 chi0_wGG[:, 0, 0] = np.dot(d_v, np.dot(chi0.chi0_Wvv[W_w], d_v).T)
83 return chi0_wGG
85 def get_coulomb_scaled_kernel(self, modified=False, Kxc_GG=None):
86 """Get the Hxc kernel rescaled by the bare Coulomb potential v(q).
88 Calculates
89 ˷
90 K(q) = v^(-1/2)(q) K_Hxc(q) v^(-1/2)(q),
92 where v(q) is the bare Coulomb potential and
94 K_Hxc(q) = V(q) + K_xc(q).
96 When using the `modified` flag, the specified Coulomb kernel will be
97 replaced with its modified analogue,
98 ˍ
99 V(q) -> V(q)
100 """
101 qpd = self.chi0.qpd
102 if self.coulomb is self.bare_coulomb:
103 v_G = self.coulomb.V(qpd) # bare Coulomb interaction
104 K_GG = np.eye(len(v_G), dtype=complex)
105 else:
106 v_G = self.bare_coulomb.V(qpd)
107 V_G = self.coulomb.V(qpd)
108 K_GG = np.diag(V_G / v_G)
109 if modified:
110 K_GG[0, 0] = 0.
111 if Kxc_GG is not None:
112 sqrtv_G = v_G**0.5
113 K_GG += Kxc_GG / sqrtv_G / sqrtv_G[:, np.newaxis]
114 return v_G, K_GG
116 @staticmethod
117 def invert_dyson_like_equation(in_wGG, K_GG, reuse_buffer=True):
118 """Generalized Dyson equation invertion.
120 Calculates
122 B(q,ω) = [1 - A(q,ω) K(q)]⁻¹ A(q,ω)
124 while possibly storing the output B(q,ω) in the input A(q,ω) buffer.
125 """
126 if reuse_buffer:
127 out_wGG = in_wGG
128 else:
129 out_wGG = np.zeros_like(in_wGG)
130 for w, in_GG in enumerate(in_wGG):
131 out_wGG[w] = DysonEquation(in_GG, in_GG @ K_GG).invert()
132 return out_wGG
134 def rpa_density_response(self, direction='x', qinf_v=None):
135 """Calculate the RPA susceptibility for (semi-)finite q.
137 Currently this is only used by the QEH code, why we don't support a top
138 level user interface.
139 """
140 # Extract χ₀(q,ω)
141 qpd = self.chi0.qpd
142 chi0_wGG = self.get_chi0_wGG(direction=direction)
143 if qpd.optical_limit:
144 # Restore the q-dependence of the head and wings in the q→0 limit
145 assert qinf_v is not None and np.linalg.norm(qinf_v) > 0.
146 d_v = self._normalize(direction)
147 chi0_wGG[:, 1:, 0] *= np.dot(qinf_v, d_v)
148 chi0_wGG[:, 0, 1:] *= np.dot(qinf_v, d_v)
149 chi0_wGG[:, 0, 0] *= np.dot(qinf_v, d_v)**2
150 # Invert Dyson equation, χ(q,ω) = [1 - χ₀(q,ω) V(q)]⁻¹ χ₀(q,ω)
151 V_GG = self.coulomb.kernel(qpd, q_v=qinf_v)
152 chi_wGG = self.invert_dyson_like_equation(chi0_wGG, V_GG)
153 return qpd, chi_wGG, self.wblocks
155 def inverse_dielectric_function(self, direction='x'):
156 """Calculate v^(1/2) χ v^(1/2), from which ε⁻¹(q,ω) is constructed."""
157 return InverseDielectricFunction.from_chi0_dyson_eqs(
158 self, *self.calculate_vchi_symm(direction=direction))
160 def calculate_vchi_symm(self, direction='x', modified=False):
161 """Calculate v^(1/2) χ v^(1/2).
163 Starting from the TDDFT Dyson equation
165 χ(q,ω) = χ₀(q,ω) + χ₀(q,ω) K_Hxc(q,ω) χ(q,ω), (1)
167 the Coulomb scaled susceptibility,
168 ˷
169 χ(q,ω) = v^(1/2)(q) χ(q,ω) v^(1/2)(q)
171 can be calculated from the Dyson-like equation
172 ˷ ˷ ˷ ˷ ˷
173 χ(q,ω) = χ₀(q,ω) + χ₀(q,ω) K(q,ω) χ(q,ω) (2)
175 where
176 ˷
177 K(q,ω) = v^(-1/2)(q) K_Hxc(q,ω) v^(-1/2)(q).
179 Here v(q) refers to the bare Coulomb potential. It should be emphasized
180 that invertion of (2) rather than (1) is not merely a rescaling
181 excercise. In the optical q → 0 limit, the Coulomb kernel v(q) diverges
182 as 1/|G+q|² while the Kohn-Sham susceptibility χ₀(q,ω) vanishes as
183 |G+q|². Treating v^(1/2)(q) χ₀(q,ω) v^(1/2)(q) as a single variable,
184 the effects of this cancellation can be treated accurately within k.p
185 perturbation theory.
186 """
187 chi0_wGG = self.get_chi0_wGG(direction=direction)
188 Kxc_GG = self.xc_kernel(self.chi0.qpd, chi0_wGG=chi0_wGG) \
189 if self.xc_kernel else None
190 v_G, K_GG = self.get_coulomb_scaled_kernel(
191 modified=modified, Kxc_GG=Kxc_GG)
192 # Calculate v^(1/2)(q) χ₀(q,ω) v^(1/2)(q)
193 sqrtv_G = v_G**0.5
194 vchi0_symm_wGG = chi0_wGG # reuse buffer
195 for w, chi0_GG in enumerate(chi0_wGG):
196 vchi0_symm_wGG[w] = chi0_GG * sqrtv_G * sqrtv_G[:, np.newaxis]
197 # Invert Dyson equation
198 vchi_symm_wGG = self.invert_dyson_like_equation(
199 vchi0_symm_wGG, K_GG, reuse_buffer=False)
200 return vchi0_symm_wGG, vchi_symm_wGG
202 def customized_dielectric_function(self, direction='x'):
203 """Calculate Ε(q,ω) = 1 - V(q) P(q,ω)."""
204 V_GG = self.coulomb.kernel(self.chi0.qpd)
205 P_wGG = self.polarizability_operator(direction=direction)
206 nG = len(V_GG)
207 eps_wGG = P_wGG # reuse buffer
208 for w, P_GG in enumerate(P_wGG):
209 eps_wGG[w] = np.eye(nG) - V_GG @ P_GG
210 return CustomizableDielectricFunction.from_chi0_dyson_eqs(
211 self, eps_wGG)
213 def bare_dielectric_function(self, direction='x'):
214 """Calculate v^(1/2) ̄χ v^(1/2), from which ̄ϵ=1-v ̄χ is constructed.
216 The unscreened susceptibility is given by the Dyson-like equation
217 ˍ ˍ ˍ
218 χ(q,ω) = P(q,ω) + P(q,ω) V(q) χ(q,ω), (3)
220 In the special case of RPA, where P(q,ω) = χ₀(q,ω), one may notice that
221 the Dyson-like equation (3) is exactly identical to the TDDFT Dyson
222 equation (1) when replacing the Hartree-exchange-correlation kernel
223 with the modified Coulomb interaction:
224 ˍ
225 K_Hxc(q) -> V(q).
226 ˍ
227 We may thus reuse that functionality to calculate v^(1/2) χ v^(1/2).
228 """
229 if self.xc_kernel:
230 raise NotImplementedError(
231 'Calculation of the bare dielectric function within TDDFT has '
232 'not yet been implemented. For TDDFT dielectric properties, '
233 'please calculate the inverse dielectric function.')
234 vP_symm_wGG, vchibar_symm_wGG = self.calculate_vchi_symm(
235 direction=direction, modified=True)
236 return BareDielectricFunction.from_chi0_dyson_eqs(
237 self, vP_symm_wGG, vchibar_symm_wGG)
239 def polarizability_operator(self, direction='x'):
240 """Calculate the polarizability operator P(q,ω).
242 Depending on the theory (RPA, TDDFT, MBPT etc.), the polarizability
243 operator is approximated in various ways see e.g.
244 [Rev. Mod. Phys. 74, 601 (2002)].
246 In RPA:
247 P(q,ω) = χ₀(q,ω)
249 In TDDFT:
250 P(q,ω) = [1 - χ₀(q,ω) K_xc(q,ω)]⁻¹ χ₀(q,ω)
251 """
252 chi0_wGG = self.get_chi0_wGG(direction=direction)
253 if not self.xc_kernel: # RPA
254 return chi0_wGG
255 # TDDFT (in adiabatic approximations to the kernel)
256 if self.chi0.qpd.optical_limit:
257 raise NotImplementedError(
258 'Calculation of the TDDFT dielectric function via the '
259 'polarizability operator has not been implemented for the '
260 'optical limit. Please calculate the inverse dielectric '
261 'function instead.')
262 # The TDDFT implementation here is invalid in the optical limit
263 # since the chi0_wGG variable already contains Coulomb effects,
264 # chi0_wGG ~ V^(1/2) χ₀ V^(1/2)/(4π).
265 # Furthermore, a direct evaluation of V(q) P(q,ω) does not seem
266 # sensible, since it does not account for the exact cancellation
267 # of the q-dependences of the two functions.
268 # In principle, one could treat the v(q) P(q,ω) product in
269 # perturbation theory, similar to the χ₀(q,ω) v(q) product in the
270 # Dyson equation for χ, but unless we need to calculate the TDDFT
271 # polarizability using truncated kernels, this isn't really
272 # necessary.
273 Kxc_GG = self.xc_kernel(self.chi0.qpd, chi0_wGG=chi0_wGG)
274 return self.invert_dyson_like_equation(chi0_wGG, Kxc_GG)
277@dataclass
278class DielectricFunctionData(ABC):
279 cd: CellDescriptor
280 qpd: SingleQPWDescriptor
281 wd: FrequencyDescriptor
282 wblocks: Blocks1D
283 coulomb: CoulombKernel
284 bare_coulomb: CoulombKernel
286 @classmethod
287 def from_chi0_dyson_eqs(cls, chi0_dyson_eqs, *args, **kwargs):
288 chi0 = chi0_dyson_eqs.chi0
289 return cls(chi0_dyson_eqs.cd, chi0.qpd,
290 chi0.wd, chi0_dyson_eqs.wblocks,
291 chi0_dyson_eqs.coulomb, chi0_dyson_eqs.bare_coulomb,
292 *args, **kwargs)
294 def _macroscopic_component(self, in_wGG):
295 return self.wblocks.all_gather(in_wGG[:, 0, 0])
297 @property
298 def v_G(self):
299 return self.bare_coulomb.V(self.qpd)
301 @abstractmethod
302 def macroscopic_dielectric_function(self) -> ScalarResponseFunctionSet:
303 """Get the macroscopic dielectric function ε_M(q,ω)."""
305 def dielectric_constant(self):
306 return self.macroscopic_dielectric_function().static_limit.real
308 def eels_spectrum(self):
309 """Get the macroscopic EELS spectrum.
311 The spectrum is defined as
313 1
314 EELS(q,ω) ≡ -Im ε⁻¹(q,ω) = -Im ‾‾‾‾‾‾‾.
315 00 ε (q,ω)
316 M
318 In addition to the many-body spectrum, we also calculate the
319 EELS spectrum in the relevant independent-particle approximation,
320 here defined as
322 1
323 EELS₀(ω) = -Im ‾‾‾‾‾‾‾.
324 IP
325 ε (q,ω)
326 M
327 """
328 _, eps0_W, eps_W = self.macroscopic_dielectric_function().arrays
329 eels0_W = -(1. / eps0_W).imag
330 eels_W = -(1. / eps_W).imag
331 return ScalarResponseFunctionSet(self.wd, eels0_W, eels_W)
333 def polarizability(self):
334 """Get the macroscopic polarizability α_M(q,ω).
336 In the special case where dielectric properties are calculated
337 solely based on the bare Coulomb interaction V(q) = v(q), the
338 macroscopic part of the electronic polarizability α(q,ω) is related
339 directly to the macroscopic dielectric function ε_M(q,ω).
341 In particular, we define the macroscroscopic polarizability so as to
342 make it independent of the nonperiodic cell vector lengths. Namely,
344 α (q,ω) ≡ Λ α (q,ω)
345 M 00
347 where Λ is the nonperiodic hypervolume of the unit cell. Thus,
349 α_M(q,ω) = Λ/(4π) (ϵ_M(q,ω) - 1) = Λ/(4π) (ε_M(q,ω) - 1),
351 where the latter equality only holds in the special case V(q) = v(q).
353 In addition to α_M(q,ω), we calculate also the macroscopic
354 polarizability in the relevant independent-particle approximation by
355 replacing ϵ_M with ε_M^(IP).
356 """
357 if self.coulomb is not self.bare_coulomb:
358 raise ValueError(
359 'When using a truncated Coulomb kernel, the polarizability '
360 'cannot be calculated based on the macroscopic dielectric '
361 'function. Please calculate the bare dielectric function '
362 'instead.')
363 _, eps0_W, eps_W = self.macroscopic_dielectric_function().arrays
364 return self._polarizability(eps0_W, eps_W)
366 def _polarizability(self, eps0_W, eps_W):
367 L = self.cd.nonperiodic_hypervolume
368 alpha0_W = L / (4 * np.pi) * (eps0_W - 1)
369 alpha_W = L / (4 * np.pi) * (eps_W - 1)
370 return ScalarResponseFunctionSet(self.wd, alpha0_W, alpha_W)
373@dataclass
374class InverseDielectricFunction(DielectricFunctionData):
375 """Data class for the inverse dielectric function ε⁻¹(q,ω).
377 The inverse dielectric function characterizes the longitudinal response
379 V (q,ω) = ε⁻¹(q,ω) V (q,ω),
380 tot ext
382 where the induced potential due to the electronic system is given by vχ,
384 ε⁻¹(q,ω) = 1 + v(q) χ(q,ω).
386 In this data class, ε⁻¹ is cast in terms if its symmetrized representation
387 ˷
388 ε⁻¹(q,ω) = v^(-1/2)(q) ε⁻¹(q,ω) v^(1/2)(q),
390 that is, in terms of v^(1/2)(q) χ(q,ω) v^(1/2)(q).
392 Please remark that v(q) here refers to the bare Coulomb potential
393 irregardless of whether χ(q,ω) was determined using a truncated analogue.
394 """
395 vchi0_symm_wGG: np.ndarray # v^(1/2)(q) χ₀(q,ω) v^(1/2)(q)
396 vchi_symm_wGG: np.ndarray
398 def macroscopic_components(self):
399 vchi0_W = self._macroscopic_component(self.vchi0_symm_wGG)
400 vchi_W = self._macroscopic_component(self.vchi_symm_wGG)
401 return vchi0_W, vchi_W
403 def macroscopic_dielectric_function(self):
404 """Get the macroscopic dielectric function ε_M(q,ω).
406 Calculates
408 1
409 ε (q,ω) = ‾‾‾‾‾‾‾‾
410 M ε⁻¹(q,ω)
411 00
413 along with the macroscopic dielectric function in the independent-
414 particle random-phase approximation [Rev. Mod. Phys. 74, 601 (2002)],
416 IPRPA
417 ε (q,ω) = 1 - v(q) χ⁰(q,ω)
418 M 00
420 that is, neglecting local-field and exchange-correlation effects.
421 """
422 vchi0_W, vchi_W = self.macroscopic_components()
423 eps0_W = 1 - vchi0_W
424 eps_W = 1 / (1 + vchi_W)
425 return ScalarResponseFunctionSet(self.wd, eps0_W, eps_W)
427 def dynamic_susceptibility(self):
428 """Get the macroscopic components of χ(q,ω) and χ₀(q,ω)."""
429 vchi0_W, vchi_W = self.macroscopic_components()
430 v0 = self.v_G[0] # Macroscopic Coulomb potential (4π/q²)
431 return ScalarResponseFunctionSet(self.wd, vchi0_W / v0, vchi_W / v0)
434@dataclass
435class CustomizableDielectricFunction(DielectricFunctionData):
436 """Data class for customized dielectric functions Ε(q,ω).
438 Ε(q,ω) is customizable in the sense that bare Coulomb interaction v(q) is
439 replaced by an arbitrary interaction V(q) in the formula for ε(q,ω),
441 Ε(q,ω) = 1 - V(q) P(q,ω),
443 where P is the polarizability operator [Rev. Mod. Phys. 74, 601 (2002)].
444 Thus, for any truncated or otherwise cusomized interaction V(q) ≠ v(q),
445 Ε(q,ω) ≠ ε(q,ω) and Ε⁻¹(q,ω) ≠ ε⁻¹(q,ω).
446 """
447 eps_wGG: np.ndarray
449 def macroscopic_customized_dielectric_function(self):
450 """Get the macroscopic customized dielectric function Ε_M(q,ω).
452 We define the macroscopic customized dielectric function as
454 1
455 Ε (q,ω) = ‾‾‾‾‾‾‾‾,
456 M Ε⁻¹(q,ω)
457 00
459 and calculate also the macroscopic dielectric function in the
460 customizable independent-particle approximation:
462 cIP
463 ε (q,ω) = 1 - VP (q,ω) = Ε (q,ω).
464 M 00 00
465 """
466 eps0_W = self._macroscopic_component(self.eps_wGG)
467 # Invert Ε(q,ω) one frequency at a time to compute Ε_M(q,ω)
468 eps_w = np.zeros((self.wblocks.nlocal,), complex)
469 for w, eps_GG in enumerate(self.eps_wGG):
470 eps_w[w] = 1 / np.linalg.inv(eps_GG)[0, 0]
471 eps_W = self.wblocks.all_gather(eps_w)
472 return ScalarResponseFunctionSet(self.wd, eps0_W, eps_W)
474 def macroscopic_dielectric_function(self):
475 """Get the macroscopic dielectric function ε_M(q,ω).
477 In the special case V(q) = v(q), Ε_M(q,ω) = ε_M(q,ω).
478 """
479 if self.coulomb is not self.bare_coulomb:
480 raise ValueError(
481 'The macroscopic dielectric function is defined in terms of '
482 'the bare Coulomb interaction. To truncate the Hartree '
483 'electron-electron correlations, please calculate the inverse '
484 'dielectric function instead.')
485 return self.macroscopic_customized_dielectric_function()
488@dataclass
489class BareDielectricFunction(DielectricFunctionData):
490 """Data class for the bare (unscreened) dielectric function.
492 The bare dielectric function is defined in terms of the unscreened
493 susceptibility,
494 ˍ ˍ
495 ϵ(q,ω) = 1 - v(q) χ(q,ω),
496 ˍ
497 here represented in terms of v^(1/2)(q) χ(q,ω) v^(1/2)(q).
498 """
499 vP_symm_wGG: np.ndarray # v^(1/2) P v^(1/2)
500 vchibar_symm_wGG: np.ndarray
502 def macroscopic_components(self):
503 vP_W = self._macroscopic_component(self.vP_symm_wGG)
504 vchibar_W = self._macroscopic_component(self.vchibar_symm_wGG)
505 return vP_W, vchibar_W
507 def macroscopic_bare_dielectric_function(self):
508 """Get the macroscopic bare dielectric function ϵ_M(q,ω).
510 Calculates
511 ˍ
512 ϵ_M(q,ω) = ϵ (q,ω)
513 00
515 along with the macroscopic dielectric function in the independent-
516 particle approximation:
518 IP
519 ε (q,ω) = 1 - vP (q,ω).
520 M 00
521 """
522 vP_W, vchibar_W = self.macroscopic_components()
523 eps0_W = 1. - vP_W
524 eps_W = 1. - vchibar_W
525 return ScalarResponseFunctionSet(self.wd, eps0_W, eps_W)
527 def macroscopic_dielectric_function(self):
528 """Get the macroscopic dielectric function ε_M(q,ω).
530 In the special case where the unscreened susceptibility is calculated
531 using the bare Coulomb interaction,
532 ˍ ˍ ˍ
533 χ(q,ω) = P(q,ω) + P(q,ω) v(q) χ(q,ω),
535 it holds identically that [Rev. Mod. Phys. 74, 601 (2002)]:
537 ε_M(q,ω) = ϵ_M(q,ω).
538 """
539 if self.coulomb is not self.bare_coulomb:
540 raise ValueError(
541 'The macroscopic dielectric function cannot be obtained from '
542 'the bare dielectric function calculated based on a truncated '
543 'Coulomb interaction. Please calculate the inverse dielectric '
544 'function instead')
545 return self.macroscopic_bare_dielectric_function()
547 def polarizability(self):
548 """Get the macroscopic polarizability α_M(q,ω).
550 In the most general case, the electronic polarizability can be defined
551 in terms of the unscreened susceptibility
552 ˍ ˍ
553 α(q,ω) ≡ - 1/(4π) v(q) χ(q,ω) = 1/(4π) (ϵ(q,ω) - 1).
555 This is especially useful in systems of reduced dimensionality, where
556 χ(q,ω) needs to be calculated using a truncated Coulomb kernel in order
557 to achieve convergence in a feasible way.
558 """
559 _, eps0_W, eps_W = self.macroscopic_bare_dielectric_function().arrays
560 return self._polarizability(eps0_W, eps_W)
563class DielectricFunctionCalculator:
564 def __init__(self, chi0calc: Chi0Calculator):
565 self.chi0calc = chi0calc
566 self.gs = chi0calc.gs
567 self.context = chi0calc.context
569 self._chi0cache: dict[tuple[str, ...], Chi0Data] = {}
571 def get_chi0(self, q_c: list | np.ndarray) -> Chi0Data:
572 """Get the Kohn-Sham susceptibility χ₀(q,ω) for input wave vector q.
574 Keeps a cache of χ₀ for the latest calculated wave vector, thus
575 allowing for investigation of multiple dielectric properties,
576 Coulomb truncations, xc kernels etc. without recalculating χ₀.
577 """
578 # As cache key, we round off and use a string representation.
579 # Not very elegant, but it should work almost always.
580 q_key = [f'{q:.10f}' for q in q_c]
581 key = tuple(q_key)
582 if key not in self._chi0cache:
583 self._chi0cache.clear()
584 self._chi0cache[key] = self.chi0calc.calculate(q_c)
585 self.context.write_timer()
586 return self._chi0cache[key]
588 def get_chi0_dyson_eqs(self,
589 q_c: list | np.ndarray = [0, 0, 0],
590 truncation: str | None = None,
591 xc: str = 'RPA',
592 **xckwargs
593 ) -> Chi0DysonEquations:
594 """Set up the Dyson equation for χ(q,ω) at given wave vector q.
596 Parameters
597 ----------
598 truncation : str or None
599 Truncation of the Hartree kernel.
600 xc : str
601 Exchange-correlation kernel for LR-TDDFT calculations.
602 If xc == 'RPA', the dielectric response is treated in the random
603 phase approximation.
604 **xckwargs
605 Additional parameters for the chosen xc kernel.
606 """
607 chi0 = self.get_chi0(q_c)
608 coulomb = CoulombKernel.from_gs(self.gs, truncation=truncation)
609 if xc == 'RPA':
610 xc_kernel = None
611 else:
612 xc_kernel = DensityXCKernel.from_functional(
613 self.gs, self.context, functional=xc, **xckwargs)
614 return Chi0DysonEquations(chi0, coulomb, xc_kernel, self.gs.cd)
616 def get_bare_dielectric_function(self, direction='x', **kwargs
617 ) -> BareDielectricFunction:
618 """Calculate the bare dielectric function ̄ϵ(q,ω) = 1 - v(q) ̄χ(q,ω).
620 Here v(q) is the bare Coulomb potential while ̄χ is the unscreened
621 susceptibility calculated based on the modified (and possibly
622 truncated) Coulomb potential.
623 """
624 return self.get_chi0_dyson_eqs(
625 **kwargs).bare_dielectric_function(direction=direction)
627 def get_literal_dielectric_function(self, direction='x', **kwargs
628 ) -> CustomizableDielectricFunction:
629 """Calculate the dielectric function ε(q,ω) = 1 - v(q) P(q,ω)."""
630 return self.get_chi0_dyson_eqs(
631 truncation=None, **kwargs).customized_dielectric_function(
632 direction=direction)
634 def get_customized_dielectric_function(self, direction='x', *,
635 truncation: str, **kwargs
636 ) -> CustomizableDielectricFunction:
637 """Calculate the customized dielectric function Ε(q,ω) = 1 -V(q)P(q,ω).
639 In comparison with the literal dielectric function, the bare Coulomb
640 interaction has here been replaced with a truncated analogue v(q)→V(q).
641 """
642 return self.get_chi0_dyson_eqs(
643 truncation=truncation, **kwargs).customized_dielectric_function(
644 direction=direction)
646 def get_inverse_dielectric_function(self, direction='x', **kwargs
647 ) -> InverseDielectricFunction:
648 """Calculate the inverse dielectric function ε⁻¹(q,ω) = v(q) χ(q,ω).
649 """
650 return self.get_chi0_dyson_eqs(
651 **kwargs).inverse_dielectric_function(direction=direction)
654class DielectricFunction(DielectricFunctionCalculator):
655 """This class defines dielectric function related physical quantities."""
657 def __init__(self, calc, *,
658 frequencies=None,
659 ecut=50,
660 hilbert=True,
661 nbands=None, eta=0.2,
662 intraband=True, nblocks=1, world=mpi.world, txt=sys.stdout,
663 truncation=None,
664 qsymmetry=True,
665 integrationmode='point integration', rate=0.0,
666 eshift: float | None = None):
667 """Creates a DielectricFunction object.
669 calc: str
670 The ground-state calculation file that the linear response
671 calculation is based on.
672 frequencies:
673 Input parameters for frequency_grid.
674 Can be an array of frequencies to evaluate the response function at
675 or dictionary of parameters for build-in nonlinear grid
676 (see :ref:`frequency grid`).
677 ecut: float | dict
678 Plane-wave cut-off or dictionary for anoptional planewave
679 descriptor. See response/qpd.py for details.
680 hilbert: bool
681 Use hilbert transform.
682 nbands: int
683 Number of bands from calculation.
684 eta: float
685 Broadening parameter.
686 intraband: bool
687 Include intraband transitions.
688 world: comm
689 mpi communicator.
690 nblocks: int
691 Split matrices in nblocks blocks and distribute them G-vectors or
692 frequencies over processes.
693 txt: str
694 Output file.
695 truncation: str or None
696 None for no truncation.
697 '2D' for standard analytical truncation scheme.
698 Non-periodic directions are determined from k-point grid
699 integrationmode: str
700 if == 'tetrahedron integration' then tetrahedron
701 integration is performed
702 if == 'point integration' then point integration is used
703 eshift: float
704 Shift unoccupied bands
705 """
706 gs, context = get_gs_and_context(calc, txt, world, timer=None)
707 wd = get_frequency_descriptor(frequencies, gs=gs, nbands=nbands)
709 if integrationmode is None:
710 raise DeprecationWarning(
711 "Please use `integrationmode='point integration'` instead")
712 integrationmode = 'point integration'
714 chi0calc = Chi0Calculator(
715 gs, context, nblocks=nblocks,
716 wd=wd,
717 ecut=ecut, nbands=nbands, eta=eta,
718 hilbert=hilbert,
719 intraband=intraband,
720 qsymmetry=qsymmetry,
721 integrationmode=integrationmode,
722 rate=rate, eshift=eshift
723 )
724 super().__init__(chi0calc)
725 self.truncation = truncation
727 def get_frequencies(self) -> np.ndarray:
728 """Return frequencies (in eV) that the χ is evaluated on."""
729 return self.chi0calc.wd.omega_w * Hartree
731 def get_dynamic_susceptibility(self, *args, xc='ALDA',
732 filename='chiM_w.csv',
733 **kwargs):
734 dynsus = self.get_inverse_dielectric_function(
735 *args, xc=xc, truncation=self.truncation,
736 **kwargs).dynamic_susceptibility()
737 if filename:
738 dynsus.write(filename)
739 return dynsus.unpack()
741 def get_dielectric_function(self, *args, filename='df.csv', **kwargs):
742 """Calculate the dielectric function.
744 Generates a file 'df.csv', unless filename is set to None.
746 Returns
747 -------
748 df_NLFC_w: np.ndarray
749 Dielectric function without local field corrections.
750 df_LFC_w: np.ndarray
751 Dielectric function with local field corrections.
752 """
753 df = self.get_inverse_dielectric_function(
754 *args, truncation=self.truncation,
755 **kwargs).macroscopic_dielectric_function()
756 if filename:
757 df.write(filename)
758 return df.unpack()
760 def get_eels_spectrum(self, *args, filename='eels.csv', **kwargs):
761 """Calculate the macroscopic EELS spectrum.
763 Generates a file 'eels.csv', unless filename is set to None.
765 Returns
766 -------
767 eels0_w: np.ndarray
768 Spectrum in the independent-particle random-phase approximation.
769 eels_w: np.ndarray
770 Fully screened EELS spectrum.
771 """
772 eels = self.get_inverse_dielectric_function(
773 *args, truncation=self.truncation, **kwargs).eels_spectrum()
774 if filename:
775 eels.write(filename)
776 return eels.unpack()
778 def get_polarizability(self, q_c: list | np.ndarray = [0, 0, 0],
779 direction='x', filename='polarizability.csv',
780 **xckwargs):
781 """Calculate the macroscopic polarizability.
783 Generate a file 'polarizability.csv', unless filename is set to None.
785 Returns:
786 --------
787 alpha0_w: np.ndarray
788 Polarizability calculated without local-field corrections
789 alpha_w: np.ndarray
790 Polarizability calculated with local-field corrections.
791 """
792 chi0_dyson_eqs = self.get_chi0_dyson_eqs(
793 q_c, truncation=self.truncation, **xckwargs)
794 if self.truncation:
795 # eps: BareDielectricFunction
796 method = chi0_dyson_eqs.bare_dielectric_function
797 else:
798 # eps: CustomizableDielectricFunction
799 method = chi0_dyson_eqs.customized_dielectric_function
800 eps = method(direction=direction)
801 pol = eps.polarizability()
802 if filename:
803 pol.write(filename)
804 return pol.unpack()
806 def get_macroscopic_dielectric_constant(self, xc='RPA', direction='x'):
807 """Calculate the macroscopic dielectric constant.
809 The macroscopic dielectric constant is defined as the real part of the
810 dielectric function in the static limit.
812 Returns:
813 --------
814 eps0: float
815 Dielectric constant without local field corrections.
816 eps: float
817 Dielectric constant with local field correction. (RPA, ALDA)
818 """
819 return self.get_inverse_dielectric_function(
820 xc=xc, direction=direction).dielectric_constant()
823# ----- Serialized dataclasses and IO ----- #
826@dataclass
827class ScalarResponseFunctionSet:
828 """A set of scalar response functions rf₀(ω) and rf(ω)."""
829 wd: FrequencyDescriptor
830 rf0_w: np.ndarray
831 rf_w: np.ndarray
833 @property
834 def arrays(self):
835 return self.wd.omega_w * Hartree, self.rf0_w, self.rf_w
837 def unpack(self):
838 # Legacy feature to support old DielectricFunction output format
839 # ... to be deprecated ...
840 return self.rf0_w, self.rf_w
842 def write(self, filename):
843 if mpi.rank == 0:
844 write_response_function(filename, *self.arrays)
846 @property
847 def static_limit(self):
848 """Return the value of the response functions in the static limit."""
849 w0 = np.argmin(np.abs(self.wd.omega_w))
850 assert abs(self.wd.omega_w[w0]) < 1e-8
851 return np.array([self.rf0_w[w0], self.rf_w[w0]])
854def write_response_function(filename, omega_w, rf0_w, rf_w):
855 with open(filename, 'w') as fd:
856 for omega, rf0, rf in zip(omega_w, rf0_w, rf_w):
857 if rf0_w.dtype == complex:
858 print('%.6f, %.6f, %.6f, %.6f, %.6f' %
859 (omega, rf0.real, rf0.imag, rf.real, rf.imag),
860 file=fd)
861 else:
862 print(f'{omega:.6f}, {rf0:.6f}, {rf:.6f}', file=fd)
865def read_response_function(filename):
866 """Read a stored response function file"""
867 d = np.loadtxt(filename, delimiter=',')
868 omega_w = np.array(d[:, 0], float)
870 if d.shape[1] == 3:
871 # Real response function
872 rf0_w = np.array(d[:, 1], float)
873 rf_w = np.array(d[:, 2], float)
874 elif d.shape[1] == 5:
875 rf0_w = np.array(d[:, 1], complex)
876 rf0_w.imag = d[:, 2]
877 rf_w = np.array(d[:, 3], complex)
878 rf_w.imag = d[:, 4]
879 else:
880 raise ValueError(f'Unexpected array dimension {d.shape}')
882 return omega_w, rf0_w, rf_w