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': 1e12})
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¶
For maximum performance on large systems, it is advicable to use
ScaLAPACK and large band parallelization with augment_grids
enabled.
This can be achieved with parallelization settings like
parallel={'sl_auto': True, 'domain': 8, 'augment_grids': True}
(see Parallelization options),
which will use groups of 8 tasks for domain parallelization and the rest for
band parallelization (for example, with a total of 144 cores this would mean
domain and band parallelizations of 8 and 18, respectively).
Instead of sl_auto
, the ScaLAPACK settings can be set by hand
as sl_default=(m, n, block)
(see ScaLAPACK,
in which case it is important that m * n`
equals
the total number of cores used by the calculator
and that max(m, n) * block < nbands
.
It is also possible to run the code without ScaLAPACK, but it is very inefficient for large systems as in that case only a single core is used for linear algebra.
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': 1e12},
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¶
In this tutorial, we demonstrate the use of
efficient parallelization settings and
calculate the photoabsorption spectrum of
an icosahedral Ag55 cluster
.
We use GLLBSC potential to significantly improve the description of d states,
pvalence basis sets to improve the description of unoccupied states, and
11electron Ag setup to reduce computational cost.
When calculating other systems, remember to check the convergence with respect to the used basis sets. Recall hints here.
The workflow is the same as in the previous examples. First, we calculate ground state (takes around 20 minutes with 36 cores):
from ase.io import read
from gpaw import GPAW, FermiDirac, Mixer, PoissonSolver
# Read the structure from the xyz file
atoms = read('Ag55.xyz')
atoms.center(vacuum=6.0)
# Increase the accuracy of density for ground state
convergence = {'density': 1e12}
# Use occupation smearing and weak mixer to facilitate convergence
occupations = FermiDirac(25e3)
mixer = Mixer(0.02, 5, 1.0)
# Parallelzation settings
parallel = {'sl_auto': True, 'domain': 2, 'augment_grids': True}
# Increase the accuracy of PoissonSolver and
# apply multipole corrections for monopole and dipoles
poissonsolver = PoissonSolver(eps=1e16, remove_moment=1 + 3)
# Groundstate calculation
calc = GPAW(mode='lcao', xc='GLLBSC', h=0.3, nbands=360,
setups={'Ag': '11'},
basis={'Ag': 'pvalence.dz', 'default': 'dpz'},
convergence=convergence, poissonsolver=poissonsolver,
occupations=occupations, mixer=mixer, parallel=parallel,
maxiter=1000,
txt='gs.out')
atoms.set_calculator(calc)
atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
Then, we carry the time propagation for 30 femtoseconds in steps of 10 attoseconds (takes around 3.5 hours with 36 cores):
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
# Parallelzation settings
parallel = {'sl_auto': True, 'domain': 2, 'augment_grids': True}
# Time propagation
td_calc = LCAOTDDFT('gs.gpw', parallel=parallel, txt='td.out')
DipoleMomentWriter(td_calc, 'dm.dat')
td_calc.absorption_kick([1e5, 0.0, 0.0])
td_calc.propagate(10, 3000)
Finally, we calculate the spectrum (takes a few seconds):
from gpaw.tddft.spectrum import photoabsorption_spectrum
# Calculate spectrum
photoabsorption_spectrum('dm.dat', 'spec.dat', width=0.1, delta_e=0.01)
The resulting spectrum shows an emerging plasmon resonance at around 4 eV:
For further insight on plasmon resonance in metal nanoparticles, see 1 and 3.
Usergenerated basis sets¶
The pvalence basis sets distributed with GPAW and used in the above tutorial have been generated from atomic PBE orbitals. Similar basis sets can be generated based on atomic GLLBSC orbitals:
from gpaw.atom.generator import Generator
from gpaw.atom.basis import BasisMaker
from gpaw.atom.configurations import parameters, parameters_extra
xc = 'GLLBSC'
name = 'my'
if 'Ag' in parameters_extra:
args = parameters_extra['Ag'] # Choose the 11electron setup
else:
args = parameters['Ag'] # Choose the 17electron setup
args.update(dict(name=name, use_restart_file=False, exx=True))
# Generate setup
generator = Generator('Ag', xc, scalarrel=True)
generator.run(write_xml=True, **args)
# Generate basis
bm = BasisMaker('Ag', name='{}.{}'.format(name, xc), xc=xc, run=False)
bm.generator.run(write_xml=False, **args)
basis = bm.generate(zetacount=2, polarizationcount=0,
jvalues=[0, 1, 2], # include d, s and p
)
basis.write_xml()
The Ag55 cluster can be calculated as in the above tutorial, once the input scripts have been modified to use the generated setup and basis set. Changes to the groundstate script:
 /tmp/gpawdocswebpage/gpaw/doc/documentation/tddft/lcaotddft_Ag55/gs.py
+++ /tmp/gpawdocswebpage/gpaw/doc/documentation/tddft/lcaotddft_Ag55/mybasis/gs.py
@@ 1,5 +1,9 @@
from ase.io import read
from gpaw import GPAW, FermiDirac, Mixer, PoissonSolver
+from gpaw import setup_paths
+
+# Insert the path to the created basis set
+setup_paths.insert(0, '.')
# Read the structure from the xyz file
atoms = read('Ag55.xyz')
@@ 21,8 +25,8 @@
# Groundstate calculation
calc = GPAW(mode='lcao', xc='GLLBSC', h=0.3, nbands=360,
 setups={'Ag': '11'},
 basis={'Ag': 'pvalence.dz', 'default': 'dpz'},
+ setups={'Ag': 'my'},
+ basis={'Ag': 'GLLBSC.dz', 'default': 'dpz'},
convergence=convergence, poissonsolver=poissonsolver,
occupations=occupations, mixer=mixer, parallel=parallel,
maxiter=1000,
Changes to the timepropagation script:
 /tmp/gpawdocswebpage/gpaw/doc/documentation/tddft/lcaotddft_Ag55/td.py
+++ /tmp/gpawdocswebpage/gpaw/doc/documentation/tddft/lcaotddft_Ag55/mybasis/td.py
@@ 1,5 +1,9 @@
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
+from gpaw import setup_paths
+
+# Insert the path to the created basis set
+setup_paths.insert(0, '.')
# Parallelzation settings
parallel = {'sl_auto': True, 'domain': 2, 'augment_grids': True}
The calculation with this generated “my” pvalence basis set results only in small differences in the spectrum in comparison to the distributed pvalence basis sets:
The spectrum with the default dzp basis sets is also shown for reference, resulting in unconverged spectrum due to the lack of diffuse functions. This demonstrates the importance of checking convergence with respect to the used basis sets. Recall hints here and see 1 and 2 for further discussion on the basis sets.
References¶
 1(1,2,3,4,5)
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,3)
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(1,2)
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