Coverage for gpaw/lcaotddft/qed.py: 91%
128 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1from __future__ import annotations
3import warnings
4from abc import ABC, abstractmethod
6import numpy as np
7from ase.units import alpha, Hartree, Bohr
8from gpaw.external import ConstantElectricField
9from gpaw.lcaotddft.hamiltonian import KickHamiltonian
12def create_environment(environment='waveguide', **kwargs):
13 if isinstance(environment, Environment):
14 assert not kwargs, 'kwargs must not be given here'
15 return environment
17 if isinstance(environment, dict):
18 assert not kwargs, 'please do not give kwargs together with dict'
19 kwargs = dict(environment)
20 environment = kwargs.pop('environment', 'waveguide')
22 if environment == 'waveguide':
23 return WaveguideEnvironment(**kwargs)
25 raise ValueError(f'Unknown environment {environment}')
28def forward_finite_difference(coefficients: list[int],
29 data_tx: np.ndarray):
30 """ Calculate derivatives by forward difference.
32 The data corresponding to the last time is repeated several times,
33 in order to obtain a result of the same shape as the original array.
35 Parameters
36 ----------
37 coefficients
38 List of coefficients.
39 Check e.g. https://en.wikipedia.org/wiki/Finite_difference_coefficient
40 data_tx
41 Array of data, the first axis being the time axis.
43 """
44 N = len(data_tx)
46 padded_shape = (N + len(coefficients) - 1, ) + data_tx.shape[1:]
47 padded_data_tx = np.zeros(padded_shape, dtype=data_tx.dtype)
48 padded_data_tx[:len(coefficients)] = data_tx[:1]
49 padded_data_tx[len(coefficients):] = data_tx[1:]
51 result_tx = np.zeros_like(data_tx)
53 for start, coefficient in enumerate(coefficients):
54 result_tx[:] += coefficient * padded_data_tx[start:start + N]
56 return result_tx
59def calculate_first_derivative(timestep, data_tx):
60 """ Calculate the first derivative by first order forward difference.
62 The data corresponding to the last time is repeated several times,
63 in order to obtain a result of the same shape as the original array.
65 Parameters
66 ----------
67 timestep
68 Time step
69 data_tx
70 Array of data, the first axis being the time axis.
72 Returns
73 -------
74 First derivative, array same shape as :attr:`data_tx`.
75 """
76 if len(data_tx) == 1 and timestep is None:
77 return np.zeros_like(data_tx)
79 # Coefficients from oldest times to newest times
80 coefficients = [-1, 1]
81 return forward_finite_difference(coefficients, data_tx) / timestep
84def calculate_third_derivative(timestep, data_tx):
85 """ Calculate the third derivative by first order forward difference.
87 The data corresponding to the last time is repeated several times,
88 in order to obtain a result of the same shape as the original array.
90 Parameters
91 ----------
92 timestep
93 Time step
94 data_tx
95 Array of data, the first axis being the time axis.
97 Returns
98 -------
99 Third derivative, array same shape as :attr:`data_tx`.
100 """
101 if len(data_tx) == 1 and timestep is None:
102 return np.zeros_like(data_tx)
104 # Coefficients from oldest times to newest times
105 coefficients = [-1, 3, -3, 1]
106 # For a second order accuracy, the coefficients would be
107 # [-5/2, 9, -12, 7, -3/2]
108 return forward_finite_difference(coefficients, data_tx) / timestep ** 3
111class RRemission:
112 r"""
113 Radiation-reaction potential according to Schaefer et al.
114 [https://doi.org/10.1103/PhysRevLett.128.156402] and
115 Schaefer [https://doi.org/10.48550/arXiv.2204.01602].
116 The potential accounts for the friction
117 forces acting on the radiating system of oscillating charges
118 emitting into a single dimension. A more elegant
119 formulation would use the current instead of the dipole.
120 Please contact christian.schaefer.physics@gmail.com if any problems
121 should appear or you would like to consider more complex systems.
122 Big thanks to Tuomas Rossi and Jakub Fojt for their help.
124 Parameters
125 ----------
126 environment
127 Environment, or dictionary with environment parameters.
129 Notes
130 -----
131 A legacy form for constructing the RRemission object is supported
132 for backwards compatibility. In this form, two parameters are given.
134 quantization_plane: float
135 value of :math:`quantization_plane` in units of Å^2
136 cavity_polarization: array
137 value of :math:`cavity_polarization`; dimensionless (directional)
138 """
140 def __init__(self, *args):
141 if len(args) == 1:
142 self.environment = create_environment(args[0])
143 elif len(args) == 2:
144 # Legacy syntax
145 self.environment = create_environment(
146 'waveguide',
147 quantization_plane=args[0],
148 cavity_polarization=args[1])
149 warnings.warn(
150 "Use RRemission({'environment': 'waveguide', "
151 "'quantization_plane': quantization_plane, "
152 "'cavity_polarization': cavity_polarization})) instead of "
153 'RRemission(quantization_plane, cavity_polarization).',
154 FutureWarning)
155 else:
156 raise TypeError('Please provide only one argument '
157 '(two for legacy syntax)')
159 # Recorded dipole moment over time
160 # The entries all correspond to unique times
161 # If there is a kick, the dipole after the kick is recorded
162 # but not the dipole just before the kick
163 self.dipole_tv = np.zeros((0, 3))
164 # Time of last record
165 self.time = 0
167 # Time step in TDDFT
168 self.timestep = None
170 def read(self, reader):
171 self.dipole_tv = reader.recorded_dipole / Bohr
172 self.timestep = reader.timestep
174 def write(self, writer):
175 writer.write(recorded_dipole=self.dipole_tv * Bohr)
176 writer.write(timestep=self.timestep)
177 self.environment.write(writer)
179 @classmethod
180 def from_reader(cls, reader):
181 parameters = reader.parameters.asdict()
183 rremission = cls(parameters)
184 rremission.read(reader)
186 return rremission
188 def initialize(self, paw):
189 self.density = paw.density
190 self.wfs = paw.wfs
191 self.hamiltonian = paw.hamiltonian
192 self.time = paw.time
193 dipole_v = paw.density.calculate_dipole_moment()
194 if len(self.dipole_tv) == 0:
195 # If we are not reading from a restart file, this will be empty
196 self.record_dipole(dipole_v)
198 # Set up three external fields summing to the polarization cavity
199 strength_v = self.environment.cavity_polarization * Hartree / Bohr
200 assert len(strength_v) == 3
201 ext_x = ConstantElectricField(strength_v[0], [1, 0, 0])
202 ext_y = ConstantElectricField(strength_v[1], [0, 1, 0])
203 ext_z = ConstantElectricField(strength_v[2], [0, 0, 1])
204 ext_i = [ext_x, ext_y, ext_z]
206 # Set up the dipole matrix
207 get_matrix = self.wfs.eigensolver.calculate_hamiltonian_matrix
208 self.V_iuMM = []
209 for ext in ext_i:
210 V_uMM = []
211 hamiltonian = KickHamiltonian(self.hamiltonian, self.density, ext)
212 for kpt in self.wfs.kpt_u:
213 V_MM = get_matrix(hamiltonian, self.wfs, kpt,
214 add_kinetic=False, root=-1)
215 V_uMM.append(V_MM)
216 self.V_iuMM.append(V_uMM)
218 def record_dipole(self, dipole_v):
219 """ Record the dipole for the next time step.
221 Parameters
222 ----------
223 dipole_v
224 The new dipole moment.
225 """
226 self.dipole_tv = np.vstack((self.dipole_tv, dipole_v))
228 def update_dipole(self, time, dipole_v):
229 """ Record the new dipole
231 If the time has changed, record the new dipole in a new entry.
233 If not, then it either means that we are iteratively
234 performing propagator steps, or that a kick has been done.
235 Either way, we overwrite the last entry in the record.
237 Parameters
238 ----------
239 time
240 The new time
241 dipole_v
242 The new dipole moment.
243 """
244 if time > self.time or len(self.dipole_tv) == 0:
245 self.record_dipole(dipole_v)
246 timestep = time - self.time
247 if (self.timestep is not None and
248 not np.isclose(self.timestep, timestep)):
249 raise NotImplementedError(
250 'Variable time step in TDDFT not supported'
251 f'{self.timestep} != {timestep}')
252 self.timestep = timestep
253 self.time = time
254 else:
255 self.dipole_tv[-1] = dipole_v
257 def get_MM(self, field_v, kpt):
258 """ Get potential matrix.
260 Parameters
261 ----------
262 field_v
263 The radiation reaction field.
264 kpt
265 kpoint
267 Returns
268 -------
269 The radiation reaction potential matrix.
270 """
271 kpt_rank, q = self.wfs.kd.get_rank_and_index(kpt.k)
272 u = q * self.wfs.nspins + kpt.s
274 Ni = len(self.V_iuMM)
275 Vrr_MM = -field_v[0] * self.V_iuMM[0][u]
276 for i in range(1, Ni):
277 Vrr_MM -= field_v[i] * self.V_iuMM[i][u]
279 return Vrr_MM
281 def vradiationreaction(self, kpt, time):
282 self.update_dipole(time, self.density.calculate_dipole_moment())
284 field_v = self.environment.radiation_reaction_field(
285 self.timestep, self.dipole_tv)
287 Vrr_MM = self.get_MM(field_v, kpt)
288 return Vrr_MM
291class Environment(ABC):
293 def write(self, writer):
294 writer.child('parameters').write(**self.todict())
296 @abstractmethod
297 def radiation_reaction_field(self, timestep, dipole_tv):
298 raise NotImplementedError
300 @abstractmethod
301 def todict(self):
302 raise NotImplementedError
305class WaveguideEnvironment(Environment):
307 r""" The radiation reaction potential in a 1D waveguide is
309 .. math::
311 -\frac{4 \pi \alpha}{A}
312 \boldsymbol{\varepsilon_c} \cdot \partial_t \boldsymbol{R}(t)
313 \boldsymbol{\varepsilon_c} \cdot \hat{\boldsymbol{r}}
315 Parameters
316 ----------
317 quantization_plane
318 Quantization plane :math:`A` in units of Å^2.
319 cavity_polarization
320 The polarization of the cavity :math:`\boldsymbol{\varepsilon_c}`.
321 Direction vector.
322 """
324 def __init__(self, quantization_plane, cavity_polarization):
325 self.quantization_plane = quantization_plane / Bohr**2
326 self.cavity_polarization = np.array(cavity_polarization)
328 def radiation_reaction_field(self, timestep, dipole_tv):
329 d1_dipole_v = calculate_first_derivative(timestep, dipole_tv)[-1]
331 field_v = (4.0 * np.pi * alpha / self.quantization_plane
332 * (d1_dipole_v @ self.cavity_polarization)
333 * self.cavity_polarization)
335 return field_v
337 def todict(self):
338 return {'environment': 'waveguide',
339 'quantization_plane': self.quantization_plane * Bohr**2,
340 'cavity_polarization': list(self.cavity_polarization)}