rhodent: A Python package for analyzing real-time TDDFT response

TL;DR

rhodent is a Python package that simplifies the analysis of real-time TDDFT simulation outputs, enabling efficient computation of spectroscopic and dynamic properties across different codes.

Contribution

It introduces a modular Python tool for analyzing rt-TDDFT results, compatible with multiple codes and capable of various spectral and dynamic analyses.

Findings

Supports calculation of hot-carrier distributions and energies

Handles analysis for Al, Ag clusters, and organic molecules

Can compute responses to laser perturbations efficiently

Abstract

Real-time time-dependent density functional theory (rt-TDDFT) is a well-established method for studying the dynamic response of matter in the femtosecond or optical range. In this method, the Kohn-Sham (KS) wave functions are propagated forward in time, and in principle, one can extract any observable at any given time. Alternatively, by taking a Fourier transform, spectroscopic quantities can be extracted. There are many publicly available codes implementing rt-TDDFT, which differ in their numeric solution of the KS equations, their available exchange-correlation functionals, and in their analysis capabilities. For users of rt-TDDFT, this is an inconvenient situation because they may need to use a numerical method that is available in one code, but an analysis method available in another. Here, we introduce rhodent, a modular Python package for processing the output of rt-TDDFT calculations. Our package can be used to calculate hot-carrier distributions, energies, induced densities, and dipole moments, and various decompositions thereof. In its current version, rhodent handles calculation results from the gpaw code, but can readily be extended to support other rt-TDDFT codes. Additionally, under the assumption of linear response, rhodent can be used to calculate the response to a narrow-band laser, from the response to a broad-band perturbation, greatly speeding up the analysis of frequency-dependent excitations. We demonstrate the capabilities of rhodent via a set of examples, for systems consisting of Al and Ag clusters and organic molecules.

Equations

$$\begin{array}{rl}& {\rho }_{n{n}^{\prime }}(t)=\\ & \sum _k{f}_k\int {\psi }_{n^{\prime }}^{\ast }(\mathbf{r},0){\psi }_k(\mathbf{r},t)\mathrm{d}\mathbf{r}\int {\psi }_k^{\ast }({\mathbf{r}}^{\prime },t){\psi }_n({\mathbf{r}}^{\prime },0)\mathrm{d}{\mathbf{r}}^{\prime },\end{array}$$

$$\frac{\delta {\rho }_{n{n}^{\prime }}(\omega )}{v(\omega )}=\frac{{\int }_0^{\infty }\delta {\rho }_{n{n}^{\prime }}(t)e^{i\omega t}\mathrm{d}t}{{\int }_0^{\infty }v(t)e^{i\omega t}\mathrm{d}t},$$

$$\delta {\rho }_{n{n}^{\prime }}(\omega )$$

$$v(\omega )$$

$$\delta {\rho }_{n{n}^{\prime }}(\omega ;\sigma )$$

$$\delta {\rho }_{ia}(t)={\int }_0^t{\hat{\chi }}_{ia}(t-t^{\prime })\hat{v}(t^{\prime })\mathrm{d}{t}^{\prime },$$

$$\delta {\rho }_{ia}(\omega )={\hat{\chi }}_{ia}(\omega )\hat{v}(\omega ),$$

$$\delta {\rho }_{ia}^{\prime }(\omega )$$

$$=\delta {\rho }_{ia}(\omega )\frac{v^{\prime }(\omega )}{v(\omega )}.$$

$$\delta {\hat{\rho }}_{k,ia}$$

$${\hat{v}}_k$$

$${\hat{v}}_k^{\prime }$$

$$\delta {\rho }_{ia}^{\prime }(t_j)=\frac{1}{N^{\prime }}\sum _{k=0}^{N-1}\frac{\delta {\hat{\rho }}_{k,ia}\cdot {\hat{v}}_k^{\prime }}{{\hat{v}}_k}{e}^{-i{\omega }_k{t}_j}.$$

$$v(t)=s_0\frac{\sin ({\omega }_{cut}(t-t_0))}{{\omega }_{cut}(t-t_0)},t>0$$

$$v(\omega )=s_0\frac{e^{i\omega t_0}}{{\omega }_{\text{cut}}}\text{rect}\left(\frac{\omega }{2\pi {\omega }_{\text{cut}}}\right),$$

$${\psi }_n(\mathbf{r},t)=\sum _{\mu }{C}_{n\mu }(t){\Phi }_{\mu }(\mathbf{r}),$$

$${\rho }_{\mu \nu }(t)=\sum _k{C}_{k\nu }^{\ast }(t)f_k{C}_{k\nu }(t).$$

$$\delta {\rho }_{n{n}^{\prime }}(t)=\sum _{\mu \nu }{P}_{n^{\prime }\nu }^{\ast }\left[{\rho }_{\mu \nu }(t)-{\rho }_{\mu \nu }(0)\right]P_{n\mu }.$$

$${\rho }_{n{n}^{\prime }}(\omega )=\sum _{\mu \nu }{P}_{n^{\prime }\nu }^{\ast }{\rho }_{\mu \nu }(\omega )P_{n\mu }.$$

$$w_{n{n}^{\prime }}^{\text{proj}}={\int }_{\text{Voronoi cell}}{\psi }_{n^{\prime }}^{\ast }(\mathbf{r},0){\psi }_n(\mathbf{r},0)\mathrm{d}\mathbf{r}.$$

