Time-propagation TDDFT

Optical photoabsorption spectrum as well as nonlinear effects can be studied using time propagation TDDFT. This approach scales better than linear response, but the prefactor is so large that for small and moderate systems linear response is significantly faster.

Ground state

To obtain the ground state for TDDFT, one has to just do a standard ground state with a larger simulation box. A proper distance from any atom to edge of the simulation box is problem dependent, but a minimum reasonable value is around 6 Ångströms and recommended between 8-10 Ång. In TDDFT, one can use larger grid spacing than for geometry optimization. For example, if you use h=0.25 for geometry optimization, try h=0.3 for TDDFT. This saves a lot of time.

A good way to start is to use too small box (vacuum=6.0), too large grid spacing (h=0.35), and too large time step (dt=16.0). Then repeat the simulation with better parameter values and compare. Probably lowest peaks are already pretty good, and far beyond the ionization limit, in the continuum, the spectrum is not going to converge anyway. The first run takes only fraction of the time of the second run.

For a parallel-over-states TDDFT calculation, you must choose the number of states so, that these can be distributed equally to processors. For example, if you have 79 occupied states and you want to use 8 processes in parallelization over states, add one unoccupied state to get 80 states in total.

Ground state example:

# Standard magic
from ase import Atoms
from gpaw import GPAW

