Coverage for gpaw/tddft/__init__.py: 79%
273 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1"""This module implements a class for (true) time-dependent density
2functional theory calculations.
4"""
5import time
6import warnings
7from math import log
9import numpy as np
11from gpaw.calculator import GPAW
12from gpaw.mixer import DummyMixer
13from gpaw.preconditioner import Preconditioner
14from gpaw.tddft.units import (attosec_to_autime, autime_to_attosec,
15 aufrequency_to_eV)
16from gpaw.tddft.utils import MultiBlas
17from gpaw.tddft.solvers import create_solver
18from gpaw.tddft.propagators import \
19 create_propagator, \
20 AbsorptionKick
21from gpaw.tddft.tdopers import \
22 TimeDependentHamiltonian, \
23 TimeDependentOverlap, \
24 TimeDependentWaveFunctions, \
25 TimeDependentDensity, \
26 AbsorptionKickHamiltonian
27from gpaw.wavefunctions.fd import FD
29from gpaw.tddft.spectrum import photoabsorption_spectrum
30from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
31from gpaw.lcaotddft.magneticmomentwriter import MagneticMomentWriter
32from gpaw.lcaotddft.restartfilewriter import RestartFileWriter
35__all__ = ['TDDFT', 'photoabsorption_spectrum',
36 'DipoleMomentWriter', 'MagneticMomentWriter',
37 'RestartFileWriter']
40# T^-1
41# Bad preconditioner
42class KineticEnergyPreconditioner:
43 def __init__(self, gd, kin, dtype):
44 self.preconditioner = Preconditioner(gd, kin, dtype)
45 self.preconditioner.allocate()
47 def apply(self, kpt, psi, psin):
48 for i in range(len(psi)):
49 psin[i][:] = self.preconditioner(psi[i], kpt.phase_cd, None, None)
52# S^-1
53class InverseOverlapPreconditioner:
54 """Preconditioner for TDDFT."""
55 def __init__(self, overlap):
56 self.overlap = overlap
58 def apply(self, kpt, psi, psin):
59 self.overlap.apply_inverse(psi, psin, kpt)
60# ^^^^^^^^^^
63class FDTDDFTMode(FD):
64 def __call__(self, *args, **kwargs):
65 reuse_wfs_method = kwargs.pop('reuse_wfs_method', None)
66 assert reuse_wfs_method is None
67 return TimeDependentWaveFunctions(self.nn, *args, **kwargs)
70class TDDFT(GPAW):
71 """Time-dependent density functional theory calculation based on GPAW.
73 This class is the core class of the time-dependent density functional
74 theory implementation and is the only class which a user has to use.
75 """
77 def __init__(self, filename: str, *,
78 td_potential: object = None,
79 calculate_energy: bool = True,
80 propagator: dict = None,
81 solver: dict = None,
82 tolerance: float = None, # deprecated
83 parallel: dict = None,
84 communicator: object = None,
85 txt: str = '-'):
86 """
87 Parameters
88 ----------
89 filename
90 File containing ground state or time-dependent state to propagate.
91 td_potential
92 Function class for the time-dependent potential. Must have a method
93 ``strength(time)`` which returns the strength of
94 the linear potential to each direction as a vector of three floats.
95 calculate_energy
96 Whether to calculate energy during propagation.
97 propagator
98 Time propagator for the Kohn-Sham wavefunctions.
99 solver
100 The iterative linear equations solver for propagator.
101 tolerance
102 Deprecated. Do not use this, but use solver dictionary instead.
103 Tolerance for the linear solver.
104 parallel
105 Parallelization options
106 communicator
107 MPI communicator
108 txt
109 Text output
110 """
111 # Default values
112 if propagator is None:
113 propagator = dict(name='SICN')
114 if solver is None:
115 solver = dict(name='CSCG', tolerance=1e-8)
117 assert filename is not None
119 # For communicating with observers
120 self.action = ''
122 self.time = 0.0
123 self.kick_strength = np.array([0.0, 0.0, 0.0], dtype=float)
124 self.kick_gauge = ''
125 self.niter = 0
126 self.dm_file = None # XXX remove and use observer instead
128 # Parallelization dictionary should default to strided bands
129 # but it does not work XXX
130 # self.default_parallel = GPAW.default_parallel.copy()
131 # self.default_parallel['stridebands'] = True
133 self.default_parameters = GPAW.default_parameters.copy()
134 self.default_parameters['mixer'] = DummyMixer()
135 self.default_parameters['experimental']['reuse_wfs_method'] = None
137 # NB: TDDFT restart files contain additional information which
138 # will override the initial settings for time/kick/niter.
139 GPAW.__init__(self, filename, parallel=parallel,
140 communicator=communicator, txt=txt)
141 if len(self.symmetry.op_scc) > 1:
142 raise ValueError('Symmetries are not allowed for TDDFT. '
143 'Run the ground state calculation with '
144 'symmetry={"point_group": False}.')
146 assert isinstance(self.wfs, TimeDependentWaveFunctions)
147 assert isinstance(self.wfs.overlap, TimeDependentOverlap)
149 self.set_positions()
151 # Don't be too strict
152 self.density.charge_eps = 1e-5
154 self.rank = self.wfs.world.rank
156 self.calculate_energy = calculate_energy
157 if self.hamiltonian.xc.name.startswith('GLLB'):
158 self.log('GLLB model potential. Not updating energy.')
159 self.calculate_energy = False
161 # Time-dependent variables and operators
162 self.td_hamiltonian = TimeDependentHamiltonian(self.wfs, self.spos_ac,
163 self.hamiltonian,
164 td_potential)
165 self.td_density = TimeDependentDensity(self)
167 # Solver for linear equations
168 if isinstance(solver, str):
169 solver = dict(name=solver)
170 if tolerance is not None:
171 warnings.warn(
172 "Please specify the solver tolerance using dictionary "
173 "solver={'name': name, 'tolerance': tolerance}. "
174 "Confirm the used tolerance from the output file. "
175 "The old syntax will throw an error in the future.",
176 FutureWarning)
177 solver.update(tolerance=tolerance)
178 self.solver = create_solver(solver)
180 # Preconditioner
181 # No preconditioner as none good found
182 self.preconditioner = None # TODO! check out SSOR preconditioning
183 # self.preconditioner = InverseOverlapPreconditioner(self.overlap)
184 # self.preconditioner = KineticEnergyPreconditioner(
185 # self.wfs.gd, self.td_hamiltonian.hamiltonian.kin, complex)
187 # Time propagator
188 if isinstance(propagator, str):
189 propagator = dict(name=propagator)
190 self.propagator = create_propagator(propagator)
192 self.hpsit = None
193 self.eps_tmp = None
194 self.mblas = MultiBlas(self.wfs.gd)
196 self.tddft_initialized = False
198 def tddft_init(self):
199 if self.tddft_initialized:
200 return
202 self.log('')
203 self.log('')
204 self.log('------------------------------------------')
205 self.log(' Time-propagation TDDFT ')
206 self.log('------------------------------------------')
207 self.log('')
209 self.log('Charge epsilon:', self.density.charge_eps)
211 # Density mixer
212 self.td_density.density.set_mixer(DummyMixer())
214 # Solver
215 self.solver.initialize(self.wfs.gd, self.timer)
216 self.log('Solver:', self.solver.todict())
218 # Preconditioner
219 self.log('Preconditioner:', self.preconditioner)
221 # Propagator
222 self.propagator.initialize(self.td_density, self.td_hamiltonian,
223 self.wfs.overlap, self.solver,
224 self.preconditioner,
225 self.wfs.gd, self.timer)
226 self.log('Propagator:', self.propagator.todict())
228 # Parallelization
229 wfs = self.wfs
230 if self.rank == 0:
231 if wfs.kd.comm.size > 1:
232 if wfs.nspins == 2:
233 self.log('Parallelization Over Spin')
235 if wfs.gd.comm.size > 1:
236 self.log('Using Domain Decomposition: %d x %d x %d' %
237 tuple(wfs.gd.parsize_c))
239 if wfs.bd.comm.size > 1:
240 self.log('Parallelization Over bands on %d Processors' %
241 wfs.bd.comm.size)
242 self.log('States per processor =', wfs.bd.mynbands)
244 # Restarting an FDTD run generates hamiltonian.fdtd_poisson, which
245 # now overwrites hamiltonian.poisson
246 if hasattr(self.hamiltonian, 'fdtd_poisson'):
247 self.hamiltonian.poisson = self.hamiltonian.fdtd_poisson
248 self.hamiltonian.poisson.set_grid_descriptor(self.density.finegd)
250 # For electrodynamics mode
251 if self.hamiltonian.poisson.get_description() == 'FDTD+TDDFT':
252 self.hamiltonian.poisson.set_density(self.density)
253 self.hamiltonian.poisson.print_messages(self.log)
254 self.log.flush()
256 # Update density and Hamiltonian
257 self.propagator.update_time_dependent_operators(self.time)
259 # XXX remove dipole moment handling and use observer instead
260 self._dm_args0 = (self.density.finegd.integrate(self.density.rhot_g),
261 self.calculate_dipole_moment())
263 # Call observers
264 self.action = 'init'
265 self.call_observers(self.niter)
267 self.tddft_initialized = True
269 def create_wave_functions(self, mode, *args, **kwargs):
270 mode = FDTDDFTMode(mode.nn, mode.interpolation, True)
271 GPAW.create_wave_functions(self, mode, *args, **kwargs)
273 def read(self, filename):
274 reader = GPAW.read(self, filename)
275 if 'tddft' in reader:
276 r = reader.tddft
277 self.time = r.time
278 self.niter = r.niter
279 self.kick_strength = r.kick_strength
281 def _write(self, writer, mode):
282 GPAW._write(self, writer, mode)
283 w = writer.child('tddft')
284 w.write(time=self.time,
285 niter=self.niter,
286 kick_strength=self.kick_strength)
288 def propagate(self, time_step: float, iterations: int,
289 dipole_moment_file: str = None, # deprecated
290 restart_file: str = None, # deprecated
291 dump_interval: int = None, # deprecated
292 ):
293 """Propagates wavefunctions.
295 Parameters
296 ----------
297 time_step
298 Time step in attoseconds (10^-18 s).
299 iterations
300 Number of iterations.
301 dipole_moment_file
302 Deprecated. Do not use this.
303 Name of the data file where to the time-dependent dipole
304 moment is saved.
305 restart_file
306 Deprecated. Do not use this.
307 Name of the restart file.
308 dump_interval
309 Deprecated. Do not use this.
310 After how many iterations restart data is dumped.
311 """
312 self.tddft_init()
314 if self.propagator.todict()['name'] in ['EFSICN', 'EFSICN_HGH']:
315 msg = ("You are using propagator for Ehrenfest dynamics. "
316 "Please use regular propagator.")
317 raise RuntimeError(msg)
319 def warn_deprecated(parameter, observer):
320 warnings.warn(
321 f"The {parameter} parameter is deprecated. "
322 f"Please use {observer} observer instead. "
323 "The old syntax will throw an error in the future.",
324 FutureWarning)
326 if dipole_moment_file is not None:
327 warn_deprecated('dipole_moment_file', 'DipoleMomentWriter')
329 if restart_file is not None:
330 warn_deprecated('restart_file', 'RestartFileWriter')
332 if dump_interval is not None:
333 warn_deprecated('dump_interval', 'RestartFileWriter')
334 else:
335 dump_interval = 100
337 if self.rank == 0:
338 self.log()
339 self.log('Starting time: %7.2f as'
340 % (self.time * autime_to_attosec))
341 self.log('Time step: %7.2f as' % time_step)
342 header = """\
343 Simulation Total log10 Iterations:
344 Time time Energy (eV) Norm Propagator"""
345 self.log()
346 self.log(header)
348 # Convert to atomic units
349 time_step = time_step * attosec_to_autime
351 if dipole_moment_file is not None:
352 self.initialize_dipole_moment_file(dipole_moment_file)
354 # Set these as class properties for use of observers
355 self.time_step = time_step
356 self.dump_interval = dump_interval # XXX remove, deprecated
357 self.restart_file = restart_file # XXX remove, deprecated
359 niterpropagator = 0
360 self.maxiter = self.niter + iterations
362 # FDTD requires extra care
363 if self.hamiltonian.poisson.get_description() == 'FDTD+TDDFT':
364 self.hamiltonian.poisson.set_time(self.time)
365 self.hamiltonian.poisson.set_time_step(self.time_step)
367 # The propagate calculation_mode causes classical part to evolve
368 # in time when self.hamiltonian.poisson.solve(...) is called
369 self.hamiltonian.poisson.set_calculation_mode('propagate')
371 # During each time step, self.hamiltonian.poisson.solve may be
372 # called several times (depending on the used propagator).
373 # Using the attached observer one ensures that actual propagation
374 # takes place only once. This is because the FDTDPoissonSolver
375 # changes the calculation_mode from propagate to
376 # something else when the propagation is finished.
377 self.attach(self.hamiltonian.poisson.set_calculation_mode, 1,
378 'propagate')
380 self.timer.start('Propagate')
381 while self.niter < self.maxiter:
382 norm = self.density.finegd.integrate(self.density.rhot_g)
384 # Write dipole moment at every iteration
385 if dipole_moment_file is not None:
386 if self._dm_args0 is not None:
387 self.update_dipole_moment_file(*self._dm_args0)
388 self._dm_args0 = None
389 else:
390 dm = self.calculate_dipole_moment()
391 self.update_dipole_moment_file(norm, dm)
393 # Propagate the Kohn-Shame wavefunctions a single timestep
394 niterpropagator = self.propagator.propagate(self.time, time_step)
395 self.time += time_step
396 self.niter += 1
398 # print output (energy etc.) every 10th iteration
399 if self.niter % 10 == 0:
400 self.get_td_energy()
402 T = time.localtime()
403 if self.rank == 0:
404 iter_text = 'iter: %3d %02d:%02d:%02d %11.2f' \
405 ' %13.6f %9.1f %10d'
406 self.log(iter_text %
407 (self.niter, T[3], T[4], T[5],
408 self.time * autime_to_attosec,
409 self.Etot * aufrequency_to_eV,
410 log(abs(norm) + 1e-16) / log(10),
411 niterpropagator))
413 self.log.flush()
415 # Call registered callback functions
416 self.action = 'propagate'
417 self.call_observers(self.niter)
419 # Write restart data
420 if restart_file is not None and self.niter % dump_interval == 0:
421 self.write(restart_file, 'all')
422 if self.rank == 0:
423 print('Wrote restart file.')
424 print(self.niter, ' iterations done. Current time is ',
425 self.time * autime_to_attosec, ' as.')
427 self.timer.stop('Propagate')
429 # Write final results and close dipole moment file
430 if dipole_moment_file is not None:
431 # TODO final iteration is propagated, but nothing is updated
432 # norm = self.density.finegd.integrate(self.density.rhot_g)
433 # self.finalize_dipole_moment_file(norm)
434 self.finalize_dipole_moment_file()
436 # Finalize FDTDPoissonSolver
437 if self.hamiltonian.poisson.get_description() == 'FDTD+TDDFT':
438 self.hamiltonian.poisson.finalize_propagation()
440 if restart_file is not None:
441 self.write(restart_file, 'all')
443 def initialize_dipole_moment_file(self, dipole_moment_file):
444 if self.rank == 0:
445 if self.dm_file is not None and not self.dm_file.closed:
446 raise RuntimeError('Dipole moment file is already open')
448 if self.time == 0.0:
449 mode = 'w'
450 else:
451 # We probably continue from restart
452 mode = 'a'
454 self.dm_file = open(dipole_moment_file, mode)
456 # If the dipole moment file is empty, add a header
457 if self.dm_file.tell() == 0:
458 header = '# Kick = [%22.12le, %22.12le, %22.12le]\n' \
459 % (self.kick_strength[0], self.kick_strength[1],
460 self.kick_strength[2])
461 header += '# %15s %15s %22s %22s %22s\n' \
462 % ('time', 'norm', 'dmx', 'dmy', 'dmz')
463 self.dm_file.write(header)
464 self.dm_file.flush()
466 def calculate_dipole_moment(self):
467 dm = self.density.finegd.calculate_dipole_moment(self.density.rhot_g)
468 if self.hamiltonian.poisson.get_description() == 'FDTD+TDDFT':
469 dm += self.hamiltonian.poisson.get_classical_dipole_moment()
470 return dm
472 def update_dipole_moment_file(self, norm, dm):
473 if self.rank == 0:
474 line = '%20.8lf %20.8le %22.12le %22.12le %22.12le\n' \
475 % (self.time, norm, dm[0], dm[1], dm[2])
476 self.dm_file.write(line)
477 self.dm_file.flush()
479 def finalize_dipole_moment_file(self, norm=None):
480 if norm is not None:
481 dm = self.calculate_dipole_moment()
482 self.update_dipole_moment_file(norm, dm)
484 if self.rank == 0:
485 self.dm_file.close()
486 self.dm_file = None
488 def update_eigenvalues(self):
490 kpt_u = self.wfs.kpt_u
491 if self.hpsit is None:
492 self.hpsit = self.wfs.gd.zeros(len(kpt_u[0].psit_nG),
493 dtype=complex)
494 if self.eps_tmp is None:
495 self.eps_tmp = np.zeros(len(kpt_u[0].eps_n),
496 dtype=complex)
498 # self.Eband = sum_i <psi_i|H|psi_j>
499 for kpt in kpt_u:
500 self.td_hamiltonian.apply(kpt, kpt.psit_nG, self.hpsit,
501 calculate_P_ani=False)
502 self.mblas.multi_zdotc(self.eps_tmp, kpt.psit_nG,
503 self.hpsit, len(kpt_u[0].psit_nG))
504 self.eps_tmp *= self.wfs.gd.dv
505 kpt.eps_n[:] = self.eps_tmp.real
507 H = self.td_hamiltonian.hamiltonian
509 # Nonlocal
510 self.Enlkin = H.xc.get_kinetic_energy_correction()
512 # PAW
513 e_band = self.wfs.calculate_band_energy()
514 self.Ekin = H.e_kinetic0 + e_band + self.Enlkin
515 self.e_coulomb = H.e_coulomb
516 self.Eext = H.e_external
517 self.Ebar = H.e_zero
518 self.Exc = H.e_xc
519 self.Etot = self.Ekin + self.e_coulomb + self.Ebar + self.Exc
521 def get_td_energy(self):
522 """Calculate the time-dependent total energy"""
523 self.tddft_init()
525 if not self.calculate_energy:
526 self.Etot = 0.0
527 return 0.0
529 self.wfs.overlap.update(self.wfs)
530 self.td_density.update()
531 self.td_hamiltonian.update(self.td_density.get_density(),
532 self.time)
533 self.update_eigenvalues()
535 return self.Etot
537 def set_absorbing_boundary(self, absorbing_boundary):
538 self.td_hamiltonian.set_absorbing_boundary(absorbing_boundary)
540 # exp(ip.r) psi
541 def absorption_kick(self, kick_strength):
542 """Delta absorption kick for photoabsorption spectrum.
544 Parameters
545 ----------
546 kick_strength
547 Strength of the kick in atomic units
548 """
549 self.tddft_init()
551 if self.rank == 0:
552 self.log('Delta kick =', kick_strength)
554 self.kick_strength = np.array(kick_strength)
556 abs_kick_hamiltonian = AbsorptionKickHamiltonian(
557 self.wfs, self.spos_ac,
558 np.array(kick_strength, float))
559 abs_kick = AbsorptionKick(self.wfs, abs_kick_hamiltonian,
560 self.wfs.overlap, self.solver,
561 self.preconditioner, self.wfs.gd, self.timer)
562 abs_kick.kick()
564 # Kick the classical part, if it is present
565 if self.hamiltonian.poisson.get_description() == 'FDTD+TDDFT':
566 self.hamiltonian.poisson.set_kick(kick=self.kick_strength)
568 # Update density and Hamiltonian
569 self.propagator.update_time_dependent_operators(self.time)
571 # Call observers after kick
572 self.action = 'kick'
573 self.call_observers(self.niter)