$$\delta n(\mathbf{r})=-\sum _{i{a}^{\prime }}{\psi }_a^{\ast }(\mathbf{r},0){\psi }_i(\mathbf{r},0)\delta {\rho }_{ia}.$$

$$\delta n(\mathbf{r})=-2\sum _{ia}^{f_i>f_a}{\psi }_a(\mathbf{r},0){\psi }_i(\mathbf{r},0)\mathrm{Re}\delta {\rho }_{ia}.$$

$$\delta n_{\text{filter}}(\mathbf{r})=-2\sum _{ia}^{\begin{array}{c}{f}_i>f_a\\ {\epsilon }_a>{\epsilon }_{\text{low}}\end{array}}{\psi }_a(\mathbf{r},0){\psi }_i(\mathbf{r},0)\mathrm{Re}\delta {\rho }_{ia}.$$

$$\delta \mathbf{\mu }=-2\sum _{ia}^{f_i>f_a}{\mathbf{\mu }}_{ia}\mathrm{Re}\delta {\rho }_{ia},$$

$${\mathbf{\mu }}_{ia}=\int {\psi }_a^{\ast }(\mathbf{r},0)\mathbf{r}{\psi }_i(\mathbf{r},0)\mathrm{d}\mathbf{r}.$$

$$\delta {\mathbf{\mu }}^{\text{occ}\to \text{unocc},ia}=-2\sum _{ia}^{f_i>f_a}{\mathbf{\mu }}_{ia}\mathrm{Re}\delta {\rho }_{ia}{w}_{ii}^{\text{proj,occ}}{w}_{aa}^{\text{proj,unocc}},$$

$$-2\sum _{ia}^{f_i>f_a}$$

$$\begin{array}{rl}-2\sum _{ia}^{f_i>f_a}& {\mathbf{\mu }}_{ia}\mathrm{Re}\delta {\rho }_{ia}{w}_{ii}^{\text{proj,occ}}{w}^{\text{proj,unocc}}\\ & G({\epsilon }_{\text{occ}}-{\epsilon }_i)G({\epsilon }_{\text{unocc}}-{\epsilon }_a).\end{array}$$

$$G(\epsilon -{\epsilon }_n)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{{\left({\epsilon }_i-\epsilon \right)}^2}{2{\sigma }^2}\right).$$

$$\hat{v}(t)=v(t)x,$$

$${\mathbf{\alpha }}_x(\omega )=\frac{{\int }_0^{\infty }\delta \mathbf{\mu }(t)e^{i\omega t}\mathrm{d}t}{v(\omega )}.$$

$$\mathbf{\alpha }(\omega )=\left[\begin{array}{ccc}{\mathbf{\alpha }}_x(\omega )& {\mathbf{\alpha }}_y(\omega )& {\mathbf{\alpha }}_z(\omega )\end{array}\right].$$

Full text

I Introduction

Real-time time-dependent density functional theory (rt-TDDFT) is a well-established method for studying the response of electronic systems [1, 2]. It has been used for spectroscopy calculations [3, 4, 5, 6, 7, 8, 9, 10], including circular dichroism [11, 12], X-rays [13], non-linear spectroscopy [14], and high-harmonic generation [15, 16, 17], as well as to study electron [18, 19, 20, 21, 22, 23] and spin dynamics [24] on femtosecond timescales. In essence, one simply takes the ground state calculated with Kohn-Sham (KS)-density functional theory (DFT) [25, 26] and numerically propagates the single-particle Schrödinger equation forward in time, under the influence of something driving the system out of the ground state, typically an external electric field. There are various codes implementing rt-TDDFT, including gpaw [27], octopus [28, 29], turbomole [30], Siesta [31], and cp2k [32]. These codes differ in implementation details and numerical representations of the KS wave functions.

There is, however, a need for modular libraries [33]. It may be desirable for some users to employ one code, due to its numerical implementation details, while using the analysis capabilities of another. Without modular libraries, the only solution to this problem is for the developer to either port the analysis capabilities or the numerical implementation to the other code. Here, we introduce rhodent – a Python package for processing the response from rt-TDDFT calculations. In its current version, calculation outputs from gpaw are supported, but the modular design of rhodent prepares it for future integration with other programs. rhodent is open source and comes with proper unit- and integration tests, easing development, even for external contributors. An extensive documentation is available online [34].

II Program overview

The rhodent package is written in Python and designed in a modular and object-oriented fashion. The key components are the Response and Calculator objects. The Response class allows reading output files generated in a previous rt-TDDFT calculation (outside rhodent), and transforms the response into the appropriate form for rhodent, either in the time or frequency domain. There are several implemented Response classes, for use with different rt-TDDFT output files. The Calculator is responsible for computing the observables at various times in the simulation, or the Fourier transform thereof. Several Calculator classes are available.

The modular structure of rhodent opens up the possibility to, with relatively little effort, extend the functionality of the code. In order to interface rhodent to another rt-TDDFT code than gpaw, which is currently supported, a new Response class needs to be added. Similarly, in order to compute a new type of observable, one in principle only needs to implement a new Calculator.

III Implementation and theory

The central quantity that rhodent works with is the time-dependent KS density matrix in the basis of ground-state KS orbitals

[TABLE]

where $\psi_{k}(\boldsymbol{r},t)$ are the time-dependent KS wave functions and $f_{k}$ their occupation numbers (including a spin degeneracy of 2). In particular, we are interested in the induced density matrix $\delta\rho_{nn^{\prime}}(t)=\rho_{nn^{\prime}}(t)-\rho_{nn^{\prime}}(0)$. This quantity contains all information about the system response to some time-dependent perturbation $\hat{v}(t)$, and is obtained from an rt-TDDFT calculation (e.g., using gpaw). We note here that, since rhodent currently is implemented for non-periodic and spin-paired systems, the only quantum number of the KS orbitals is $n$. We are restricted to fixed-atom calculations, so the time-dependence comes only from electronic degrees of freedom.

III.1 Obtaining the response

In rhodent, we can work with the response in the time domain, which is obtained directly from the rt-TDDFT calculation, or in the frequency domain. For the latter, we define the normalized Fourier transform

[TABLE]

where $v(t)$ is the scalar amplitude of the perturbation, and the lower bound on the integrals can be taken to be zero, because the integrands are zero before this time. This quantity is related to the Casida eigenvectors and gives similar information as the solution of the Casida equation [35].

III.1.1 In the frequency domain

In practice, the rt-TDDFT calculation results in samples of $\delta\rho(t)$ on a finite grid of $N$ times $t_{j}=j\Delta t$. We then approximate $\delta\rho_{nn^{\prime}}(\omega)$ by the discrete Fourier transform

[TABLE]

A common choice is to perform the rt-TDDFT calculation using a so-called $\delta$-kick, where $v(t)=K\delta(t)$ and the Fourier transform $v(\omega)=K$ is exactly constant. For all other kinds of perturbations, we approximate

[TABLE]

The finite simulation length results in a convolution of the true Fourier transform $\delta\rho_{nn^{\prime}}(\omega)$ by a sinc-shape, leading to a noisy spectrum. A common remedy is to artificially dampen the response function, forcing it to zero before the end of the simulation. In rhodent, Gaussian broadening is implemented, where the delta-kick response is multiplied by a Gaussian envelope of width $\sigma$

[TABLE]

This is equivalent to convoluting the noisy $\delta\rho_{ia}(\omega)$ with a Gaussian $e^{-\omega^{2}/2\sigma^{2}}$ in the frequency domain.

III.1.2 In the time domain through the convolution trick

The electron-hole part of the induced density matrix $\rho_{ia}$, where $f_{i}>f_{a}$ is linear in the perturbation [36], meaning that for sufficiently weak perturbations $\hat{v}$ there is a linear response regime where

[TABLE]

or equivalently, in the frequency domain,

[TABLE]

where $\hat{\chi}_{ia}$ is a response function describing the response of matrix element $\delta\rho_{ia}$ to $\hat{v}$. By convention, we use indices $i$ and $a$ to denote occupied and unoccupied ground state orbitals (equivalent to holes and electrons), respectively. Using rhodent, we can exploit this linearity to compute the response $\delta\rho_{ia}^{\prime}$ to perturbation $\hat{v}^{\prime}$ knowing the response to $\delta\rho_{ia}$ to perturbation $\hat{v}$, without performing another rt-TDDFT calculation.

[TABLE]

We are restricted to perturbations of the same spatial shape, where the onset of $v^{\prime}(t)$ is not earlier in time than $v(t)$, and where the spectrum of $v(\omega)$ covers $v^{\prime}(\omega)$ entirely. The response in the time domain can be obtained by carrying out an inverse Fourier transformation of Eq. (9).

In practice, quantities are available on a finite grid of $N$ times $t_{j}=j\Delta t$. We then take a discrete Fourier transform of the induced density matrix and both perturbations on the grid of $N$ time instances

[TABLE]

The chosen frequencies $\omega_{k}=2\pi k/(\Delta tN^{\prime})$ correspond to a grid of $N^{\prime}$ times ($N^{\prime}$ is at least $2N$), which is equivalent to padding the data with zeros after the end of the simulation. This is necessary to prevent circular correlation (see A). We obtain the induced density matrix to the new perturbation as the inverse discrete Fourier transform of Eq. (9)

[TABLE]

When constructing the density matrix from the time-dependent wave functions file, it is important to write the wave functions at every time step, if a $\delta$-kick was used. Otherwise, if the wave functions are written at a sparse interval $\Delta t$ there will be aliasing in the Fourier spectrum above the Nyquist frequency $1/2\Delta t$.

This behavior can be alleviated by using different temporal shapes. For example, a sinc-shaped pulse (Fig. 1)

[TABLE]

has the Fourier transform

[TABLE]

where the rectangular function $\text{rect}(x)$ is equal to 1 for $0<x<1$ and 0 otherwise. In principle, this pulse does not induce any response above the cutoff frequency and allows saving the time-dependent wave functions file using a sparse interval.

However, in finite simulations, $v(\omega)$ is effectively convoluted with a sinc function leading to a smooth cutoff as $v(t)$ does not fit in the simulation window in its entirety. Choosing $t_{0}$, one should strike a balance between fitting as much as possible of $v(t)$ in the simulation window (leading to a sharper cut-off) and perturbing the system with a sufficiently strong field early on in the perturbation (leading to better numerical stability). In particular, it is a judicious choice to let the offset be an integer and a quarter offset of the time between oscillations $t_{0}=(n+1/4)/(2\pi\omega_{\text{cut}})$ so that the pulse starts at a local maximum.

III.1.3 Constructing the KS density matrix

Currently, rhodent provides an interface to the linear combination of atomic orbitals (LCAO) implementation of rt-TDDFT in gpaw [27, 4, 37]. In this implementation, the KS wave functions are represented as the linear combination

[TABLE]

