Linear response TDDFT

Ground state

The linear response TDDFT calculation needs a converged ground state calculation with a set of unoccupied states. The standard 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(nbands=10, h=0.3)
# Attach calculator to atoms
atoms.set_calculator(calc)

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

# converge also the empty states (the density is converged already)
calc.set(convergence={'bands': 8},
         fixdensity=True,
         eigensolver='cg')
atoms.get_potential_energy()

# 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', 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', 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('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

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('lr.dat.gz')
lr.diagonalize(energy_range=2.)

This can be automated by using the check_convergence function:

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

lr = LrTDDFT('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('lr.dat.gz')
lr.diagonalize()

# analyse transition 1
lr.analyse(1)

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

Relaxation in the excited state

This example shows how to relax in the B excited state of the sodium dimer:

from ase import Atom, io, optimize
from gpaw import GPAW, FermiDirac
from gpaw.cluster import Cluster
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

s = Cluster([Atom('Na'), Atom('Na', [0, 0, R])])
s.minimal_box(box, h=h)

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

ex = ExcitedState(lr, iex, d=d)
s.set_calculator(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)

Quick reference

Parameters for LrTDDFT:

keyword type default value description
calculator GPAW   Calculator object of ground state calculation
filename string   read the state of LrTDDFT calculation (i.e. omega matrix, excitations) from filename
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])
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
eps float 0.001 Minimal occupation difference for a transition