Linear response TDDFT

Ground state

The linear response TDDFT calculation needs a converged ground state calculation with a set of unoccupied states. Note, that the eigensolver ‘rmm-diis’ should not be used for the calculation of unoccupied states, better use ‘dav’ or ‘cg’:

from ase import Atoms
from gpaw import GPAW

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

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

# Create GPAW calculator
calc = GPAW(mode='fd', nbands=10, h=0.3)
# Attach calculator to atoms
atoms.calc = calc

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

# converge also the empty states (the density is converged already)
calc = calc.fixed_density(
    convergence={'bands': 8})

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

Calculating the Omega Matrix

The next step is to calculate the Omega Matrix from the ground state orbitals:

from gpaw import GPAW
from gpaw.lrtddft import LrTDDFT

c = GPAW('Be_gs_8bands.gpw')

istart = 0  # band index of the first occ. band to consider
jend = 8  # band index of the last unocc. band to consider
lr = LrTDDFT(c, xc='LDA', restrict={'istart': istart, 'jend': jend},
             nspins=2)  # force the calculation of triplet excitations also
lr.write('lr.dat.gz')

alternatively one can also restrict the number of transitions by their energy:

from gpaw import GPAW
from gpaw.lrtddft import LrTDDFT

c = GPAW('Be_gs_8bands.gpw')

dE = 10  # maximal Kohn-Sham transition energy to consider in eV
lr = LrTDDFT(c, xc='LDA', restrict={'energy_range': dE})
lr.write('lr_dE.dat.gz')

Note, that parallelization over spin does not work here. As a workaround, domain decomposition only (parallel={'domain': world.size}, see Domain decomposition) has to be used for spin polarised calculations in parallel.

Extracting the spectrum

The dipole spectrum can be evaluated from the Omega matrix and written to a file:

from gpaw.lrtddft import LrTDDFT
from gpaw.lrtddft import photoabsorption_spectrum

lr = LrTDDFT.read('lr.dat.gz')
lr.diagonalize()

# write the spectrum to the data file
photoabsorption_spectrum(lr,
                         'spectrum_w.05eV.dat',  # data file name
                         width=0.05)             # width in eV

The spectrum may be also extracted and plotted in energy or wavelength directly:

import matplotlib.pyplot as plt
from gpaw.lrtddft import LrTDDFT
from gpaw.lrtddft.spectrum import get_folded_spectrum

lr = LrTDDFT.read('lr.dat.gz')

plt.subplot(121)
# spectrum in energy
x, y = get_folded_spectrum(lr, width=0.05)
plt.plot(x, y[:, 0])
plt.xlabel('energy [eV]')
plt.ylabel('folded oscillator strength')

plt.subplot(122)
# spectrum in wavelengths
x, y = get_folded_spectrum(lr, energyunit='nm', width=10)
plt.plot(x, y[:, 0])
plt.xlabel('wavelength [nm]')

plt.show()

Testing convergence

You can test the convergence of the Kohn-Sham transition basis size by restricting the basis in the diagonalisation step, e.g.:

from gpaw.lrtddft import LrTDDFT

lr = LrTDDFT.read('lr.dat.gz')
lr.diagonalize(restrict={'energy_range':6.})

This can be automated by using the check_convergence function:

from gpaw.lrtddft import LrTDDFT
from gpaw.lrtddft.convergence import check_convergence

lr = LrTDDFT.read('lr.dat.gz')
check_convergence(lr,
                  'linear_response',
                  'my plot title',
                   dE=.2,
                   emax=6.)

which will create a directory ‘linear_response’. In this directory there will be a file ‘conv.gpl’ for gnuplot that compares the spectra varying the basis size.

Analysing the transitions

The single transitions (or a list of transitions) can be analysed as follows (output printed):

from gpaw.lrtddft import LrTDDFT

lr = LrTDDFT.read('lr.dat.gz')
lr.diagonalize()

# analyse transition 1
lr.analyse(1)

# analyse transition 0-10
lr.analyse(range(11))

Relaxation in the excited state

Despite that we do not have analytical gradients in linear response TDDFT, we may use finite differences for relaxation in the excited state. This example shows how to relax the sodium dimer in the B excited state:

from ase import Atoms, io, optimize
from gpaw import GPAW, FermiDirac
from gpaw.utilities.adjust_cell import adjust_cell
from gpaw.lrtddft import LrTDDFT
from gpaw.lrtddft.excited_state import ExcitedState

box = 5.      # box dimension
h = 0.25      # grid spacing
width = 0.01  # Fermi width
nbands = 6    # bands in GS calculation
nconv = 4     # bands in GS calculation to converge
R = 2.99      # starting distance
iex = 1       # excited state index
d = 0.01      # step for numerical force evaluation
exc = 'LDA'   # xc for the linear response TDDFT kernel
parallel = 4  # number of independent parrallel workers for force evaluation

s = Atoms('Na2', [[0, 0, 0], [0, 0, R]])
adjust_cell(s, box, h=h)

c = GPAW(mode='fd', h=h, nbands=nbands, eigensolver='cg',
         occupations=FermiDirac(width=width),
         setups={'Na': '1'},
         convergence={'bands': nconv})
c.calculate(s)
lr = LrTDDFT(c, xc=exc, restrict={'eps': 0.1, 'jend': nconv - 1})

ex = ExcitedState(lr, iex, d=d, parallel=parallel)
s.calc = ex

ftraj = 'relax_ex' + str(iex)
ftraj += '_box' + str(box) + '_h' + str(h)
ftraj += '_d' + str(d) + '.traj'
traj = io.Trajectory(ftraj, 'w', s)
dyn = optimize.FIRE(s)
dyn.attach(traj.write)
dyn.run(fmax=0.05)

The example runs on a single core. If started on 8 cores, it will split world into 4 independent parts (2 cores each) that can calculate 4 displacements in parallel at the same time.

Quick reference

Parameters for LrTDDFT:

keyword	type	default value	description
`calculator`	`GPAW`		Calculator object of ground state calculation
`nspins`	`int`	1	number of excited state spins, i.e. singlet-triplet transitions are calculated with `nspins=2`. Effective only if ground state is spin-compensated
`xc`	`string`	xc of calculator	Exchange-correlation for LrTDDFT, can differ from ground state value. Use “RPA” for the random phase approximation.
`restrict`	`dict`	{}	Restrictions collected as dict, see below

The keyword restrict can also be used to restrict the already pre-calculated linear response TDDFT object as:

from gpaw.lrtddft import LrTDDFT


lr = LrTDDFT.read('lr.dat.gz')
lr.diagonalize()
print('Full matrix gives', len(lr), 'excitations')

# now diagonalize using a selectron of orbitals
restrict = {'energy_range': 5}  # eV
lr.diagonalize(restrict=restrict)
print('Restriction by', restrict, 'gives', len(lr), 'excitation')

restrict = {'to': [3, 4]}
lr.diagonalize(restrict=restrict)
print('Restriction by', restrict, 'gives', len(lr), 'excitations')

The parameters of restrict are:

keyword	type	default value	description
`eps	`float`	0.001	Minimal occupation difference for a transition
`istart`	`int`	0	first occupied state to consider
`jend`	`int`	number of bands	last unoccupied state to consider
`energy_range`	`float`	None	Energy range to consider in the involved Kohn-Sham orbitals (replaces `[istart,jend]`)
`from`	`List[int]`	None	List of indices allowed for occupied Kohn-Sham orbitals
`to`	`List[int]`	None	List of indices allowed for unoccupied Kohn-Sham orbitals