Time-propagation TDDFT with LCAO

This page documents the use of time-propagation TDDFT in LCAO mode. The implementation is described in Ref. [1].

Real-time propagation of LCAO functions

In the real-time LCAO-TDDFT approach, the time-dependent wave functions are represented using the localized basis functions \(\tilde{\phi}_{\mu}(\mathbf r)\) as

\[\tilde{\psi}_n(\mathbf{r},t) = \sum_{\mu} \tilde{\phi}_{\mu}(\mathbf{r}-\mathbf{R}^\mu) c_{\mu n}(t) .\]

The time-dependent Kohn–Sham equation in the PAW formalism can be written as

\[\left[ \widehat T^\dagger \left( -i \frac{{\rm d}}{{\rm d}t} + \hat H_{\rm KS}(t) \right) \widehat T \right] \tilde{\Psi(\mathbf{r},t)} = 0.\]

From this, the following matrix equation can be derived for the LCAO wave function coefficients:

\[{\rm i}\mathbf{S} \frac{{\rm d}\mathbf{C}(t)}{{\rm d}t} = \mathbf{H}(t) \mathbf{C}(t).\]

In the current implementation, \(\mathbf C\), \(\mathbf S\) and \(\mathbf H\) are dense matrices which are distributed using ScaLAPACK/BLACS. Currently, a semi-implicit 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

\[\left( \mathbf{S} + {\rm i} \mathbf{H}(t) {\rm d}t / 2 \right) \mathbf{C}'(t+{\rm d}t) = \left( \mathbf{S} - {\rm i} \mathbf{H}(t) {\rm d}t / 2 \right) \mathbf{C}(t).\]

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

\[\mathbf{H}(t+{\rm d}t/2) = (\mathbf{H}(t) + \mathbf{H}'(t+{\rm d}t)) / 2\]

With the improved Hamiltonian, the wave functions are again propagated from \(t\) to \(t+\mathrm dt\) by solving

\[\left( \mathbf{S} + {\rm i} \mathbf{H}(t+{\rm d}t/2) {\rm d}t / 2 \right) \mathbf{C}(t+{\rm d}t) = \left( \mathbf{S} - {\rm i} \mathbf{H}(t+{\rm d}t/2) {\rm d}t / 2 \right) \mathbf{C}(t).\]

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 ground-state 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=1e-20, remove_moment=1 + 3 + 5)

# Ground-state calculation
from gpaw import GPAW
calc = GPAW(mode='lcao', h=0.3, basis='dzp',
            setups={'Na': '1'},
            poissonsolver=poissonsolver,
            convergence={'density': 1e-8})
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 ground-state calculations.
  • The convergence tolerance of the Poisson solver should be at least 1e-16, but 1e-20 does not hurt (note that this is the quadratic error).
  • One should use multipole-corrected 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:

# Time-propagation calculation
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
# Read converged ground-state 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, 1e-5])
# 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 real-time 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 time-propagation LCAO-TDDFT, it is much more important to think about the basis sets compared to ground-state 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 time-propagation 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, p-valence basis sets and Completeness-optimized basis sets.

p-valence basis sets

The so-called p-valence basis sets are constructed for transition metals by replacing the p-type polarization function of the default basis sets with a bound unoccupied p-type orbital and its split-valence complement. Such basis sets correspond to the ones used in Ref. [1]. These basis sets significantly improve density of states of unoccupied states.

The p-valence basis sets can be easily obtained for the appropriate elements with the gpaw install-data tool using the following options:

$ gpaw install-data {<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.

Completeness-optimized basis sets

A systematic approach for improving the basis sets can be obtained with the so-called completeness-optimization 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 install-data {<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

LCAO-TDDFT 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, where 64 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 ground-state 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 time-dependent dipole moment data collected with DipoleMomentWriter observer. However, any (also user-written) analysis tools can be attached as a separate observers in the general time-propagation framework.

There are two ways to perform analysis:
  1. Perform analysis on-the-fly during the propagation:

    # Read starting point
    td_calc = LCAOTDDFT('gs.gpw')
    
    # Attach analysis tools
    MyObserver(td_calc, ...)
    
    # Kick and propagate
    td_calc.absorption_kick([1e-5, 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.

  2. 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 time-propagation 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 ground-state 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 = 1e-16
ps = ExtraVacuumPoissonSolver(gpts=(512, 256, 256),
                              poissonsolver_large=PoissonSolver(eps=eps),
                              coarses=2,
                              poissonsolver_small=PoissonSolver(eps=eps))

# Ground-state calculation
calc = GPAW(mode='lcao', h=0.3, basis='pvalence.dz', xc='LDA', nbands=6,
            setups={'Na': '1'},
            poissonsolver=ps,
            convergence={'density': 1e-8},
            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 p-valence basis sets.

Recording the wave functions and replaying the time propagation

We can record the time-dependent 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 ground-state 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([1e-5, 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 time-dependent 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 ground-state 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:

update variables updated during replay
'all' Hamiltonian and density
'density' density
'none' 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 electron-hole 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 electon-hole 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)
../../_images/spec.png

The two main resonances are analyzed in the following.

Frequency-space 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 ground-state 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 electron-hole basis

First, we construct the Kohn–Sham electron-hole basis. For that we need to calculate unoccupied Kohn–Sham states, which is conveniently done by restarting from the earlier ground-state 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 electron-hole 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 electron-hole 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
../../_images/tcm_1.12.png
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
../../_images/tcm_2.48.png

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:

../../_images/ind_1.12.png ../../_images/ind_2.48.png

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 Kohn-Sham 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 double-zeta basis set with 5p orbital in valence for silver for GLLB-SC potential. Note that the extra 5p valence state effectively improves on the ordinary polarization function, so this basis set is better than the default double-zeta polarized one. We will use the 11-electron Ag setup, since the semi-core 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': 1e-8}

# Increase the accuracy of PoissonSolver and
# apply multipole corrections for monopole and dipoles
poissonsolver = PoissonSolver(eps=1e-16, 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([1e-5, 0.0, 0.0])
td_calc.propagate(20, 500)

photoabsorption_spectrum('ag55.dm', 'ag55.spec', width=0.2)

Code runs for approximately two hours wall-clock time. The resulting spectrum shows already emerging plasmonic excitation around 4 eV. For more details, see [1].

../../_images/fig1.png

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 first-principles time-dependent 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 completeness-optimized, 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 Real-Time Time-Dependent Density-Functional Theory: An Efficient Tool for Analyzing Plasmonic Excitations, J. Chem. Theory Comput. 13, 4779 (2017). doi:10.1021/acs.jctc.7b00589