Circular dichroism with real-time TDDFT
In this tutorial, we calculate the rotatory strength spectrum of
(R)-methyloxirane molecule
.
The equations underlying the implementation are described in Ref. [1].
LCAO mode
In this example, we use real-time TDDFT LCAO mode. We recall that the LCAO calculations are sensitive to the used basis sets, and we also demonstrate the construction of augmented basis sets.
We augment the default dzp basis sets with numerical Gaussian-type orbitals (NGTOs):
from gpaw.lcao.generate_ngto_augmented import \
create_GTO_dictionary as GTO, generate_nao_ngto_basis
gtos_atom = {'H': [GTO('s', 0.02974),
GTO('p', 0.14100)],
'C': [GTO('s', 0.04690),
GTO('p', 0.04041),
GTO('d', 0.15100)],
'O': [GTO('s', 0.07896),
GTO('p', 0.06856),
GTO('d', 0.33200)],
}
for atom in ['H', 'C', 'O']:
generate_nao_ngto_basis(atom, xc='PBE', name='aug',
nao='dzp', gtos=gtos_atom[atom])
Here, the Gaussian parameters correspond to diffuse functions in aug-cc-pvdz basis sets as obtained from Basis Set Exchange.
Let’s do the ground-state calculation (note the basis
keyword):
from ase.io import read
from gpaw import setup_paths
from gpaw import GPAW
# Insert the path to the created basis set
setup_paths.insert(0, '.')
atoms = read('r-methyloxirane.xyz')
atoms.center(vacuum=6.0)
calc = GPAW(mode='lcao',
basis='aug.dzp',
h=0.3,
xc='PBE',
nbands=16,
poissonsolver={
'name': 'MomentCorrectionPoissonSolver',
'poissonsolver': 'fast',
'moment_corrections': 1 + 3 + 5},
convergence={'density': 1e-12},
txt='gs.out',
symmetry={'point_group': False})
atoms.calc = calc
atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
Then, we carry the time propagation as usual in
real-time TDDFT LCAO mode,
but we attach MagneticMomentWriter
to record the time-dependent magnetic moment.
In this script, we wrap the time propagation code
inside main()
function to make the same script reusable
with different kick directions:
from gpaw import setup_paths
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.lcaotddft.magneticmomentwriter import MagneticMomentWriter
# Insert the path to the created basis set
setup_paths.insert(0, '.')
def main(kick):
kick_strength = [0., 0., 0.]
kick_strength['xyz'.index(kick)] = 1e-5
td_calc = LCAOTDDFT('gs.gpw', txt=f'td-{kick}.out')
DipoleMomentWriter(td_calc, f'dm-{kick}.dat')
MagneticMomentWriter(td_calc, f'mm-{kick}.dat')
td_calc.absorption_kick(kick_strength)
td_calc.propagate(10, 2000)
td_calc.write(f'td-{kick}.gpw', mode='all')
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('--kick', default='x')
kwargs = vars(parser.parse_args())
main(**kwargs)
After repeating the calculation for kicks in x, y, and z directions, we calculate the rotatory strength spectrum from the magnetic moments:
from gpaw.tddft.spectrum import rotatory_strength_spectrum
rotatory_strength_spectrum(['mm-x.dat', 'mm-y.dat', 'mm-z.dat'],
'rot_spec.dat',
folding='Gauss', width=0.2,
e_min=0.0, e_max=10.0, delta_e=0.01)
Comparing the resulting spectrum to one calculated with plain dzp basis sets shows the importance of augmented basis sets:
FD mode
In this example, we use real-time TDDFT FD mode.
Let’s do the ground-state calculation (note that FD mode typically requires larger vacuum and smaller grid spacing than LCAO mode; convergence with respect to these parameters need to be checked for both modes separately):
from ase.io import read
from gpaw import GPAW
atoms = read('r-methyloxirane.xyz')
atoms.center(vacuum=8.0)
calc = GPAW(mode='fd',
h=0.2,
xc='PBE',
nbands=16,
poissonsolver={
'name': 'MomentCorrectionPoissonSolver',
'poissonsolver': 'fast',
'moment_corrections': 1 + 3 + 5},
convergence={'density': 1e-12},
txt='gs.out',
symmetry={'point_group': False})
atoms.calc = calc
atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
Then, similarly to LCAO mode, we carry the time propagation as usual but attach
MagneticMomentWriter
to record the time-dependent magnetic moment
(note the tolerance
parameter for the iterative solver;
smaller values might be required when signal is weak):
from gpaw.tddft import TDDFT, DipoleMomentWriter, MagneticMomentWriter
def main(kick):
kick_strength = [0., 0., 0.]
kick_strength['xyz'.index(kick)] = 1e-5
td_calc = TDDFT('gs.gpw',
solver=dict(name='CSCG', tolerance=1e-8),
txt=f'td-{kick}.out')
DipoleMomentWriter(td_calc, f'dm-{kick}.dat')
MagneticMomentWriter(td_calc, f'mm-{kick}.dat')
td_calc.absorption_kick(kick_strength)
td_calc.propagate(10, 2000)
td_calc.write(f'td-{kick}.gpw', mode='all')
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('--kick', default='x')
kwargs = vars(parser.parse_args())
main(**kwargs)
After repeating the calculation for kicks in x, y, and z directions, we calculate the rotatory strength spectrum from the magnetic moments:
from gpaw.tddft.spectrum import rotatory_strength_spectrum
rotatory_strength_spectrum(['mm-x.dat', 'mm-y.dat', 'mm-z.dat'],
'rot_spec.dat',
folding='Gauss', width=0.2,
e_min=0.0, e_max=10.0, delta_e=0.01)
The resulting spectrum:
The spectrum compares well with the one obtained with LCAO mode and augmented basis sets.
Origin dependence
The circular dichroism spectra obtained with the present implementation
depend on the choice of origin.
See the documentation of MagneticMomentWriter
for parameters controlling the origin location.
The magnetic moment data can be written at multiple different origins during a single propagation as demonstrated in this script:
from gpaw import setup_paths
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.lcaotddft.magneticmomentwriter import MagneticMomentWriter
# Insert the path to the created basis set
setup_paths.insert(0, '.')
def main(kick):
kick_strength = [0., 0., 0.]
kick_strength['xyz'.index(kick)] = 1e-5
td_calc = LCAOTDDFT('gs.gpw', txt=f'td-{kick}.out')
DipoleMomentWriter(td_calc, f'dm-{kick}.dat')
# Origin: center of mass
MagneticMomentWriter(td_calc, f'mm-COM-{kick}.dat',
origin='COM')
# Origin: center of mass + 5 Å shift
for shift_axis in 'xyz':
origin_shift = [0, 0, 0]
origin_shift['xyz'.index(shift_axis)] = 5
MagneticMomentWriter(td_calc, f'mm-COM+{shift_axis}-{kick}.dat',
origin='COM', origin_shift=origin_shift)
# Origin: arbitrary coordinate
MagneticMomentWriter(td_calc, f'mm-123-{kick}.dat',
origin='zero', origin_shift=[1, 2, 3])
td_calc.absorption_kick(kick_strength)
td_calc.propagate(10, 2000)
td_calc.write(f'td-{kick}.gpw', mode='all')
if __name__ == '__main__':
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('--kick', default='x')
kwargs = vars(parser.parse_args())
main(**kwargs)
The resulting spectra: