Coverage for gpaw/response/chi0_drude.py: 98%
87 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1from __future__ import annotations
2from time import ctime
4from typing import TYPE_CHECKING
5import numpy as np
6from ase.units import Ha
8from gpaw.response.symmetrize import HeadSymmetryOperators
9from gpaw.response.integrators import Integrand, HilbertTetrahedron, Intraband
10from gpaw.response.chi0_base import Chi0ComponentCalculator
11from gpaw.response.chi0_data import Chi0DrudeData
12from gpaw.response.frequencies import FrequencyGridDescriptor
14if TYPE_CHECKING:
15 from gpaw.response.kpoints import KPointDomainGenerator
18class Chi0DrudeCalculator(Chi0ComponentCalculator):
19 """Class for calculating the plasma frequency contribution to Chi0,
20 that is, the contribution from intraband transitions inside of metallic
21 bands. This corresponds directly to the dielectric function in the Drude
22 model."""
24 def __init__(self, *args, **kwargs):
25 # Serial block distribution
26 super().__init__(*args, nblocks=1, **kwargs)
28 # task: IntegralTask from gpaw.response.integrators
29 # wd: FrequencyDescriptor from gpaw.response.frequencies
30 self.task, self.wd = self.construct_integral_task_and_wd()
32 def calculate(self, wd, rate) -> Chi0DrudeData:
33 """Calculate the Drude dielectric response.
35 Parameters
36 ----------
37 wd : FrequencyDescriptor from gpaw.response.frequencies
38 Frequencies to evaluate the reponse function at.
39 rate : float
40 Plasma frequency decay rate (in eV), corresponding to the
41 imaginary part of the complex frequency.
42 """
43 self.print_info(wd, rate)
45 chi0_drude = Chi0DrudeData.from_frequency_descriptor(wd, rate)
46 self._calculate(chi0_drude)
48 return chi0_drude
50 def _calculate(self, chi0_drude: Chi0DrudeData):
51 """In-place calculation of the Drude dielectric response function,
52 based on the free-space plasma frequency of the intraband transitions.
53 """
54 q_c = [0., 0., 0.]
55 # symmetries: QSymmetries from gpaw.response.symmetry
56 # generator: KPointDomainGenerator from gpaw.response.kpoints
57 # domain: Domain from from gpaw.response.integrators
58 symmetries, generator, domain, prefactor = self.get_integration_domain(
59 q_c=q_c, spins=range(self.gs.nspins))
61 # The plasma frequency integral is special in the way that only
62 # the spectral part is needed
63 integrand = PlasmaFrequencyIntegrand(
64 self, generator, self.gs.gd.cell_cv)
66 # Integrate using temporary array
67 tmp_plasmafreq_wvv = np.zeros((1,) + chi0_drude.vv_shape, complex)
69 # integrator: Integrator from gpaw.response.integrators (or child)
70 self.integrator.integrate(task=self.task,
71 domain=domain, # Integration domain
72 integrand=integrand,
73 wd=self.wd,
74 out_wxx=tmp_plasmafreq_wvv) # Output array
75 tmp_plasmafreq_wvv *= prefactor
77 # Symmetrize the plasma frequency
78 operators = HeadSymmetryOperators(symmetries, self.gs.gd)
79 plasmafreq_vv = tmp_plasmafreq_wvv[0].copy()
80 operators.symmetrize_wvv(plasmafreq_vv[np.newaxis])
82 # Store and print the plasma frequency
83 chi0_drude.plasmafreq_vv += 4 * np.pi * plasmafreq_vv
84 self.context.print('Plasma frequency:', flush=False)
85 self.context.print((chi0_drude.plasmafreq_vv**0.5 * Ha).round(2))
87 # Calculate the Drude dielectric response function from the
88 # free-space plasma frequency
89 # χ_D(ω+iη) = ω_p^2 / (ω+iη)^2
91 # zd: ComplexFrequencyDescriptor from gpaw.response.frequencies
92 assert chi0_drude.zd.upper_half_plane
93 chi0_drude.chi_Zvv += plasmafreq_vv[np.newaxis] \
94 / chi0_drude.zd.hz_z[:, np.newaxis, np.newaxis]**2
96 def construct_integral_task_and_wd(self):
97 if self.integrationmode == 'tetrahedron integration':
98 # Calculate intraband transitions at T=0
99 fermi_level = self.gs.fermi_level
100 wd = FrequencyGridDescriptor([-fermi_level])
101 task = HilbertTetrahedron(self.integrator.blockcomm)
102 else:
103 task = Intraband()
105 # We want to pass None for frequency descriptor, but
106 # if that goes wrong we'll get TypeError which is unhelpful.
107 # This dummy class will give us error messages that allow finding
108 # this spot in the code.
109 class NotAFrequencyDescriptor:
110 pass
112 wd = NotAFrequencyDescriptor()
113 return task, wd
115 def print_info(self, wd, rate):
116 isl = ['',
117 f'{ctime()}',
118 'Calculating drude chi0 with:',
119 f' Number of frequency points:{len(wd)}',
120 f' Plasma frequency decay rate: {rate} eV',
121 '',
122 self.get_gs_info_string(tab=' ')]
124 # context: ResponseContext from gpaw.response.context
125 self.context.print('\n'.join(isl))
128class PlasmaFrequencyIntegrand(Integrand):
129 def __init__(self, chi0drudecalc: Chi0DrudeCalculator,
130 generator: KPointDomainGenerator,
131 cell_cv: np.ndarray):
132 self._drude = chi0drudecalc
133 self.generator = generator
134 self.cell_cv = cell_cv
136 def _band_summation(self):
137 # Intraband response needs only integrate partially unoccupied bands.
138 # gs: ResponseGroundStateAdapter from gpaw.response.groundstate
139 return self._drude.gs.nocc1, self._drude.gs.nocc2
141 def matrix_element(self, point):
142 """NB: In dire need of documentation! XXX."""
143 k_v = point.kpt_c # XXX _v vs _c discrepancy
144 n1, n2 = self._band_summation()
145 k_c = np.dot(self.cell_cv, k_v) / (2 * np.pi)
147 # kptpair_factory: KPointPairFactory from gpaw.response.pair
148 kptpair_factory = self._drude.kptpair_factory
150 K0 = kptpair_factory.gs.kpoints.kptfinder.find(k_c) # XXX
152 kpt1 = kptpair_factory.get_k_point(point.spin, K0, n1, n2)
153 n_n = np.arange(n1, n2)
155 pair_calc = kptpair_factory.pair_calculator(
156 blockcomm=self._drude.blockcomm)
157 vel_nv = pair_calc.intraband_pair_density(kpt1, n_n)
159 if self._drude.integrationmode == 'point integration':
160 f_n = kpt1.f_n
161 width = self._drude.gs.get_occupations_width()
162 if width > 1e-15:
163 dfde_n = - 1. / width * (f_n - f_n**2.0)
164 else:
165 dfde_n = np.zeros_like(f_n)
166 vel_nv *= np.sqrt(-dfde_n[:, np.newaxis])
167 weight = np.sqrt(self.generator.get_kpoint_weight(k_c) /
168 self.generator.how_many_symmetries())
169 vel_nv *= weight
171 return vel_nv
173 def eigenvalues(self, point):
174 """A function that can return the intraband eigenvalues.
176 A method describing the integrand of
177 the response function which gives an output that
178 is compatible with the gpaw k-point integration
179 routines."""
180 n1, n2 = self._band_summation()
181 # gs: ResponseGroundStateAdapter from gpaw.response.groundstate
182 gs = self._drude.gs
183 # kd: KPointDescriptor object from gpaw.kpt_descriptor
184 kd = gs.kd
185 k_v = point.kpt_c # XXX v/c discrepancy
186 # gd: GridDescriptor from gpaw.grid_descriptor
187 k_c = np.dot(self.cell_cv, k_v) / (2 * np.pi)
188 K1 = gs.kpoints.kptfinder.find(k_c)
189 ik = kd.bz2ibz_k[K1]
190 kpt1 = gs.kpt_qs[ik][point.spin]
191 assert gs.kd.comm.size == 1
193 return kpt1.eps_n[n1:n2]