where $\Phi_{\mu}(\boldsymbol{r})$ are atom-centered orbitals (typically a few orbitals for each atom) and only the coefficients $C_{n\mu}$ are time-dependent.

III.1.4 Time-dependent wave functions file

In gpaw, the time-dependent coefficients $C_{n\mu}$ can be written to file periodically during time propagation. From this file, we can construct the KS density matrix in the LCAO basis

[TABLE]

Combining Eq. (1) and Eq. (17), we find the basis transformation matrix $P_{n\mu}=\sum_{\nu}C_{n\mu}(0)S_{\nu\mu}$ with the overlap integral $S_{\mu\nu}=\int\Phi^{*}_{\nu}(\boldsymbol{r})\Phi_{\mu}(\boldsymbol{r})\mathrm{d}{\boldsymbol{r}}$. Then the induced density matrix in the basis of ground state KS orbitals is

[TABLE]

III.1.5 Fourier transform of induced KS density matrix.

gpaw also supports building the Fourier transform $\delta\rho_{\mu\nu}(\omega)=\int_{0}^{\infty}[\rho_{\mu\nu}(t)-\rho_{\mu\nu}(0)]e^{i\omega t}\mathrm{d}t$ on the fly during time propagation on a grid of predefined frequencies [35]. Having saved this quantity to file, rhodent can read it and transform the density matrix to the basis of ground state KS orbitals as

[TABLE]

III.1.6 Reading and writing numpy binary files

The construction of the density matrix from different sources, and transformations between time and frequency domain, are time and memory intensive. In particular, memory usage can be limiting in parallel execution as, for example, the transformation matrix $P_{n\mu}$ is duplicated on every process. Therefore, rhodent allows writing the density matrix in the KS basis to a NumPy binary file on disk after every transformation. The files can be read in order to continue with another transformation (for example, the pulse convolution outlined in the next section), or simply to compute observables.

III.2 Calculating observables

In rhodent Calculators are used to evaluate observables, which can be computed using the full response $\delta\rho_{ia}$, including:

• •

The density in the time domain, or the Fourier transform thereof.

• •

The dipole moment in the time domain, or its Fourier transform (the polarizability in the frequency domain).

• •

Stored energy in the system in the time domain.

• •

Electron and hole occupations, i.e., hot carrier (HC) distributions, in the time domain.

Additionally, a density of states calculator, which only relies on ground state information, and a spectrum calculator, which only needs a time-dependent dipole moment file, are implemented. For analysis purposes, the observables can be decomposed either by the energies or spatial confinement of occupied and unoccupied states $i$ and $a$. For the latter, we use the Voronoi cell of an atom, which is the set of all points in space closer to that atom than to any other, to define the Voronoi weights

[TABLE]

In the following, we use Hartree atomic units.

Density

The induced charge density is obtained from $\delta\rho_{ia}$ as

[TABLE]

Because the electron-electron and hole-hole parts of the KS density matrix ($\delta\rho_{nn^{\prime}}$ where $f_{n}=f_{n^{\prime}}$) are quadratic in perturbation, the electron-hole part dominates. Using that $\rho_{ia}=\rho_{ai}^{*}$ and that the KS orbitals are real, we can then write

[TABLE]

An energy filter can be supplied to rhodent to include only a subset of transitions in the sum. For example, one could want to include only transitions to electrons above $\varepsilon_{\text{low}}$ in energy

[TABLE]

Dipole moment

The induced dipole moment (dipole moment minus the static part) is

[TABLE]

where

[TABLE]

We can use rhodent to calculate the induced dipole Eq. (24), the induced dipole projected on Voronoi weights of occupied and unoccupied states for a projection on one or several atoms

[TABLE]

or the dipole or projected dipole as a transition contribution map (TCM) [38]

[TABLE]

The TCM enables illustrative analyses of matrices in the basis of ground state KS orbitals, by broadening them onto energy axes using Gaussians

[TABLE]

For a perturbation corresponding to a spatially constant electric field polarized in the $x$ direction

[TABLE]

the normalized Fourier transform of the induced dipole moment gives one column of the polarizability

[TABLE]

In order to construct the full polarizability tensor, three rt-TDDFT calculations with orthogonal polarization directions along $x$, $y$, and $z$ are needed

[TABLE]

The optical absorption spectrum, expressed as the oscillator strength function, is related to the polarizability as

[TABLE]

This quantity obeys the sum rule that its integral $\int_{0}^{\infty}S(\omega)\mathrm{d}\omega$ equals the number of electrons in the system. The same kind of decomposition is available in the frequency domain as in the time domain. For example, we can construct the TCM [38] for the photoabsorption cross section

[TABLE]

Stored energy

In KS DFT, the total energy $E_{\text{tot}}$ can be partitioned into the kinetic energy $T$, the potential energy due to the Hartree and exchange-correlation (XC) potential $E_{\text{Hxc}}$, and the external energy. Using perturbation expansions up to second order [36] we can write the decomposition of the energy stored in the system

[TABLE]

where the first term is the total energy of the ground state and the last term is the external energy due to the electric field of the perturbation

[TABLE]

Here, we assume that the electric field of the perturbation is spatially constant and polarized in $\hat{\boldsymbol{e}}$ direction. The change in kinetic and Hartree-exchange-correlation (Hxc) energies can be decomposed into contributions from all electron-hole transitions

[TABLE]

where

[TABLE]

are calculated from the following quantities and their time derivatives

[TABLE]

We can also calculate the rate of energy change using higher derivatives

[TABLE]

