# 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)],
pbc=False)

# Add 6.0 ang vacuum around the atom
atoms.center(vacuum=6.0)

# Create GPAW calculator
calc = GPAW(nbands=1, h=0.3)
# Attach calculator to atoms
atoms.set_calculator(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.

Example:

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

td_calc = TDDFT('be_gs.gpw')

# Kick with a delta pulse to z-direction
td_calc.absorption_kick(kick_strength=kick_strength)

# Propagate, save the time-dependent dipole moment to 'be_dm.dat',
# and use 'be_td.gpw' as restart file
td_calc.propagate(time_step, iterations, 'be_dm.dat', 'be_td.gpw')

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


Note

Make sure to number of iterations is divisible by the dump interval such that the last iteration will be stored in the restart file. Otherwise append td_calc.write(‘be_td.gpw’, mode=’all’) to the script.

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.

Example:

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')

# Propagate more, appending the time-dependent dipole moment to the
# already existing 'be_dm.dat' and use 'be_td2.gpw' as restart file
td_calc.propagate(time_step, iterations, 'be_dm.dat', 'be_td2.gpw')

# 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¶

If one merely wishes to record the time-evolution of the dipole moment and analyze the resulting spectrum, TDDFT offers little advantage over the much faster LrTDDFT. Further, one must bear in mind that only excitations induced by the absorption kick will show up in the spectrum.

However, propagating a slightly perturbed ground state density may offer much more structural information, starting with the ability to distinguish which spectral peaks correspond to which principal directions in a lattice.

Since the dipole moment is generated by displacements in the charge density, most strong peaks in the optical photoabsorption spectrum signify nearly harmonic oscillations herein. Therefore, taking Fourier transforms of the time-evolution of the density at the resonant frequencies is a great way of analyzing the spatial extent of the oscillating modes.

The discrete moving-average Fourier transform of the pseudo-electron density $$\tilde{n}(\mathbf{r},t)$$ is defined:

$F_N(\mathbf{r},\omega) = \frac{1}{\sqrt{\pi}} \sum_{j=0}^N \big( \tilde{n}(\mathbf{r},t_j)-\overline{n}_N(\mathbf{r})\big) \mathrm{e}^{-\textstyle\frac{1}{2}t_j^2\sigma^2} \mathrm{e}^{\displaystyle\mathrm{i}\omega t_j} \Delta t_j$

, where we allow for variable time-step $$\Delta t_j$$ along the $$N$$ propagation steps in the time-series $$j=0,1,\ldots,N$$. With a total propagation time of $$t_N$$, the Fourier transforms are taken relative to the time-average $$\overline{n}_N(\mathbf{r})$$ of the pseudo density:

$\overline{n}_N(\mathbf{r}) = \frac{1}{t_{N+1}} \sum_{j=0}^N \tilde{n}(\mathbf{r},t_j) \Delta t_j \qquad, t_N = \sum_{j=0}^{N-1}\Delta t_j$

Regrettably, having arrived at time $$t_N$$ will not enable us to perform the above summations because recording $$N\sim 10^4$$ sets of grid data is completely intractable. Instead, an iterative cumulation scheme is implemented, which only requires data from one time-step at a time.

XXX more on this later

The class gpaw.tddft.fourier.DensityFourierTransform is used to calculate and maintain Fourier transforms of the pseudo electron density. It functions by attaching itself to a TDDFT instance, which in turn notifies it after each time-step and allows it to update the density Fourier transforms.

Important

An incontestable restriction of the iterative approach is the requirement that the frequencies must be given upon initialization (i.e. time zero). To avoid wasted effort, getting the peak frequencies right is essential.

It is recommended to use either LrTDDFT or a somewhat cruder time-propagation to estimate which frequencies could be of interest. In the latter case, applying a weak kick [1e-3, 1e-3, 1e-3] will probably be sufficient to excite and detect all the relevant modes in a short time-span. For quick estimates, using the ECN propagator and the CSCG eigensolver with a tolerance around 1e-4 works reasonably well for timesteps of 5-10 as.

Tip

Using a finite width $$\sigma$$ around 0.1 eV will make any ballpark figure a much safer bet. Be aware that peaks found using LrTDDFT may shift slightly.

Example:

from gpaw.tddft import TDDFT
from gpaw.tddft.fourier import DensityFourierTransform

time_step = 4.0                  # 1 attoseconds = 0.041341 autime
iterations = 5000                # 5000 x 4 as => 20 fs
kick_strength = [0.0,5e-3,0.0]   # Kick to y-direction
frequencies = [4.26,6.27,13.0, \
16.9,18.1,19.9]   # Pre-determined peak frequencies in eV
sigma = 0.05                     # Width of Gaussian envelope in eV

td_calc = TDDFT('bda_gs.gpw')

# Kick with a delta pulse to y-direction
td_calc.absorption_kick(kick_strength=kick_strength)

# Create and attach Fourier transform observer
obs = DensityFourierTransform(timestep, frequencies, sigma)
obs.initialize(td_calc)

# Propagate, save the time-dependent dipole moment to 'bda_dm.dat',
# (just for comparison) and use 'bda_td.gpw' as restart file
td_calc.propagate(time_step, iterations, 'bda_dm.dat', 'bda_td.gpw')

# Save result of the Fourier transformations to a .ftd file
obs.write('bda_fourier.ftd')


You can now resume adding data to both the dipole moment file and the density fourier transform if the spectrum is not sufficiently evolved because the total simulation time was too short.

Example:

from gpaw.tddft import TDDFT
from gpaw.tddft.fourier import DensityFourierTransform

time_step = 4.0                  # 1 attoseconds = 0.041341 autime
iterations = 5000                # 5000 x 4 as => 20 fs
frequencies = [4.26,6.27,13.0, \
16.9,18.1,19.9]   # Pre-determined peak frequencies in eV
sigma = 0.05                     # Width of Gaussian envelope in

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

# Create and attach Fourier transform observer
obs = DensityFourierTransform(timestep, frequencies, sigma)
obs.initialize(td_calc)

# Read previous result of the corresponding Fourier transformations

# Propagate more, appending the time-dependent dipole moment to the
# already existing 'bda_dm.dat' and use 'bda_td2.gpw' as restart file
td_calc.propagate(time_step, iterations, 'bda_dm.dat', 'bda_td2.gpw')

# Save result of the improved Fourier transformations to an .ftd file
obs.write('bda_fourier2.ftd')


## 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¶

The TDDFT class and keywords:

Keyword

Type

Default

Description

ground_state_file

string

Name of the ground state file

td_potential

TDPotential

None

Time-dependent external potential

propagator

string

'SICN'

Time-propagator ('ECN'/'SICN'/'SITE'/'SIKE')

solver

string

'CSCG'

Linear equation solver ('CSCG'/'BiCGStab')

tolerance

float

1e-8

Tolerance for linear solver

Keywords for absorption_kick():

Keyword

Type

Default

Description

kick_strength

float[3]

[0,0,1e-3]

Kick strength

Keywords for propagate():

Keyword

Type

Default

Description

time_step

float

Time step in attoseconds (1 autime = 24.188 as)

iterations

integer

Iterations

dipole_moment_file

string

None

Name of the dipole moment file

restart_file

string

None

Name of the restart file

dump_interal

integer

500

How often restart file is written

Keywords for gpaw.tddft.photoabsorption_spectrum():

Keyword

Type

Default

Description

dipole_moment_file

string

Name of the dipole moment file

spectrum_file

string

Name of the spectrum file

folding

string

Gauss

Gaussian folding (or Lorentzian in future)

width

float

0.2123

Width of the Gaussian/Lorentzian (in eV)

e_min

float

0.0

Lowest energy shown in spectrum (in eV)

e_max

float

30.0

Highest energy shown in spectrum (in eV)

delta_e

float

0.05

Resolution of energy in spectrum (in eV)

class gpaw.tddft.TDDFT(filename, td_potential=None, propagator='SICN', calculate_energy=True, propagator_kwargs=None, solver='CSCG', tolerance=1e-08, **kwargs)[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.

Create TDDFT-object.

Parameters:

filename: string

File containing ground state or time-dependent state to propagate

td_potential: class, optional

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.

propagator: {‘SICN’,’ETRSCN’,’ECN’,’SITE’,’SIKE4’,’SIKE5’,’SIKE6’}

Name of the time propagator for the Kohn-Sham wavefunctions

solver: {‘CSCG’,’BiCGStab’}

Name of the iterative linear equations solver for time propagation

tolerance: float

Tolerance for the linear solver

The following parameters can be used: $$txt$$, $$parallel$$, $$communicator$$ $$mixer$$ and $$dtype$$. The internal parameters $$mixer$$ and $$dtype$$ are strictly used to specify a dummy mixer and complex type respectively.

absorption_kick(kick_strength)[source]

Delta absorption kick for photoabsorption spectrum.

Parameters:

kick_strength: [float, float, float]

Strength of the kick, e.g., [0.0, 0.0, 1e-3]

get_td_energy()[source]

Calculate the time-dependent total energy

initialize(reading=False)[source]

Inexpensive initialization.

propagate(time_step, iterations, dipole_moment_file=None, restart_file=None, dump_interval=100)[source]

Propagates wavefunctions.

Parameters:

time_step: float

Time step in attoseconds (10^-18 s), e.g., 4.0 or 8.0

iterations: integer

Iterations, e.g., 20 000 as / 4.0 as = 5000

dipole_moment_file: string, optional

Name of the data file where to the time-dependent dipole moment is saved

restart_file: string, optional

Name of the restart file

dump_interval: integer

After how many iterations restart data is dumped

read(filename)[source]

Read atoms, parameters and calculated properties from output file.

Read result from self.label file. Raise ReadError if the file is not there. If the file is corrupted or contains an error message from the calculation, a ReadError should also be raised. In case of succes, these attributes must set:

atoms: Atoms object

The state of the atoms from last calculation.

parameters: Parameters object

The parameter dictionary.

results: dict

Calculated properties like energy and forces.

The FileIOCalculator.read() method will typically read atoms and parameters and get the results dict by calling the read_results() method.

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

Calculates photoabsorption spectrum from the time-dependent dipole moment.

Parameters:

dipole_moment_file: string

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

spectrum_file: string

Name of the spectrum file

folding: ‘Gauss’ or ‘Lorentz’

Whether to use Gaussian or Lorentzian folding

width: float

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: float

Minimum energy shown in the spectrum (eV)

e_max: float

Maximum energy shown in the spectrum (eV)

delta_e: float

Energy resolution (eV)

class gpaw.tddft.fourier.DensityFourierTransform(timestep, frequencies, width=None, interval=1)[source]

Parameters

timestep: float

Time step in attoseconds (10^-18 s), e.g., 4.0 or 8.0

frequencies: NumPy array or list of floats

Frequencies in eV for Fourier transforms

width: float or None

Width of Gaussian envelope in eV, otherwise no envelope

interval: int

Number of timesteps between calls (used when attaching)