Timepropagation TDDFT with LCAO¶
This page documents the use of timepropagation TDDFT in LCAO mode. The implementation is described in Ref. 1.
Realtime propagation of LCAO functions¶
In the realtime LCAOTDDFT approach, the timedependent wave functions are represented using the localized basis functions \(\tilde{\phi}_{\mu}(\mathbf r)\) as
The timedependent Kohn–Sham equation in the PAW formalism can be written as
From this, the following matrix equation can be derived for the LCAO wave function coefficients:
In the current implementation, \(\mathbf C\), \(\mathbf S\) and \(\mathbf H\) are dense matrices which are distributed using ScaLAPACK/BLACS. Currently, a semiimplicit Crank–Nicolson method (SICN) is used to propagate the wave functions. For wave functions at time \(t\), one propagates the system forward using \(\mathbf H(t)\) and solving the linear equation
Using the predicted wave functions \(C'(t+\mathrm dt)\), the Hamiltonian \(H'(t+\mathrm dt)\) is calculated and the Hamiltonian at middle of the time step is estimated as
With the improved Hamiltonian, the wave functions are again propagated from \(t\) to \(t+\mathrm dt\) by solving
This procedure is repeated using 500–2000 time steps of 5–40 as to calculate the time evolution of the electrons.
Example usage¶
First we do a standard groundstate calculation with the GPAW
calculator:
# Sodium dimer
from ase.build import molecule
atoms = molecule('Na2')
atoms.center(vacuum=6.0)
# Poisson solver with increased accuracy and multipole corrections up to l=2
from gpaw import PoissonSolver
poissonsolver = PoissonSolver(eps=1e20, remove_moment=1 + 3 + 5)
# Groundstate calculation
from gpaw import GPAW
calc = GPAW(mode='lcao', h=0.3, basis='dzp',
setups={'Na': '1'},
poissonsolver=poissonsolver,
convergence={'density': 1e8})
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
Some important points are:
The grid spacing is only used to calculate the Hamiltonian matrix and therefore a coarser grid than usual can be used.
Completely unoccupied bands should be left out of the calculation, since they are not needed.
The density convergence criterion should be a few orders of magnitude more accurate than in usual groundstate calculations.
The convergence tolerance of the Poisson solver should be at least
1e16
, but1e20
does not hurt (note that this is the quadratic error).One should use multipolecorrected Poisson solvers or other advanced Poisson solvers in any TDDFT run in order to guarantee the convergence of the potential with respect to the vacuum size. See the documentation on Advanced Poisson solvers.
Next we kick the system in the z direction and propagate 3000 steps of 10 as:
# Timepropagation calculation
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
# Read converged groundstate file
td_calc = LCAOTDDFT('gs.gpw')
# Attach any data recording or analysis tools
DipoleMomentWriter(td_calc, 'dm.dat')
# Kick
td_calc.absorption_kick([0.0, 0.0, 1e5])
# Propagate
td_calc.propagate(10, 3000)
# Save the state for restarting later
td_calc.write('td.gpw', mode='all')
After the time propagation, the spectrum can be calculated:
# Analyze the results
from gpaw.tddft.spectrum import photoabsorption_spectrum
photoabsorption_spectrum('dm.dat', 'spec.dat')
This example input script can be downloaded here
.
Extending the functionality¶
The realtime propagation LCAOTDDFT provides very modular framework for calculating many properties during the propagation. See Modular analysis tools for tutorial how to extend the analysis capabilities.
General notes about basis sets¶
In timepropagation LCAOTDDFT, it is much more important to think about the basis sets compared to groundstate LCAO calculations. It is required that the basis set can represent both the occupied (holes) and relevant unoccupied states (electrons) adequately. Custom basis sets for the time propagation should be generated according to one’s need, and then benchmarked.
Irrespective of the basis sets you choose always benchmark LCAO results with respect to the FD timepropagation code on the largest system possible. For example, one can create a prototype system which consists of similar atomic species with similar roles as in the parent system, but small enough to calculate with grid propagation mode. See Figs. 4 and 5 of Ref. 1 for example benchmarks.
After these remarks, we describe two packages of basis sets that can be used as a starting point for choosing a suitable basis set for your needs. Namely, pvalence basis sets and Completenessoptimized basis sets.
pvalence basis sets¶
The socalled pvalence basis sets are constructed for transition metals by replacing the ptype polarization function of the default basis sets with a bound unoccupied ptype orbital and its splitvalence complement. Such basis sets correspond to the ones used in Ref. 1. These basis sets significantly improve density of states of unoccupied states.
The pvalence basis sets can be easily obtained for the appropriate elements with the gpaw installdata tool using the following options:
$ gpaw installdata {<directory>} basis version=pvalence
See Installation of PAW datasets for more information on basis set installation. It is again reminded that these basis sets are not thoroughly tested and it is essential to benchmark the performance of the basis sets for your application.
Completenessoptimized basis sets¶
A systematic approach for improving the basis sets can be obtained with the socalled completenessoptimization approach. This approach is used in Ref. 2 to generate basis set series for TDDFT calculations of copper, silver, and gold clusters.
For further details of the basis sets, as well as their construction and performance, see 2. For convenience, these basis sets can be easily obtained with:
$ gpaw installdata {<directory>} basis version=coopt
See Installation of PAW datasets for basis set installation. Finally, it is again emphasized that when using the basis sets, it is essential to benchmark their suitability for your application.
Parallelization¶
LCAOTDDFT is parallelized using ScaLAPACK. It runs without ScaLAPACK, but in this case only a single core is used for linear alrebra.
Use
parallel={'sl_default':(N, M, 64)}
; See Parallelization options.It is necessary that N*M equals the total number of cores used by the calculator, and
max(N,M)*64 < nbands
, where64
is the used block size. The block size can be changed to, e.g., 16 if necessary.Apart from parallelization of linear algrebra, normal domain and band parallelizations can be used. As in groundstate LCAO calculations, use band parallelization to reduce memory consumption.
Modular analysis tools¶
In Example usage it was demonstrated how to calculate photoabsorption
spectrum from the timedependent dipole moment data collected with
DipoleMomentWriter
observer.
However, any (also userwritten) analysis tools can be attached
as a separate observers in the general timepropagation framework.
 There are two ways to perform analysis:
Perform analysis onthefly during the propagation:
# Read starting point td_calc = LCAOTDDFT('gs.gpw') # Attach analysis tools MyObserver(td_calc, ...) # Kick and propagate td_calc.absorption_kick([1e5, 0., 0.]) td_calc.propagate(10, 3000)
For example, the analysis tools can be
DipoleMomentWriter
observer for spectrum or Fourier transform of density at specific frequencies etc.Record the wave functions during the first propagation and perform the analysis later by replaying the propagation:
# Read starting point td_calc = LCAOTDDFT('gs.gpw') # Attach analysis tools MyObserver(td_calc, ...) # Replay propagation from a stored file td_calc.replay(name='wf.ulm', update='all')
From the perspective of the attached observers the replaying is identical to actual propagation.
The latter method is recommended, because one might not know beforehand what to analyze. For example, the interesting resonance frequencies are often not know before the timepropagation calculation.
In the following we give an example how to utilize the replaying capability in practice and describe some analysis tools available in GPAW.
Example¶
We use a finite sodium atom chain as an example system. First, let’s do the groundstate calculation:
from ase import Atoms
from gpaw import GPAW
from gpaw.poisson import PoissonSolver
from gpaw.poisson_extravacuum import ExtraVacuumPoissonSolver
# Sodium atom chain
atoms = Atoms('Na8',
positions=[[i * 3.0, 0, 0] for i in range(8)],
cell=[33.6, 12.0, 12.0])
atoms.center()
# Use an advanced Poisson solver
eps = 1e16
ps = ExtraVacuumPoissonSolver(gpts=(512, 256, 256),
poissonsolver_large=PoissonSolver(eps=eps),
coarses=2,
poissonsolver_small=PoissonSolver(eps=eps))
# Groundstate calculation
calc = GPAW(mode='lcao', h=0.3, basis='pvalence.dz', xc='LDA', nbands=6,
setups={'Na': '1'},
poissonsolver=ps,
convergence={'density': 1e8},
txt='gs.out')
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
Note the recommended use of Advanced Poisson solvers and pvalence basis sets.
Recording the wave functions and replaying the time propagation¶
We can record the timedependent wave functions during the propagation
with WaveFunctionWriter
observer:
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.lcaotddft.wfwriter import WaveFunctionWriter
# Read the groundstate file
td_calc = LCAOTDDFT('gs.gpw', txt='td.out')
# Attach any data recording or analysis tools
DipoleMomentWriter(td_calc, 'dm.dat')
WaveFunctionWriter(td_calc, 'wf.ulm')
# Kick and propagate
td_calc.absorption_kick([1e5, 0., 0.])
td_calc.propagate(20, 1500)
# Save the state for restarting later
td_calc.write('td.gpw', mode='all')
Tip
The time propagation can be continued in the same manner from the restart file:
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.lcaotddft.wfwriter import WaveFunctionWriter
# Read the restart file
td_calc = LCAOTDDFT('td.gpw', txt='tdc.out')
# Attach the data recording for appending the new data
DipoleMomentWriter(td_calc, 'dm.dat')
WaveFunctionWriter(td_calc, 'wf.ulm')
# Continue propagation
td_calc.propagate(20, 500)
# Save the state for restarting later
td_calc.write('td.gpw', mode='all')
The created wf.ulm
file contains the timedependent wave functions
\(C_{\mu n}(t)\) that define the state of the system at each time.
We can use the file to replay the time propagation:
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
# Read the groundstate file
td_calc = LCAOTDDFT('gs.gpw')
# Attach analysis tools
DipoleMomentWriter(td_calc, 'dm_replayed.dat')
# Replay the propagation
td_calc.replay(name='wf.ulm', update='all')
The update
keyword in replay()
has following options:

variables updated during replay 


Hamiltonian and density 

density 

nothing 
Tip
The wave functions can be written in separate files
by using split=True
:
WaveFunctionWriter(td_calc, 'wf.ulm', split=True)
This creates additional wf*.ulm
files containing the wave functions.
The replay functionality works as in the above example
even with splitted files.
Kohn–Sham decomposition of density matrix¶
Kohn–Sham decomposition is an illustrative way of analyzing electronic excitations in Kohn–Sham electronhole basis. See Ref. 3 for the description of the GPAW implementation and demonstration.
Here we demonstrate how to construct the photoabsorption decomposition at a specific frequency in Kohn–Sham electonhole basis.
First, let’s calculate and plot
the spectrum:
from gpaw.tddft.spectrum import photoabsorption_spectrum
photoabsorption_spectrum('dm.dat', 'spec.dat',
folding='Gauss', width=0.1,
e_min=0.0, e_max=10.0, delta_e=0.01)
The two main resonances are analyzed in the following.
Frequencyspace density matrix¶
We generate the density matrix for the frequencies of interest:
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.densitymatrix import DensityMatrix
from gpaw.lcaotddft.frequencydensitymatrix import FrequencyDensityMatrix
from gpaw.tddft.folding import frequencies
# Read the groundstate file
td_calc = LCAOTDDFT('gs.gpw')
# Attach analysis tools
dmat = DensityMatrix(td_calc)
freqs = frequencies([1.12, 2.48], 'Gauss', 0.1)
fdm = FrequencyDensityMatrix(td_calc, dmat, frequencies=freqs)
# Replay the propagation
td_calc.replay(name='wf.ulm', update='none')
# Store the density matrix
fdm.write('fdm.ulm')
Transform the density matrix to Kohn–Sham electronhole basis¶
First, we construct the Kohn–Sham electronhole basis. For that we need to calculate unoccupied Kohn–Sham states, which is conveniently done by restarting from the earlier groundstate file:
from gpaw import GPAW
from gpaw.lcaotddft.ksdecomposition import KohnShamDecomposition
# Calculate ground state with full unoccupied space
calc = GPAW('gs.gpw', nbands='nao', fixdensity=True,
txt='unocc.out')
atoms = calc.get_atoms()
energy = atoms.get_potential_energy()
calc.write('unocc.gpw', mode='all')
# Construct KS electronhole basis
ksd = KohnShamDecomposition(calc)
ksd.initialize(calc)
ksd.write('ksd.ulm')
Next, we can use the created objects to transform the LCAO density matrix to the Kohn–Sham electronhole basis and visualize the photoabsorption decomposition as a transition contribution map (TCM):
# Creates: tcm_1.12.png, tcm_2.48.png, table_1.12.txt, table_2.48.txt
import numpy as np
from matplotlib import pyplot as plt
from gpaw import GPAW
from gpaw.tddft.units import au_to_eV
from gpaw.lcaotddft.ksdecomposition import KohnShamDecomposition
from gpaw.lcaotddft.densitymatrix import DensityMatrix
from gpaw.lcaotddft.frequencydensitymatrix import FrequencyDensityMatrix
# Load the objects
calc = GPAW('unocc.gpw', txt=None)
ksd = KohnShamDecomposition(calc, 'ksd.ulm')
dmat = DensityMatrix(calc)
fdm = FrequencyDensityMatrix(calc, dmat, 'fdm.ulm')
def do(w):
# Select the frequency and the density matrix
rho_uMM = fdm.FReDrho_wuMM[w]
freq = fdm.freq_w[w]
frequency = freq.freq * au_to_eV
print('Frequency: %.2f eV' % frequency)
print('Folding: %s' % freq.folding)
# Transform the LCAO density matrix to KS basis
rho_up = ksd.transform(rho_uMM)
# Photoabsorption decomposition
dmrho_vp = ksd.get_dipole_moment_contributions(rho_up)
weight_p = 2 * freq.freq / np.pi * dmrho_vp[0].imag / au_to_eV * 1e5
print('Total absorption: %.2f eV^1' % np.sum(weight_p))
# Print contributions as a table
table = ksd.get_contributions_table(weight_p, minweight=0.1)
print(table)
with open('table_%.2f.txt' % frequency, 'w') as f:
f.write('Frequency: %.2f eV\n' % frequency)
f.write('Folding: %s\n' % freq.folding)
f.write('Total absorption: %.2f eV^1\n' % np.sum(weight_p))
f.write(table)
# Plot the decomposition as a TCM
r = ksd.plot_TCM(weight_p,
occ_energy_min=3, occ_energy_max=0.1,
unocc_energy_min=0.1, unocc_energy_max=3,
delta_energy=0.01, sigma=0.1
)
(ax_tcm, ax_occ_dos, ax_unocc_dos, ax_spec) = r
# Plot diagonal line at the analysis frequency
x = np.array([4, 1])
y = x + freq.freq * au_to_eV
ax_tcm.plot(x, y, c='k')
ax_occ_dos.set_title('Photoabsorption TCM of Na8 at %.2f eV' % frequency)
# Save the plot
plt.savefig('tcm_%.2f.png' % frequency)
do(0)
do(1)
Note that the sum over the decomposition (the printed total absorption) equals to the value of the photoabsorption spectrum at the particular frequency in question.
We obtain the following transition contributions for the resonances (presented both as tables and TCMs):
Frequency: 1.12 eV
Folding: Gauss(0.10000 eV)
Total absorption: 30.11 eV^1
# p i( eV) a( eV) Ediff (eV) weight
0 3( 0.214) > 4( 0.214): 0.4279 +29.2952
6 3( 0.214) > 6( 1.434): 1.6479 +0.5588
36 3( 0.214) > 16( 2.456): 2.6696 +0.1294
5 2( 0.593) > 5( 0.801): 1.3943 +0.1198
rest: +0.0075
total: +30.1106
Frequency: 2.48 eV
Folding: Gauss(0.10000 eV)
Total absorption: 1.53 eV^1
# p i( eV) a( eV) Ediff (eV) weight
6 3( 0.214) > 6( 1.434): 1.6479 +1.0450
0 3( 0.214) > 4( 0.214): 0.4279 0.5812
5 2( 0.593) > 5( 0.801): 1.3943 +0.4450
3 1( 0.871) > 4( 0.214): 1.0853 +0.3866
rest: +0.2332
total: +1.5286
Induced density¶
The density matrix gives access to any other quantities. For instance, the induced density can be conveniently obtained from the density matrix:
import numpy as np
from ase.io import write
from gpaw import GPAW
from gpaw.tddft.units import au_to_eV
from gpaw.lcaotddft.densitymatrix import DensityMatrix
from gpaw.lcaotddft.frequencydensitymatrix import FrequencyDensityMatrix
# Load the objects
calc = GPAW('unocc.gpw', txt=None)
calc.initialize_positions() # Initialize in order to calculate density
dmat = DensityMatrix(calc)
fdm = FrequencyDensityMatrix(calc, dmat, 'fdm.ulm')
def do(w):
# Select the frequency and the density matrix
rho_MM = fdm.FReDrho_wuMM[w][0]
freq = fdm.freq_w[w]
frequency = freq.freq * au_to_eV
print('Frequency: %.2f eV' % frequency)
print('Folding: %s' % freq.folding)
# Induced density
rho_g = dmat.get_density([rho_MM.imag])
# Save as a cube file
write('ind_%.2f.cube' % frequency, calc.atoms, data=rho_g)
# Calculate dipole moment for reference
dm_v = dmat.density.finegd.calculate_dipole_moment(rho_g, center=True)
absorption = 2 * freq.freq / np.pi * dm_v[0] / au_to_eV * 1e5
print('Total absorption: %.2f eV^1' % absorption)
do(0)
do(1)
The resulting cube files can be visualized, for example, with
this script
producing
the figures:
Advanced tutorials¶
Plasmon resonance of silver cluster¶
One should think about what type of transitions of interest are present, and make sure that the basis set can represent such KohnSham electron and hole wave functions. The first transitions in silver clusters will be \(5s \rightarrow 5p\) like. We require 5p orbitals in the basis set, and thus, we must generate a custom basis set.
Here is how to generate a doublezeta basis set with 5p orbital in valence for silver for GLLBSC potential. Note that the extra 5p valence state effectively improves on the ordinary polarization function, so this basis set is better than the default doublezeta polarized one. We will use the 11electron Ag setup, since the semicore p states included in the default setup are not relevant here.
from gpaw.atom.generator import Generator
from gpaw.atom.basis import BasisMaker
args = {'core': '[Kr]', 'rcut': 2.45}
generator = Generator('Ag', 'GLLBSC', scalarrel=True)
generator.N *= 2 # Increase grid resolution
generator.run(**args)
bm = BasisMaker(generator, name='GLLBSC', run=False)
basis = bm.generate(zetacount=2, polarizationcount=0,
energysplit=0.07,
jvalues=[0, 1, 2], # include d, s and p
rcutmax=12.0)
basis.write_xml()
We calculate the icosahedral Ag55 cluster: ag55.xyz
This code uses ScaLAPACK parallelization with 48 cores.
from ase.io import read
from gpaw import GPAW
from gpaw import PoissonSolver
from gpaw.occupations import FermiDirac
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.tddft.spectrum import photoabsorption_spectrum
from gpaw import setup_paths
# Set the path for the created basis set
setup_paths.insert(0, '.')
# Read the clusterm from the xyz file
atoms = read('ag55.xyz')
atoms.center(vacuum=5.0)
# Increase the accuracy of density for ground state
convergence = {'density': 1e8}
# Increase the accuracy of PoissonSolver and
# apply multipole corrections for monopole and dipoles
poissonsolver = PoissonSolver(eps=1e16, remove_moment=1 + 3)
# Calculate ground state in LCAO mode
calc = GPAW(xc='GLLBSC', basis='GLLBSC.dz', h=0.3, nbands=352, mode='lcao',
convergence=convergence, poissonsolver=poissonsolver,
occupations=FermiDirac(0.1),
parallel={'sl_default': (8, 6, 32), 'band': 2})
atoms.set_calculator(calc)
# Relax the ground state
atoms.get_potential_energy()
# Save the intermediate ground state result to a file
calc.write('ag55_gs.gpw', mode='all')
# Time propagation
td_calc = LCAOTDDFT('ag55_gs.gpw',
parallel={'sl_default': (8, 6, 32), 'band': 2})
DipoleMomentWriter(td_calc, 'ag55.dm')
td_calc.absorption_kick([1e5, 0.0, 0.0])
td_calc.propagate(20, 500)
photoabsorption_spectrum('ag55.dm', 'ag55.spec', width=0.2)
Code runs for approximately two hours wallclock time. The resulting spectrum shows already emerging plasmonic excitation around 4 eV. For more details, see 1.
References¶
 1(1,2,3,4)
M. Kuisma, A. Sakko, T. P. Rossi, A. H. Larsen, J. Enkovaara, L. Lehtovaara, and T. T. Rantala, Localized surface plasmon resonance in silver nanoparticles: Atomistic firstprinciples timedependent density functional theory calculations, Phys. Rev. B 69, 245419 (2004). doi:10.1103/PhysRevB.91.115431
 2(1,2)
T. P. Rossi, S. Lehtola, A. Sakko, M. J. Puska, and R. M. Nieminen, Nanoplasmonics simulations at the basis set limit through completenessoptimized, local numerical basis sets, J. Chem. Phys. 142, 094114 (2015). doi:10.1063/1.4913739
 3
T. P. Rossi, M. Kuisma, M. J. Puska, R. M. Nieminen, and P. Erhart, Kohn–Sham Decomposition in RealTime TimeDependent DensityFunctional Theory: An Efficient Tool for Analyzing Plasmonic Excitations, J. Chem. Theory Comput. 13, 4779 (2017). doi:10.1021/acs.jctc.7b00589