For energies, we can compute Voronoi projections analogously to Eq. (26), and perform energy filtering based on the energy of pairs $ia$. For example, we can compute the total and Hxc energy above a threshold $\varepsilon_{\text{low}}$

[TABLE]

Similarly, we can compute TCMs analogously to Eq. (27) and Eq. (28).

Hot carriers

The equal-occupation-numbers part ($f_{n}=f_{n^{\prime}}$) of the induced density matrix is quadratic in the perturbation [36], and can be obtained from the electron-hole part

[TABLE]

where

[TABLE]

The purely electron-electron ($\delta\rho_{aa^{\prime}}$, $f_{a}=f_{a^{\prime}}=2$) and hole-hole parts ($\delta\rho_{ii^{\prime}}$, $f_{i}=f_{i^{\prime}}=0$) of the KS density matrix are equal to $\delta\rho^{\text{electrons}}_{aa^{\prime}}$ and $\delta\rho^{\text{holes}}_{aa^{\prime}}$, respectively, but for the partially occupied states $0<f_{n}=f_{n^{\prime}}<2$ both terms in Eq. (48) are non-zero.

Hot-carrier distributions are calculated by convolution of $\delta\rho_{ii^{\prime}}$ (holes) and $\delta\rho_{aa^{\prime}}$ (electrons) with Gaussian broadening functions on the energy grid.

[TABLE]

We can also project hot-carrier distributions according to the Voronoi weights

[TABLE]

Optical spectrum and densities of state

Optical response calculations are a routine task in rt-TDDFT since the calculations by Yabana and Bertsch [3]. For spectroscopic applications, the dipolar response of the matter system due to a dipolar external field is of interest, as it is the only non-vanishing contribution in the far field.

[TABLE]

or, equivalently, in the time domain

[TABLE]

Finally, we can compute the total density of states (DOS)

[TABLE]

and the projected density of states (PDOS) as

[TABLE]

where $w^{\text{proj}}_{kk}$ are the Voronoi weights given by Eq. (20).

IV Examples

In this section, we study three plasmonic systems to demonstrate the capabilities of rhodent. The calculation workflow is as follows: For each system, we first compute the electronic ground state within DFT using LCAO basis sets in gpaw. Then we perform an rt-TDDFT calculation using gpaw, propagating the ground state in time under an external perturbation. We simulate the dynamics for \qty30\femto, which is enough to resolve the formation and subsequent dephasing of the localized surface plasmon resonance (LSPR). A propagation time step of \qty10\atto yields converged results. We apply a sufficiently weak perturbation to only probe the linear-response domain. By choosing a sinc-shaped (Eq. (14)) perturbation in the time domain, we avoid exciting high frequency components of the response, which allows us to write a wave function trajectory to file at a relatively sparse interval, without risking aliasing. This file is read by rhodent to calculate the response. Further information concerning the DFT and rt-TDDFT calculations can be found in Sect. VI.

IV.1 Aluminum nanoparticle

As an initial example, we consider a plasmonic 201-atom Al nanoparticle (NP), with a diameter of \qty1.6\nano. An in-depth analysis of this system can be found in Ref. [10]. We perform one rt-TDDFT calculation, using a sinc-pulse with a cutoff of \qty16, and write the wave function trajectory to disk at an interval of \qty100\atto (Nyquist frequency of \qty20.7). Due to the symmetry of the NP, its response is isotropic, and one polarization direction is enough to probe the entire response.

Analysis in the frequency domain

Using the dipole moment file from the rt-TDDFT calculation, and Eq. (33), we construct a photoabsorption spectrum (Fig. 2a) for the system. The main feature in the spectrum is a peak at \qty7.6. We attribute the peak to the LSPR, which is a collective electronic excitation [39, 40]. We can gain insight about the microscopic origin of this feature by constructing the TCM to the absorption, according to Eq. (34). This calculation requires the full response of the system, i.e., the wave function trajectory. Visualizing the contributions as a TCM (Fig. 2b), we see that many individual electron-hole pairs near the Fermi level collectively contribute to the absorption. We also compute the DOS, using Eq. (57), and plot it together with the TCM. Despite the energy difference between most involved electron and hole states being less than 3 to \qty4 or so, they are excited at \qty7.6, thanks to their collective coupling through Coulomb interactions.

Analysis in the time domain

Now we consider the response of the Al NP to a Gaussian laser pulse. We consider a spatially constant electric field

[TABLE]

choosing the strength $|\boldsymbol{E}_{0}|=\qty{51}{\mu/}$, an envelope centered at $t_{0}=\qty{10}{\femto}$, and $\tau=\qty{2.19}{\femto}$, giving a full width at half-maximum (FWHM) of about \qty5.2\femto in the time domain and \qty0.7 in the frequency domain. Thanks to the pulse convolution functionality in rhodent (see Sect. III.1.2), we do not need to perform any new rt-TDDFT calculations, which is the most expensive part of the simulation workflow. By appropriately setting up the Response object, the pulse convolution is performed before calculating observables.

We excite the system at the LSPR ($\hbar\omega_{0}=\qty{7.6}{})$, and evaluate the time-dependence of the energy stored in the system (Fig. 3). As the pulse is turned on, the system responds with a dipole moment $\mu$ that lags roughly $\pi/2$ in phase after the electric field, corresponding to resonant response (absorption), and peaks in amplitude a few oscillations after the maximum of the pulse. Therefore, the energy contribution associated with the pulse (Eq. (36)) oscillates around zero at twice the pulse frequency.