# Beryllium atom
atoms = Atoms(symbols='Be',
              positions=[(0, 0, 0)],

# Add 6.0 ang vacuum around the atom

# Create GPAW calculator
calc = GPAW(nbands=1, h=0.3)
# Attach calculator to atoms
atoms.calc = calc

# Calculate the ground state
energy = atoms.get_potential_energy()

# Save the ground state
calc.write('be_gs.gpw', 'all')

Time-propagation TDDFT is also available in LCAO mode.

Optical photoabsorption spectrum

Optical photoabsorption spectrum can be obtained by applying a weak delta pulse of dipole electric field, and then letting the system evolve freely while recording the dipole moment. A time-step around 4.0-8.0 attoseconds is reasonable. The total simulation time should be few tens of femtoseconds depending on the desired resolution.


from gpaw.tddft import *

time_step = 8.0                  # 1 attoseconds = 0.041341 autime
iterations = 2500                # 2500 x 8 as => 20 fs
kick_strength = [0.0,0.0,1e-3]   # Kick to z-direction

# Read ground state
td_calc = TDDFT('be_gs.gpw')

# Save the time-dependent dipole moment to 'be_dm.dat'
DipoleMomentWriter(td_calc, 'be_dm.dat')

# Use 'be_td.gpw' as restart file
RestartFileWriter(td_calc, 'be_td.gpw')

# Kick with a delta pulse to z-direction

# Propagate
td_calc.propagate(time_step, iterations)

# Save end result to 'be_td.gpw'
td_calc.write('be_td.gpw', mode='all')

# Calculate photoabsorption spectrum and write it to 'be_spectrum_z.dat'
photoabsorption_spectrum('be_dm.dat', 'be_spectrum_z.dat')

When propagating after an absorption kick has been applied, it is a good idea to periodically write the time-evolution state to a restart file. This ensures that you can resume adding data to the dipole moment file if you experience artificial oscillations in the spectrum because the total simulation time was too short.


from gpaw.tddft import *

time_step = 8.0                  # 1 attoseconds = 0.041341 autime
iterations = 2500                # 2500 x 8 as => 20 fs

# Read restart file with result of previous propagation
td_calc = TDDFT('be_td.gpw')

# Append the time-dependent dipole moment to the already existing 'be_dm.dat'
DipoleMomentWriter(td_calc, 'be_dm.dat')

# Use 'be_td2.gpw' as restart file
RestartFileWriter(td_calc, 'be_td2.gpw')

# Propagate more
td_calc.propagate(time_step, iterations)

# Save end result to 'be_td2.gpw'
td_calc.write('be_td2.gpw', mode='all')

# Recalculate photoabsorption spectrum and write it to 'be_spectrum_z2.dat'
photoabsorption_spectrum('be_dm.dat', 'be_spectrum_z2.dat')

Typically in experiments, the spherically averaged spectrum is measured. To obtain this, one must repeat the time-propagation to each Cartesian direction and average over the Fourier transformed dipole moments.

Fourier transformed density and electric near field

To analyze the resonances seen in the photoabsorption spectrum, see Induced density oscillations and electric near field from TDDFT.

Time propagation

Since the total CPU time also depends on the number of iterations performed by the linear solvers in each time-step, smaller time-steps around 2.0-4.0 attoseconds might prove to be faster with the ECN and SICN propagators because they have an embedded Euler step in each predictor step:

\[\tilde{\psi}_n(t+\Delta t) \approx (1 - i \hat{S}^{\;-1}_\mathrm{approx.}(t) \tilde{H}(t) \Delta t)\tilde{\psi}_n(t)\]

, where \(\hat{S}^{\;-1}_\mathrm{approx.}\) is an inexpensive operation which approximates the inverse of the overlap operator \(\hat{S}\). See the Developers Guide for details.

Therefore, as a rule-of-thumb, choose a time-step small enough to minimize the number of iterations performed by the linear solvers in each time-step, but large enough to minimize the number of time-steps required to arrive at the desired total simulation time.

TDDFT reference manual

class gpaw.tddft.TDDFT(filename: str, *, td_potential: object = None, calculate_energy: bool = True, propagator: Optional[dict] = None, solver: Optional[dict] = None, tolerance: Optional[float] = None, parallel: Optional[dict] = None, communicator: object = None, txt: str = '-')[source]

Time-dependent density functional theory calculation based on GPAW.

This class is the core class of the time-dependent density functional theory implementation and is the only class which a user has to use.

  • filename – File containing ground state or time-dependent state to propagate.

  • td_potential – Function class for the time-dependent potential. Must have a method strength(time) which returns the strength of the linear potential to each direction as a vector of three floats.

  • calculate_energy – Whether to calculate energy during propagation.

  • propagator – Time propagator for the Kohn-Sham wavefunctions.

  • solver – The iterative linear equations solver for propagator.

  • tolerance – Deprecated. Do not use this, but use solver dictionary instead. Tolerance for the linear solver.

  • parallel – Parallelization options

  • communicator – MPI communicator

  • txt – Text output


Delta absorption kick for photoabsorption spectrum.


kick_strength – Strength of the kick in atomic units


Calculate the time-dependent total energy

propagate(time_step: float, iterations: int, dipole_moment_file: Optional[str] = None, restart_file: Optional[str] = None, dump_interval: Optional[int] = None)[source]

Propagates wavefunctions.

  • time_step – Time step in attoseconds (10^-18 s).

  • iterations – Number of iterations.

  • dipole_moment_file – Deprecated. Do not use this. Name of the data file where to the time-dependent dipole moment is saved.

  • restart_file – Deprecated. Do not use this. Name of the restart file.

  • dump_interval – Deprecated. Do not use this. After how many iterations restart data is dumped.

class gpaw.tddft.DipoleMomentWriter(paw, filename, center=False, density='comp', force_new_file=False, interval=1)[source]
class gpaw.tddft.RestartFileWriter(paw, restart_filename, interval=100)[source]

Observer for writing restart files periodically.

At the given interval, the calculator restart file is written and write_restart() of every attached observer is called.

The observer attaches to the TDDFT calculator during creation.

  • paw – TDDFT calculator

  • restart_filename – File for writing the calculator object

  • interval – Update interval. The restart files are written every that many propagation steps.

gpaw.tddft.photoabsorption_spectrum(dipole_moment_file: str, spectrum_file: str, folding: str = 'Gauss', width: float = 0.2123, e_min: float = 0.0, e_max: float = 30.0, delta_e: float = 0.05)[source]

Calculates photoabsorption spectrum from the time-dependent dipole moment.

  • dipole_moment_file – Name of the time-dependent dipole moment file from which the spectrum is calculated

  • spectrum_file – Name of the spectrum file

  • folding – Gaussian ('Gauss') or Lorentzian ('Lorentz') folding

  • width – Width of the Gaussian (sigma) or Lorentzian (Gamma) Gaussian = 1/(sigma sqrt(2pi)) exp(-(1/2)(omega/sigma)^2) Lorentzian = (1/pi) (1/2) Gamma / [omega^2 + ((1/2) Gamma)^2]

  • e_min – Minimum energy shown in the spectrum (eV)

  • e_max – Maximum energy shown in the spectrum (eV)

  • delta_e – Energy resolution (eV)