The total energy increases during the entire duration of the pulse until it reaches a steady value. Having access to the full response, we can partition the total energy into kinetic and Hxc contributions, $\delta T(t)$ and $\delta E_{\text{Hxc}}(t)$, see Eq. (37) and Eq. (38). We can see how energy is periodically redistributed between the former and the latter, at a frequency matching the pulse. Increasingly more energy is stored in kinetic contributions, while the amount of Hxc energy periodically reaches zero. The Hxc contributions persist for a few \qty\femto after the disappearance of the pulse, so we interpret this as the lifetime of the LSPR. Alternatively, we can partition the total energy into resonant and non-resonant contributions, by the distinction that the eigenvalue difference of resonant electron-hole pairs should be $|\hbar\omega_{ia}-\hbar\omega_{0}|<\qty{1.4}{}$. This tells us that the number of non-resonant transitions rises and falls smoothly, with only small oscillations, during the lifetime of the LSPR.

For even deeper insight into the data, we can construct TCMs for the energy contributions at selected times (Fig. 4). At the time instance marked (1) in Fig. 4a, when the Hxc contributions are at a maximum, we see both resonant and non-resonant transitions in the total energy (Fig. 4b), and only non-resonant transitions in the Hxc contribution (Fig. 4c). After one pulse cycle, at time instance (2), there are no Hxc contributions, but the same non-resonant transitions persist in the total energy (Fig. 4d). We should note here that these are the same transitions that are visible in the absorption TCM (Fig. 2b). It is now clear that these transitions make up the LSPR, and that they oscillate between kinetic and potential energy during its lifetime. After the decay of the pulse, the plasmon dephases into resonant transitions, known as HCs. At time instance (3) only such transitions remain (Fig. 4e).

IV.2 Aluminum nanoparticle with benzene molecules

The next example system is another structure from Ref. [10], made up of an Al NP identical to the previous structure, and two benzene molecules. The molecules are placed on opposite sides at a distance of \qty3 from the faces of the NP. The molecules are oriented such that the line through the molecules and the NP is parallel with the transition dipoles of the lowest bright excitations of the molecules. The system is no longer isotropic, so we constrain our analysis to excitation polarized along this axis. We perform a new rt-TDDFT calculation, with the same parameters as for the bare Al NP.

Analysis in the frequency domain

We calculate the photoabsorption spectrum (Fig. 5a), and find two peaks. One peak at \qty7.6 is very similar to the LSPR peak of the bare Al NP, but the other at \qty6.9 is both much larger in prominence and redshifted by about \qty0.2 compared to the peak of the bare molecules. This is a characteristic of strong coupling, where the underlying excitations (LSPR of the NP and the molecular excitation) exchange energy with each other at a faster or comparable rate to the lifetime of the excitations [41]. The two resonances are known as the LP and UP.

We decompose the absorption TCMs into NP-like and molecular transitions, using the Voronoi decomposition Eq. (26). In particular, we attribute transitions from occupied NP-states into any unoccupied state as the former, and transitions from occupied molecular states to any unoccupied state as the latter, because unoccupied states tend to be more of mixed character compared to occupied states. In Fig. 5b–e, we plot the projected TCMs together with the DOS (Eq. (57)) on the unoccupied energy axes, and the PDOS (Eq. (58)) for the NP/molecule on the occupied energy axes. At the UP frequency $\hbar\omega_{0}=\qty{7.6}{}$, we find that the NP decomposition of the TCM (Fig. 5b) consists of the same excitations near the Fermi level that we attributed to the LSPR in the bare Al NP. Meanwhile, the molecular decomposition of the TCM (Fig. 5c) consists of a single excitation from the benzene highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO). The molecular transition is contributing destructively to the absorption at the frequency of the UP, meaning that the molecular dipole is directed opposite to the NP dipole. At the LP frequency $\hbar\omega_{0}=\qty{6.9}{}$, the TCM also consists of plasmonic transitions in the NP decomposition (Fig. 5d) and the single molecular excitation in the molecular decomposition (Fig. 5e). The difference is that for the LP, the LSPR and molecular excitation couple constructively. A configuration where the dipoles are parallel intuitively leads to a lower energy. The LP has significant molecular character, as \qty14.3% of the absorption comes from the molecular transition, while the UP is almost entirely LSPR-like (the destructive contribution of the molecule is only \qty2.0%).

Analysis in the time domain

Now we consider the time evolution of the dipole moment and energies after excitation by a Gaussian laser pulse. We excite the system with a pulse covering both peaks, with parameters $\hbar\omega_{0}=\qty{7.25}{}$, $t_{0}=\qty{5}{\femto}$, and $\tau=\qty{1.10}{\femto}$ (giving a FWHM of about \qty2.6\femto in the time domain and \qty1.4 in the frequency domain). We decompose the energies and dipoles in the same way as for the earlier frequency domain analysis, and find that a strong oscillating dipole first forms in the NP, then in the molecules, and then in the NP again (Fig. 6). These are Rabi oscillations, with 2-\qty3\femto in between the maxima. Oscillations in the Hxc energy are visible in the NP and molecules during the same times as oscillations in the dipole. After \qty20\femto of simulation time, only kinetic energy remains (in the form of HCs), with a 92.3/\qty7.7% division between the NP and molecules.

IV.3 Silver nanoparticle with CO molecule

As a final example, we consider a metal with valence d-electrons. We set up a 201-atom Ag NP (nearly the same diameter and shape as the Al NP), and place a CO molecule at a distance of \qty3 from the (111) on-top site. Detailed analysis of this system can be found in Ref. [23] We perform yet another rt-TDDFT calculation, with a polarization direction along the bond axis of the molecule. As for the Al systems, we use a time step of \qty10\atto and a total simulation length of \qty30\femto. Because the relevant resonances are lower for this system compared to the Al NP, we use a sinc-pulse with a cutoff of \qty8. This allows us to write the wave function trajectory to disk at an interval of \qty200\atto (Nyquist frequency of \qty10.3).

Analysis in the frequency domain

We calculate the photoabsorption spectrum (Fig. 7a) and identify the LSPR peak at \qty3.8. Unlike the Al–benzene system we studied previously, the first excitation of the CO molecule is far from the LSPR and its presence has a negligible effect on the absorption spectrum. We construct a TCM for the absorption spectrum (Fig. 7b) and see many excitations near the Fermi level constructively coupling at the LSPR frequency. A new feature, not found in the Al systems, are the destructive contributions involving transitions from the d-band edge (about \qty4 below the Fermi level) to unoccupied levels. These transitions are excited at the same time as the LSPR, with their dipoles opposite. Thus they screen the plasmon, contributing to a much lower resonance energy (\qty3.8) compared to the Al NP (\qty7.6).

Analysis in the time domain

Once again, we calculate the response to a Gaussian laser pulse. We choose the parameters $\hbar\omega_{0}=\qty{3.8}{}$, $t_{0}=\qty{10}{\femto}$, and $\tau=\qty{2.19}{\femto}$ (FWHM \qty7.3\femto in the time domain and \qty0.7 in the frequency domain) in order to cover the LSPR peak. For this system, we focus on the time after the dephasing of the plasmon. We compute the average of the hot hole (HH) and hot electron (HE) distributions (Fig. 8a) during the last \qty5\femto of the simulation, using Eq. (51) and Eq. (52). We see that HCs are created in the NP, with energies spanning from 0 to the pulse frequency (or negative pulse frequency, for holes). As we saw in the energy TCM for the Al-system (Fig. 4e), HCs present after the dephasing of the plasmon are found on the diagonal $\varepsilon_{a}-\varepsilon_{i}=\hbar\omega_{0}$, where $\varepsilon_{a}$ and $\varepsilon_{i}$ are the energies of the electron and hole respectively. We also see this in that the HE distribution looks like the HH distribution, shifted upwards in energy by the pulse energy, except for some minor shifts due to the finite width of the pulse. At the Fermi level, there are both excited electrons and holes, due to the presence of partially occupied states in the ground state calculation.

For the Ag NP + CO molecule system, the molecular levels are hybridized with the metal states. The LUMO is split into several levels, around 1 to \qty2 above the Fermi level (Fig. 8b; top axes), while the HOMO is too far below the Fermi level to allow for excitations at the present pulse energy. Therefore, we can expect some HEs to form in the molecule but not any HHs. Indeed, computing the HC distribution projected on the molecule (Eq. (53) and Eq. (54)) confirms that electrons are excited in the molecule in all branches of the hybridized LUMO (Fig. 8b). Interestingly, there is a delicate interplay between the coupling of each possible HC to the LSPR and the energetic alignment to the pulse that determines in which hybridized level most electrons are excited. This has been exploited in [23] to increase the amount of charge transfer.

Using rhodent, we can, without much additional effort, compute the HC distributions for many different pulses at once. In Fig. 8c, we plot the total number of excited electrons in the system, as well as the number of excited electrons in the molecule, by pulse frequency. The total number of electrons is related to the amount of energy absorbed, which can be deduced from the absorption spectrum pulse frequency roughly like the absorption spectrum, as all the energy absorbed turns into HCs. For the number of electrons in the molecule, there are, however, several local maxima.

V Conclusions and outlook

rhodent is a Python package for post-processing the response of rt-TDDFT calculations. It can be used to obtain, e.g., induced dipole moments and densities, the stored energy and HC distributions, as well as spatial and energetic projections thereof. Induced dipole moments and densities can also be computed in the frequency domain (the former being related to optical spectra). Additionally, since rt-TDDFT probes the linear response function when the external perturbation is sufficiently weak, rhodent can quickly compute the response to any time-dependent perturbation in the linear response regime, from just one rt-TDDFT calculation. In the linear response regime, we can also exploit the structure of the density matrix to compute energies and HC distributions in the system.

Currently, rhodent supports reading output files from calculations done using gpaw. However, thanks to the modular structure of the rhodent code, with a clear separation between the construction of response and the computation of observables, support for other rt-TDDFT softwares can be readily added. On a similar note, only calculations done in the LCAO basis are currently supported. For LCAO calculations, we can transform the KS density matrix to the basis of ground state KS orbitals without any loss of information, and still have a quantity of tractable size. For, e.g., uniform grid or plane-wave based calculations, the number of basis functions is much larger, and the basis of KS orbitals would have to be truncated.

VI Computational details

The open-source gpaw [27] code package was used for the DFT and rt-TDDFT calculations. The calculations were done in the projector augmented wave (PAW) [42] formalism using LCAO basis sets [37]; the pvalence [4] basis set, which is optimized to represent bound unoccupied states, was used for the Ag species, and the double-zeta polarized (dzp) basis set for all other species. The PBE [43, 44] (for the Al and Al + benzene calculations) and Gritsenko-van Leeuwen-van Lenthe-Baerends-solid-correlation (GLLB-sc) [45, 46] (for the Ag + CO calculations) XC-functionals, utilizing the Libxc [47] library, were used in gpaw. For the Al and Al + benzene calculations, a simulation cell of \qtyproduct28.8 x 28.8 x 43.2 was used. For the Ag + CO calculations, a simulation cell of \qtyproduct32.0 x 32.0 x 35.2 was used. In the simulation cell, wave functions were represented with a grid spacing of \qty0.2, and XC and Coulomb potentials with a grid spacing of \qty0.1. Additional analytic moment corrections [48] centered at the NP were added to the Coulomb potential. Fermi-Dirac occupation number smearing with width \qty0.05 was used. The self-consistent loop was stopped when the integral of the difference between two subsequent densities was less than \qty1e-12. Pulay [49]-mixing was used to accelerate the ground state convergence.

Data availability

The data generated for the examples in this article is available on Zenodo at https://doi.org/10.5281/zenodo.16746428. The rhodent releases are accessible on Zenodo at https://doi.org/10.5281/zenodo.13332634 while the documentation can be accessed at https://rhodent.materialsmodeling.org.

Acknowledgments

We acknowledge funding from the Knut and Alice Wallenberg foundation (Grant No. 2019.0140; J. F. and P. E.) and the Swedish Research Council (No. 2020-04935; J. F. and P. E.). The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at NSC, PDC and C3SE partially funded by the Swedish Research Council through grant agreement no. 2022-06725.

Appendix A The convolution trick in finite simulations

In the linear response regime, we have that (dropping indices $ia$ from the notation in the main text)

[TABLE]

where we assume that the perturbation $v(t)$ is zero before time zero, and the response is causal $\chi(t)=0$ for $t<0$. We can then work in the frequency domain

[TABLE]

where the Fourier transform is defined

[TABLE]

and likewise for $\chi(\omega)$ and $v(\omega)$. As the Fourier transform is defined as an integration with infinite limits, and our simulations are finite, we need to reformulate the above relations in terms of finite integrations. Here, we derive a frequency domain formula for quantities sampled in the time window $0<t<T$.

First, we define the zero-padded quantity

[TABLE]

which we take to be periodic with the period $T^{\prime}$. Next, the crucial assumption is that $T^{\prime}$ is at least twice the length of $T$. Then, for $0<t<T$, we can swap the perturbation by the zero-padded perturbation, and extend the limits of the convolution Eq. (60)

[TABLE]

This holds because the added integral is zero

[TABLE]

We can now expand the zero-padded perturbation in a Fourier series of periodicity $T^{\prime}$

[TABLE]

where $\omega_{k}=2\pi k/T^{\prime}$ and the Fourier coefficients

[TABLE]

only require integration over the true perturbation in the window $0<t<T$. Inserting the expansion of $v^{(0)}$ into Eq. (64), we get

[TABLE]

where we have defined the Fourier coefficients of the zero-padded response function

[TABLE]

We identify Eq. (69) as a Fourier series for $\delta\rho(t)$, valid in the simulation time window $0<t<T$, with the Fourier coefficients

[TABLE]

These coefficients are not, in general, equal to the coefficients of the zero-padded induced density matrix

[TABLE]

but the Fourier series using $\delta\hat{\rho}_{k}$ and $\delta\hat{\rho}^{(0)}_{k}$ are equal inside the simulation window. The Fourier series Eq. (69) is not, in general, zero outside the simulation time window $0<t<T$.

Now, we consider a different pulse $v^{\prime}(t)$ leading to a different response (new quantities denoted by primes)

[TABLE]

with Fourier coefficients defined

[TABLE]

and, as before, the Fourier coefficients $\delta\hat{\rho}^{\prime}_{k}=\hat{\chi}^{(0)}_{k}\hat{v}^{\prime(0)}_{k}$ reproduce $\delta\rho^{\prime}(t)$ in the simulation time window $0<t<T$.

Now we require that the Fourier coefficients of the new pulse are non-zero for every non-zero coefficient of the old pulse, so that we can write

[TABLE]

which is equivalent (according to the circular correlation theorem) to the pulses being related by the convolution

[TABLE]

where $K^{(0)}(t)$ is some zero-padded function and $K^{(0)}_{k}$ its Fourier coefficients. This means that the onset of the new pulse must be no later in time than the onset of the old pulse and introduce no new frequencies.

Then the Fourier series of the new response can be written

[TABLE]

which is also equivalent to a convolution with $K^{(0)}$. Because we are only interested in the response $\delta\rho^{\prime}(t)$ during the simulation time $0<t<T$, and $K^{(0)}$ is zero-padded, we can thus swap $\delta{\hat{\rho}}_{k}$ for $\delta{\hat{\rho}}^{(0)}_{k}$ and calculate the new response as the Fourier series

[TABLE]

We have used that $K_{k}=\hat{v}^{\prime(0)}_{k}/\hat{v}^{(0)}_{k}$ for every non-zero value of $\hat{v}^{\prime(0)}_{k}$.

In practice, we have sampled $\delta\rho(t)$ and $v(t)$ on a grid of $N$ times $t_{j}=j\Delta t$ with time step $\Delta t=T/N^{\prime}$. Then the Fourier coefficients are approximated as

[TABLE]

and

[TABLE]

where the terms inside the brackets are the discrete Fourier transform and inverse discrete Fourier transform respectively, and $N^{\prime}\geq 2N$. The approximation is good if the perturbation does not have any response above the Nyquist frequency $\omega=2\pi/(2\Delta t)$; otherwise, aliasing effects are seen.