<figure><img src="image_22.png" /><figcaption>FIG. 10. Calculated components (a) $\epsilon_{xx}$ and (b) $\epsilon_{xy}$ of the dielectric tensor of silicon for different frequencies of light $\Omega_0$ (in eV) and the magnetic field of 1 T along the z axis (with and without account of excitonic effects). The linewidth $\delta = 0.1$ eV is used. The blue shift by 0.7 eV is applied to take into account the GW correction to the bandgap.<sup>130,131</sup> The experimental data for $\epsilon_{xx}$<sup>132</sup> and $\epsilon_{xy}$<sup>110</sup> (scaled by 1/2) are indicated by symbols.</figcaption></figure>

accounting for excitonic effects show a better agreement with the experimental results [Figs. 10(a) and 10(b)]. Although the magnitudes of the peaks in the magneto-optical spectra are about a factor of two smaller than in the magneto-optical measurements,<sup>110</sup> they can be corrected by reducing the linewidth $\delta$ assumed in the calculations.

In the limit of a large supercell, this formalism becomes equivalent to the simpler standard formulation for finite systems, which we have also implemented for reference. In this case, the Liouville equation for the density matrix is written in the Coulomb gauge as

$$ \Omega\rho + [H_0, \rho] = [\mathbf{d} \cdot \mathbf{E} + \mathbf{m} \cdot \mathbf{B} - \delta H, \rho], \quad (43) $$

where $\mathbf{d} = -\mathbf{r}$ is the electric dipole moment, $\mathbf{m} = -\mathbf{r} \times \mathbf{V}/2c$ is the orbital magnetic dipole moment, and $\delta H$ describes local-field effects. The change in the polarizability $\alpha_{v\mu,\gamma}B_\gamma$ in the presence of the magnetic field is calculated from the dipole response

$$ \alpha_{v\mu,\gamma}(\Omega) = \mathrm{Tr}[d_v\rho^{(E_\mu B_\gamma)}(\Omega)]. \quad (44) $$

As in the periodic case, due to the “$2n + 1$” theorem,<sup>127,128</sup> the explicit calculation of the second-order derivative $\rho^{(E_\mu B_\gamma)}$ is avoided. Instead, a supplementary electric field at frequency $-\Omega$, $E_v^\prime$, is introduced.<sup>116</sup> In the case of real wavefunctions, $|\eta_v^{(E_\mu)}(-\Omega)\rangle = (|\eta_v^{(E_\mu)}(-\Omega^*)\rangle)^*$, and there is no need to solve additionally the Liouville equation for this supplementary electric field.

The efficiency of the present scheme to compute magneto-optics is comparable to standard linear-response calculations of simple optical polarizability in the absence of the magnetic field.<sup>116</sup> When local-field effects are included self-consistently, the calculations of magneto-optical spectra for solids take only twice as long as those of polarizability. For finite systems, the computational effort is the same as for the simple optical polarizability.

## X. TIME-DEPENDENT ANGULAR RESOLVED PHOTOELECTRON SPECTROSCOPY

The real-space representation of the dynamics of the electronic structure allows for a seamless and straightforward description of dynamical processes outside the material, i.e., processes where electrons are excited into vacuum. Beyond simply describing the ionization process, Octopus has routines implemented that compute photoelectron spectroscopy in different flavors. Photoelectron spectroscopy is particularly ubiquitous for the characterization of the electronic structure in solids because it provides a direct observable of the energy and momentum distribution of electronic states, known as the bandstructure. In contrast to other electronic structure methods that compute quantities linked to the photoelectron spectrum, most commonly the single-quasiparticle spectral function obtained through the GW approximation to the many-body self-energy, the approach described here does not consider unit cells of bulk materials but instead computes directly the energy and momentum resolved ionization probability.<sup>107,133</sup> Besides the study of the electronic structure of solids, photoelectron spectra of atoms and molecules are of equal interest and the implementation in Octopus is capable to describe accurately all such systems on an equal footing.

The most detailed quantity available in the experiments is the momentum-resolved photoelectron probability $P(\mathbf{p})$, i.e., the probability to detect an electron with a given momentum $\mathbf{p}$. In some cases, the experimental setup offers the possibility to measure the spin polarization along a given axis. The formalism we are going to outline can easily accommodate a non-collinear spin structure and, therefore, calculate spin-resolved quantities, but for the sake of simplicity, in the following, we are going to restrict ourselves to closed shell systems with collinear spins, i.e., where all the orbitals are doubly occupied with electrons having opposite spins. The reader interested in the most general case can find a detailed description in Ref. 134. From $P(\mathbf{p})$, one can obtain other derived quantities by simple manipulation. For instance, the energy-resolved spectrum $P(E)$, used to identify the occupied energy levels, can be obtained with a change of variable using the free electron energy dispersion relation $E = \rho^2/2$. The angle-resolved photoelectron spectrum (ARPES), normally employed to measure the quasiparticle bandstructure, is simply obtained by taking the energy resolved spectrum as a function of the electron momentum parallel to the surface $P(\mathbf{p}_{\parallel}, E)$.

The t-SURFF method was first proposed by Tao and Scrinzi<sup>135</sup> for one-electron systems and later extended to many electrons with TDDFT for periodic<sup>134</sup> and non-periodic systems.<sup>136</sup> It is based on the assumption that the Kohn-Sham Hamiltonian describing the full experimental process, i.e., including the ionization and detection, can be decomposed into the sum of two Hamiltonians acting into complementary spatial regions, *inner* and *outer*, and that can be approximated in different ways. In the *inner* region surrounding the system, the electron dynamics is governed by the interacting Kohn-Sham Hamiltonian $\hat{H}_{\mathrm{KS}}(\mathbf{r}, t)$, while in the *outer* region, electrons are free from the Coulomb tails of the parent system and behave as independent particles driven by an external field and, therefore, are described by the Volkov Hamiltonian $\hat{H}_{\mathrm{V}}(\mathbf{r}, t) = 1/2(-i\nabla - \mathbf{A}(t)/c)^2$. We express the field with a time dependent vector potential $\mathbf{A}(t)$ in the dipole approximation, i.e., by discarding the spatial dependence of the field $\mathbf{A}(\mathbf{r}, t) \approx \mathbf{A}(t)$. The
advantage of this approach is provided by the fact that the time-dependent Schrödinger equation associated with the Volkov Hamiltonian can be solved analytically and that the solutions, the Volkov waves,

$$ \chi_{\mathbf{p}}(\mathbf{r}, t) = \frac{1}{(2\pi)^{\frac{3}{2}}} e^{i\mathbf{p}\cdot\mathbf{r}} e^{i\Phi(\mathbf{p},t)} \quad (45) $$

with

$$ \Phi(\mathbf{p}, t) = \int_{0}^{t} d\tau \left( \mathbf{p} - \frac{\mathbf{A}(t)}{c} \right)^{2}, \quad (46) $$

are eigenstates of the momentum operator. We can, therefore, expand the Kohn–Sham orbitals $\phi_i$ as a superposition of detector states in the form of Volkov waves

$$ \phi_i(\mathbf{r}, t) = \int d\mathbf{p} b_i(\mathbf{p}, t) \chi_{\mathbf{p}}(\mathbf{r}, t) \quad (47) $$

and obtain the photoelectron probability in terms of the expansion coefficients by summing up the contribution of all the orbitals: $\mathcal{P}(\mathbf{p}) = \lim_{t \to \infty} 2/N \sum_{i=1}^{N/2} |b_i(\mathbf{p}, t)|^2$.

Using the continuity equation, t-SURFF allows us to express the coefficients $b_i$, and thus $\mathcal{P}(\mathbf{p})$, as a time integral of the photo-current flux through the surface S separating the domains of the two Hamiltonians. More specifically,

$$ b_i(\mathbf{p}, t) = - \int_0^t d\tau \oint_S ds \cdot \langle \chi_p(\tau) | \hat{\mathbf{j}} | \phi_i(\tau) \rangle \quad (48) $$

with the single-particle current density operator matrix element given by

$$ \langle \chi_{\mathbf{p}}(\tau) | \hat{\mathbf{j}} | \phi_i(\tau) \rangle = \frac{1}{2} \left\{ i\phi_i(\mathbf{r}, \tau) \nabla \chi_{\mathbf{p}}^*(\mathbf{r}, \tau) - i\chi_{\mathbf{p}}^*(\mathbf{r}, \tau) \nabla \phi_i(\mathbf{r}, \tau) - 2 \frac{\mathbf{A}(t)}{c} \chi_{\mathbf{p}}^*(\mathbf{r}, \tau) \phi_i(\mathbf{r}, \tau) \right\}. \quad (49) $$

The flexibility offered by the definition of the boundary surface S allows us to easily adapt t-SURFF to periodic and non-periodic systems. For periodic systems, we choose S as a plane parallel to the material surface, while for the non-periodic case, the most natural choice is a sphere, as shown in Figs. 11(a) and 11(b).

For periodic systems, one can apply the Bloch theorem to the Volkov Hamiltonian and further reduce Eq. (48) to an expression where only the periodic part of the Bloch waves is employed.¹³⁴

The present computation of the photoelectron spectrum relies on the condition that the Volkov Hamiltonian can be used to describe electrons in a region of space and entails the following approximations: (i) the ionized electrons become non-interacting at some distance from the system and (ii) they are free with respect to the tail of the Coulomb potential of the system, i.e., it has to be possible to neglect it at some distance. Both conditions are controlled by the placement of the analyzing surface S for the evaluation of the flux, and in practice, this has to be converged by placing it at

<figure><img src="image_21.png" /><figcaption>FIG. 11. Calculated photoelectron spectrum with t-SURFF, as implemented in Octopus.¹³⁶ The geometries employed to calculate the flux of the photoelectron current are depicted in (a) for periodic and (b) for non-periodic systems. (c) Pump-probe ARPES spectrum of monolayer h-BN driven by a laser field resonant with the gap at K, adapted from Ref. 134. (d) Photoelectron angular distribution obtained by strong-field ionization of $C_{60}$ with an IR field polarized along x, adapted from Ref. 136. All the electrons rescattering at the same time in one period of the field end up with final momenta forming a ring (dashed) centered at the value of $\mathbf{A}(t_r)$ at the rescattering time, $t_r$. The two arrows represent the graphical decomposition of the final momentum of the photoelectron in the vector potential at the moment of rescattering (horizontal arrow) plus the rescattering momentum.</figcaption></figure>
advantage of this approach is provided by the fact that the time-dependent Schrödinger equation associated with the Volkov Hamiltonian can be solved analytically and that the solutions, the Volkov waves,

$$ \chi_{\mathbf{p}}(\mathbf{r}, t) = \frac{1}{(2\pi)^{\frac{3}{2}}} e^{i\mathbf{p}\cdot\mathbf{r}} e^{i\Phi(\mathbf{p},t)} \quad (45) $$

with

$$ \Phi(\mathbf{p}, t) = \int_{0}^{t} d\tau \left( \mathbf{p} - \frac{\mathbf{A}(t)}{c} \right)^{2}, \quad (46) $$

are eigenstates of the momentum operator. We can, therefore, expand the Kohn–Sham orbitals $\phi_i$ as a superposition of detector states in the form of Volkov waves

$$ \phi_i(\mathbf{r}, t) = \int d\mathbf{p} b_i(\mathbf{p}, t) \chi_{\mathbf{p}}(\mathbf{r}, t) \quad (47) $$

and obtain the photoelectron probability in terms of the expansion coefficients by summing up the contribution of all the orbitals: $\mathcal{P}(\mathbf{p}) = \lim_{t \to \infty} 2/N \sum_{i=1}^{N/2} |b_i(\mathbf{p}, t)|^2$.

Using the continuity equation, t-SURFF allows us to express the coefficients $b_i$, and thus $\mathcal{P}(\mathbf{p})$, as a time integral of the photo-current flux through the surface S separating the domains of the two Hamiltonians. More specifically,

$$ b_i(\mathbf{p}, t) = - \int_0^t d\tau \oint_S ds \cdot \langle \chi_{\mathbf{p}}(\tau) | \hat{\mathbf{j}} | \phi_i(\tau) \rangle \quad (48) $$

with the single-particle current density operator matrix element given by

$$ \langle \chi_{\mathbf{p}}(\tau) | \hat{\mathbf{j}} | \phi_i(\tau) \rangle = \frac{1}{2} \left\{ i\phi_i(\mathbf{r}, \tau) \nabla \chi_{\mathbf{p}}^*(\mathbf{r}, \tau) - i\chi_{\mathbf{p}}^*(\mathbf{r}, \tau) \nabla \phi_i(\mathbf{r}, \tau) - 2 \frac{\mathbf{A}(t)}{c} \chi_{\mathbf{p}}^*(\mathbf{r}, \tau) \phi_i(\mathbf{r}, \tau) \right\}. \quad (49) $$

The flexibility offered by the definition of the boundary surface S allows us to easily adapt t-SURFF to periodic and non-periodic systems. For periodic systems, we choose S as a plane parallel to the material surface, while for the non-periodic case, the most natural choice is a sphere, as shown in Figs. 11(a) and 11(b).

For periodic systems, one can apply the Bloch theorem to the Volkov Hamiltonian and further reduce Eq. (48) to an expression where only the periodic part of the Bloch waves is employed.¹³⁴

The present computation of the photoelectron spectrum relies on the condition that the Volkov Hamiltonian can be used to describe electrons in a region of space and entails the following approximations: (i) the ionized electrons become non-interacting at some distance from the system and (ii) they are free with respect to the tail of the Coulomb potential of the system, i.e., it has to be possible to neglect it at some distance. Both conditions are controlled by the placement of the analyzing surface S for the evaluation of the flux, and in practice, this has to be converged by placing it at

<figure><img src="image_23.png" /><figcaption>FIG. 11. Calculated photoelectron spectrum with t-SURFF, as implemented in Octopus.¹³⁶ The geometries employed to calculate the flux of the photoelectron current are depicted in (a) for periodic and (b) for non-periodic systems. (c) Pump-probe ARPES spectrum of monolayer h-BN driven by a laser field resonant with the gap at K, adapted from Ref. 134. (d) Photoelectron angular distribution obtained by strong-field ionization of C₆₀ with an IR field polarized along x, adapted from Ref. 136. All the electrons rescattering at the same time in one period of the field end up with final momenta forming a ring (dashed) centered at the value of $\mathbf{A}(t_r)$ at the rescattering time, $t_r$. The two arrows represent the graphical decomposition of the final momentum of the photoelectron in the vector potential at the moment of rescattering (horizontal arrow) plus the rescattering momentum.</figcaption></figure>
successively larger distances from the system. This method needs an explicit description of the vacuum around the probed system, which for solids requires an explicit construction of the electronic structure of the surface layer, and if bulk properties are probed, one needs to converge a slab of material. Compared to unit cell calculations of solids, this puts the present method at a disadvantage. However, when aiming at simulating specific experiments, it natively includes signals coming from the surface layers and can thus capture more of the experimental reality of the spectroscopic process compared to standard perturbative many-body approaches that are directly performed in the bulk unit cell.

In practical calculations, in order to avoid spurious reflection of electrons from the boundaries of the simulation box, one has to employ absorbing boundaries, and depending on the condition of their transparency, this can add up to the size of the simulation box.¹³⁷ Other than the increased size of the computational box, this method is not computationally costly, as it requires the evaluation of the gradient operator on surface points S only and straightforwardly supports the existing k-point, state, and mesh parallelizations.

Pump-probe configurations can be simulated naturally since there is no restriction on the functional form of the vector potential A(t), and therefore, we can accommodate any linear combination of pulses. As an example, in Fig. 11(c), we show the result of a simulated pump-probe time-resolved ARPES measurement of monolayer hexagonal boron nitride (hBN) driven by a field resonant with the gap at the K-point in the Brillouin zone. As it is apparent from the simulation, the excited state population transfer to the conduction band is well observed in the resulting ARPES spectrum.

The method is not limited to simple photoelectron spectroscopy but can be used to simulate complex experimental techniques, such as the reconstruction of attosecond beating by interference of two-photon transitions (RABBITT).¹³⁸ Furthermore, due to the versatility of the real-time description in Octopus, the underlying excitation does not need to be originating from an optical pulse (i.e., an external vector potential), and since the ions can move according to TDDFT-Ehrenfest dynamics,¹³⁹ one can also simulate features coming from the ionic or lattice motion, specifically electron-phonon coupling signature in ARPES.¹⁴⁰ Pump-probe simulations are not limited to the dynamical process of excitations but can also be employed to study steady-state modifications of driven electronic states.¹⁴¹ This opens the possibility of studying Floquet physics from an ab initio perspective and to simulate directly the effect that a periodic force would have on the dressed electronic structure—an important aspect to underline, given the growing interest in the field of Floquet engineering¹⁴² and Floquet analysis.¹⁴³

Finally, t-SURFF is particularly suited to simulating strong field ionization processes such as in laser-induced photoelectron diffraction (LIED) experiments, where the ionization takes place by direct tunneling into the continuum and the field is so strong as to drive electrons in trajectories recolliding with the parent system. In fact, since in Eq. (48), we accumulate the flux through S over the time of the propagation, with t-SURFF, electrons can seamlessly flow back and forth through the surface without producing any artifact in the final spectrum. As an example, in Fig. 11(d), we present the photoelectron angular distribution of C₆₀ ionized by a strong IR pulse capable of inducing rescattering dynamics, which then gets imprinted in the photoelectron spectrum as characteristic rings centered in the value of the vector potential at the instant of

rescattering.¹⁴⁴ In this regime, the result of simulations obtained with Octopus are in excellent agreement with the experiments.¹⁴⁵,¹⁴⁶

## XI. ELECTRIC AND THERMAL CONDUCTIVITIES

Electrical conductivity in real materials can be described to a first approximation by Ohm's law, which is a linear relationship between the applied electric field and the current density generated in response. The standard method to study electrical and thermal conductivity is to carry out a DFT calculation at equilibrium and apply Kubo-Greenwood linear-response theory to evaluate the conductivities. Although this method has been used successfully in the past to calculate the transport properties,¹⁴⁷,¹⁴⁸ the application is limited to systems in the linear response regime. It is known that in many cases, especially in the presence of strong external fields, Ohm's law is no longer valid and the system exhibits non-linear behavior. In order to capture this complex behavior in materials, one must go beyond simple, linear approximations and describe electron interactions in a non-trivial manner.

One promising route to go beyond linear response is to use density functional methods to study directly thermoelectric transport. This topic has received less attention but has recently been implemented in Octopus and applied to liquid aluminum.¹⁴⁹ For detailed descriptions and derivations of non-equilibrium thermoelectric phenomena using density based methods, we refer the reader to Ref. 150.

In Octopus, we are able to calculate the current density and heat current density at each step during a time-dependent simulation. A current is induced in this study by applying an electric field of the form E(t) = E₀δ(t), where E₀ describes the magnitude and direction of the field. The electric field is induced through a time-dependent vector potential [E(t) = -c⁻¹∂ᵢA(t), where c is the speed of light] in the Hamiltonian. This vector potential satisfies the periodic boundary conditions of an extended system. It should be noted that there is no limitation on the form of the applied electric field in general. We are able to evaluate the macroscopic current density as

$$J(t) = -\frac{i}{\Omega} \int dr \sum_{j}^{N} \varphi^{\dagger}(\mathbf{r}, t)[\hat{H}(t), \hat{r}]\varphi_{i}(\mathbf{r}, t), \quad (50)$$

where Ω is the volume of the unit cell. The energy current density, $\hat{J}_n(\mathbf{r}, t)$, is expressed as

$$\hat{J}_n(\mathbf{r}, t) = \hat{J}_t(\mathbf{r}, t) + \hat{J}_v(\mathbf{r}, t) + \hat{J}_u(\mathbf{r}, t) + \hat{J}_f(\mathbf{r}, t), \quad (51)$$

where $\hat{J}_t(\mathbf{r})$ is the kinetic energy contribution given as

$$[\hat{J}_t(\mathbf{r}, t)]_i = \frac{i}{8} \left( [\partial_i \hat{\varphi}^\dagger] [\nabla^2 \hat{\varphi}] - [\nabla^2 \hat{\varphi}^\dagger] [\partial_i \hat{\varphi}] - [\partial_i \nabla \hat{\varphi}^\dagger] \cdot [\nabla \hat{\varphi}] - [\nabla \hat{\varphi}^\dagger] \cdot [\partial_i \nabla \hat{\varphi}] \right) \quad (52)$$

and $\varphi = \varphi(\mathbf{r}, t)$ are understood to be the time-dependent states. The output of the time-dependent simulation can be further analyzed to yield the frequency-dependent conductivity. The electrical (or thermal) conductivity $\sigma$ can be found by Fourier transforming the corresponding current as

$$\sigma_{ij}(\omega) = \frac{1}{[\mathbf{E}_0]_i} \int_0^\infty dt e^{-i\omega t} \mathbf{J}_j(t). \quad (53)$$

In general, it is possible to study the conductivity of an extended system in the presence of an applied electric field. Here, we
successively larger distances from the system. This method needs an explicit description of the vacuum around the probed system, which for solids requires an explicit construction of the electronic structure of the surface layer, and if bulk properties are probed, one needs to converge a slab of material. Compared to unit cell calculations of solids, this puts the present method at a disadvantage. However, when aiming at simulating specific experiments, it natively includes signals coming from the surface layers and can thus capture more of the experimental reality of the spectroscopic process compared to standard perturbative many-body approaches that are directly performed in the bulk unit cell.

In practical calculations, in order to avoid spurious reflection of electrons from the boundaries of the simulation box, one has to employ absorbing boundaries, and depending on the condition of their transparency, this can add up to the size of the simulation box.¹³⁷ Other than the increased size of the computational box, this method is not computationally costly, as it requires the evaluation of the gradient operator on surface points S only and straightforwardly supports the existing k-point, state, and mesh parallelizations.

Pump-probe configurations can be simulated naturally since there is no restriction on the functional form of the vector potential A(t), and therefore, we can accommodate any linear combination of pulses. As an example, in Fig. 11(c), we show the result of a simulated pump-probe time-resolved ARPES measurement of monolayer hexagonal boron nitride (hBN) driven by a field resonant with the gap at the K-point in the Brillouin zone. As it is apparent from the simulation, the excited state population transfer to the conduction band is well observed in the resulting ARPES spectrum.

The method is not limited to simple photoelectron spectroscopy but can be used to simulate complex experimental techniques, such as the reconstruction of attosecond beating by interference of two-photon transitions (RABBITT).¹³⁸ Furthermore, due to the versatility of the real-time description in Octopus, the underlying excitation does not need to be originating from an optical pulse (i.e., an external vector potential), and since the ions can move according to TDDFT-Ehrenfest dynamics,¹³⁹ one can also simulate features coming from the ionic or lattice motion, specifically electron-phonon coupling signature in ARPES.¹⁴⁰ Pump-probe simulations are not limited to the dynamical process of excitations but can also be employed to study steady-state modifications of driven electronic states.¹⁴¹ This opens the possibility of studying Floquet physics from an ab initio perspective and to simulate directly the effect that a periodic force would have on the dressed electronic structure—an important aspect to underline, given the growing interest in the field of Floquet engineering¹⁴² and Floquet analysis.¹⁴³

Finally, t-SURFF is particularly suited to simulating strong field ionization processes such as in laser-induced photoelectron diffraction (LIED) experiments, where the ionization takes place by direct tunneling into the continuum and the field is so strong as to drive electrons in trajectories recolliding with the parent system. In fact, since in Eq. (48), we accumulate the flux through S over the time of the propagation, with t-SURFF, electrons can seamlessly flow back and forth through the surface without producing any artifact in the final spectrum. As an example, in Fig. 11(d), we present the photoelectron angular distribution of C₆₀ ionized by a strong IR pulse capable of inducing rescattering dynamics, which then gets imprinted in the photoelectron spectrum as characteristic rings centered in the value of the vector potential at the instant of

rescattering.¹⁴⁴ In this regime, the result of simulations obtained with Octopus are in excellent agreement with the experiments.¹⁴⁵,¹⁴⁶

## XI. ELECTRIC AND THERMAL CONDUCTIVITIES

Electrical conductivity in real materials can be described to a first approximation by Ohm's law, which is a linear relationship between the applied electric field and the current density generated in response. The standard method to study electrical and thermal conductivity is to carry out a DFT calculation at equilibrium and apply Kubo-Greenwood linear-response theory to evaluate the conductivities. Although this method has been used successfully in the past to calculate the transport properties,¹⁴⁷,¹⁴⁸ the application is limited to systems in the linear response regime. It is known that in many cases, especially in the presence of strong external fields, Ohm's law is no longer valid and the system exhibits non-linear behavior. In order to capture this complex behavior in materials, one must go beyond simple, linear approximations and describe electron interactions in a non-trivial manner.

One promising route to go beyond linear response is to use density functional methods to study directly thermoelectric transport. This topic has received less attention but has recently been implemented in Octopus and applied to liquid aluminum.¹⁴⁹ For detailed descriptions and derivations of non-equilibrium thermoelectric phenomena using density based methods, we refer the reader to Ref. 150.

In Octopus, we are able to calculate the current density and heat current density at each step during a time-dependent simulation. A current is induced in this study by applying an electric field of the form $E(t) = E_0\delta(t)$, where $E_0$ describes the magnitude and direction of the field. The electric field is induced through a time-dependent vector potential $[E(t) = -c^{-1}\partial_t A(t)$, where c is the speed of light] in the Hamiltonian. This vector potential satisfies the periodic boundary conditions of an extended system. It should be noted that there is no limitation on the form of the applied electric field in general. We are able to evaluate the macroscopic current density as

$$J(t) = -\frac{i}{\Omega} \int dr \sum_{j}^{N} \varphi^{\dagger}(\mathbf{r}, t) [\hat{H}(t), \hat{r}] \varphi_{i}(\mathbf{r}, t), \quad (50)$$

where $\Omega$ is the volume of the unit cell. The energy current density, $\hat{J}_n(\mathbf{r}, t)$, is expressed as

$$\hat{J}_n(\mathbf{r}, t) = \hat{J}_t(\mathbf{r}, t) + \hat{J}_v(\mathbf{r}, t) + \hat{J}_u(\mathbf{r}, t) + \hat{J}_f(\mathbf{r}, t), \quad (51)$$

where $\hat{J}_t(\mathbf{r})$ is the kinetic energy contribution given as

$$[\hat{J}_t(\mathbf{r}, t)]_i = \frac{i}{8} \left( [\partial_t \hat{\phi}^\dagger] [\nabla^2 \hat{\phi}] - [\nabla^2 \hat{\phi}^\dagger] [\partial_t \hat{\phi}] - [\partial_t \nabla \hat{\phi}^\dagger] \cdot [\nabla \hat{\phi}] - [\nabla \hat{\phi}^\dagger] \cdot [\partial_t \nabla \hat{\phi}] \right) \quad (52)$$

and $\varphi = \varphi(\mathbf{r}, t)$ are understood to be the time-dependent states. The output of the time-dependent simulation can be further analyzed to yield the frequency-dependent conductivity. The electrical (or thermal) conductivity $\sigma$ can be found by Fourier transforming the corresponding current as

$$\sigma_{ij}(\omega) = \frac{1}{[\mathbf{E}_0]_i} \int_0^\infty dt e^{-i\omega t} \mathbf{J}_j(t). \quad (53)$$

In general, it is possible to study the conductivity of an extended system in the presence of an applied electric field. Here, we
<figure><img src="image_22.png" /><figcaption>FIG. 12. (Left) Time-dependent current density for one ionic snapshot of a supercell of 128 atoms of hydrogen in the liquid phase at 1400 K and 400 GPa with an initial electric field $E_0$ of 0.1 a.u. The blue line represents the current density as given in Eq. (50), and the red line is the heat current evaluated from Eq. (52). Both currents are decaying to zero. (Right) The Fourier transform of the respective currents gives the frequency-dependent conductivities.</figcaption></figure>

illustrate the use of this method on hydrogen at 1400 K and 400 GPa. At this temperature and pressure, hydrogen is in its liquid metallic phase. A molecular dynamics simulation was carried out with VASP¹⁵¹ to produce several ionic snapshots of the system. For each ionic snapshot, we performed a TDDFT simulation with a time step of 0.05 a.u. (0.001 21 fs) for a total time of 3 fs. The initial electric field $E_0$ was 0.1 a.u. The calculation was carried out on a 3 × 3 × 3 k-point grid to ensure that the current density decays to zero.

In Fig. 12, we plot the time-dependent current density and heat current density as a function of time after an initial electric field is applied at time zero. Both currents are shown to decay to zero by the end of the simulation. It should be noted that the heat current corresponds to the kinetic energy contribution given in Eq. (52). In general, this heat current is related to the Peltier coefficient, since it describes the heat generated from the response to an electric current. The remaining contributions in Eq. (51) are not yet implemented in Octopus. Both $J_u$ and $J_f$ arise from electron-electron interactions and $J_v$ is the potential energy. In future work, it would be interesting to evaluate these remaining contributions to determine if they are large or can reasonably be ignored for some systems. The electrical conductivity and the heat conductivity can be evaluated using Eq. (53). These frequency-dependent conductivities are shown in Fig. 12 (right). One may also calculate the DC conductivity by taking $\omega = 0$.

This suite of tools implemented in Octopus will allow for the study of time-dependent thermoelectric phenomena using TDDFT. This approach to study current and conductivity is more general and widely applicable than the standard Kubo-Greenwood approach. While Kubo-Greenwood is a linear response theory, this method can be applied to study materials where non-linear conduction effects can appear. In a recent study, this TDDFT method was used to illustrate non-linear conductivity effects in liquid aluminum,¹⁴⁹ which would not be possible by applying standard linear response theory. In the future, we plan to extend these tools by implementing a thermal vector potential that can induce a heat current in an extended system during a time-dependent simulation.

## XII. LOCAL DOMAIN CONTRIBUTION TO PHYSICAL OBSERVABLES

Usually, we are interested in studying systems consisting of different atoms, molecules, regions, or domains. In such cases, we might want to understand the different contributions from different parts of the system toward a specific observable. Based on the

Hohenberg-Kohn theorems for the DFT,¹⁵² and its extension for TDDFT, the Runge-Gross theorem,¹⁵³ we can state that any physical observable of the system is a functional of the electronic density, either static (ground state) or time dependent. In other words, the expectation value of an operator $\hat{O}$ can be expressed as a functional of the electronic density

$$ \mathcal{O}[n] = \langle \Psi[n] | \hat{O} | \Psi[n] \rangle. \qquad (54) $$

Let us now consider the case where we write the total electronic density as a sum of $N$ densities

$$ n(\mathbf{r}) = \sum_{i}^{N} n_{i}(\mathbf{r}), \qquad (55) $$

and no further constraints are imposed on $n$ or $n_i$. Such partitioning of the density is at the basis of subsystem DFT, which allows us to divide a system into several Kohn-Sham subsystems that interact with each other in a theoretically well justified manner (see Ref. 154 and references therein). However, in our case, we are only interested in this kind of partition for post-processing purposes where we assign a density $n_i$ to a specific part of the system in order to identify how it contributes to a given observable. We are, thus, interested in operators for which the following condition is true:

$$ \mathcal{O}[n] = \sum_{i}^{N} \mathcal{O}[n_i]. \qquad (56) $$

Such operators are said to be additive. One example of an additive operator that is particularly relevant in the context of TDDFT is the time-dependent dipole

$$ d(t) = \int \mathbf{r}n(\mathbf{r}, t)d\mathbf{r} = \sum_{i}^{N} \int \mathbf{r}n_{i}(\mathbf{r}, t)d\mathbf{r}, \qquad (57) $$

as this is the relevant observable for obtaining the optical absorption cross section of finite systems.

The range of operators that fulfill Eq. (56) can be further expanded if we now consider non-overlapping densities, that is, densities $n_i$ that are non-zero only in a given domain $V_i$ and that the domains do not overlap. We refer to this type of partitioning as a local domain partitioning. In such cases, any (semi-)local operator that depends (semi-)locally on the density should fulfill Eq. (56). In this context, an example of a relevant observable is the exchange-correlation energy within the LDA (or GGA) approximation
<figure><img src="image_24.png" /><figcaption>FIG. 12. (Left) Time-dependent current density for one ionic snapshot of a supercell of 128 atoms of hydrogen in the liquid phase at 1400 K and 400 GPa with an initial electric field $E_0$ of 0.1 a.u. The blue line represents the current density as given in Eq. (50), and the red line is the heat current evaluated from Eq. (52). Both currents are decaying to zero. (Right) The Fourier transform of the respective currents gives the frequency-dependent conductivities.</figcaption></figure>

illustrate the use of this method on hydrogen at 1400 K and 400 GPa. At this temperature and pressure, hydrogen is in its liquid metallic phase. A molecular dynamics simulation was carried out with VASP¹⁵¹ to produce several ionic snapshots of the system. For each ionic snapshot, we performed a TDDFT simulation with a time step of 0.05 a.u. (0.001 21 fs) for a total time of 3 fs. The initial electric field $E_0$ was 0.1 a.u. The calculation was carried out on a 3 × 3 × 3 k-point grid to ensure that the current density decays to zero.

In Fig. 12, we plot the time-dependent current density and heat current density as a function of time after an initial electric field is applied at time zero. Both currents are shown to decay to zero by the end of the simulation. It should be noted that the heat current corresponds to the kinetic energy contribution given in Eq. (52). In general, this heat current is related to the Peltier coefficient, since it describes the heat generated from the response to an electric current. The remaining contributions in Eq. (51) are not yet implemented in Octopus. Both $J_u$ and $J_f$ arise from electron-electron interactions and $J_v$ is the potential energy. In future work, it would be interesting to evaluate these remaining contributions to determine if they are large or can reasonably be ignored for some systems. The electrical conductivity and the heat conductivity can be evaluated using Eq. (53). These frequency-dependent conductivities are shown in Fig. 12 (right). One may also calculate the DC conductivity by taking $\omega = 0$.

This suite of tools implemented in Octopus will allow for the study of time-dependent thermoelectric phenomena using TDDFT. This approach to study current and conductivity is more general and widely applicable than the standard Kubo-Greenwood approach. While Kubo-Greenwood is a linear response theory, this method can be applied to study materials where non-linear conduction effects can appear. In a recent study, this TDDFT method was used to illustrate non-linear conductivity effects in liquid aluminum,¹⁴⁹ which would not be possible by applying standard linear response theory. In the future, we plan to extend these tools by implementing a thermal vector potential that can induce a heat current in an extended system during a time-dependent simulation.

## XII. LOCAL DOMAIN CONTRIBUTION TO PHYSICAL OBSERVABLES

Usually, we are interested in studying systems consisting of different atoms, molecules, regions, or domains. In such cases, we might want to understand the different contributions from different parts of the system toward a specific observable. Based on the

Hohenberg-Kohn theorems for the DFT,¹⁵² and its extension for TDDFT, the Runge-Gross theorem,¹⁵³ we can state that any physical observable of the system is a functional of the electronic density, either static (ground state) or time dependent. In other words, the expectation value of an operator $\hat{O}$ can be expressed as a functional of the electronic density

$$ \mathcal{O}[n] = \langle \Psi[n] | \hat{O} | \Psi[n] \rangle. \qquad (54) $$

Let us now consider the case where we write the total electronic density as a sum of $N$ densities

$$ n(\mathbf{r}) = \sum_{i}^{N} n_{i}(\mathbf{r}), \qquad (55) $$

and no further constraints are imposed on $n$ or $n_i$. Such partitioning of the density is at the basis of subsystem DFT, which allows us to divide a system into several Kohn-Sham subsystems that interact with each other in a theoretically well justified manner (see Ref. 154 and references therein). However, in our case, we are only interested in this kind of partition for post-processing purposes where we assign a density $n_i$ to a specific part of the system in order to identify how it contributes to a given observable. We are, thus, interested in operators for which the following condition is true:

$$ \mathcal{O}[n] = \sum_{i}^{N} \mathcal{O}[n_i]. \qquad (56) $$

Such operators are said to be additive. One example of an additive operator that is particularly relevant in the context of TDDFT is the time-dependent dipole

$$ d(t) = \int \mathbf{r}n(\mathbf{r}, t)d\mathbf{r} = \sum_{i}^{N} \int \mathbf{r}n_{i}(\mathbf{r}, t)d\mathbf{r}, \qquad (57) $$

as this is the relevant observable for obtaining the optical absorption cross section of finite systems.

The range of operators that fulfill Eq. (56) can be further expanded if we now consider non-overlapping densities, that is, densities $n_i$ that are non-zero only in a given domain $V_i$ and that the domains do not overlap. We refer to this type of partitioning as a local domain partitioning. In such cases, any (semi-)local operator that depends (semi-)locally on the density should fulfill Eq. (56). In this context, an example of a relevant observable is the exchange-correlation energy within the LDA (or GGA) approximation
$$E_{xc}^{\text{LDA}}[n] = \int n(\mathbf{r})e_{xc}^{\text{LDA}}(n(\mathbf{r}))d\mathbf{r} = \sum_{i}^{N} \int_{V_i} n_i(\mathbf{r})e_{xc}^{\text{LDA}}(n_i(\mathbf{r}))d\mathbf{r}, \quad (58)$$

where $e_{xc}^{\text{LDA}}$ is the exchange-correlation energy per unit particle.

Currently, several different options can be used in Octopus to define specific regions of the simulation box. These options include simple geometric shapes, such as spheres centered on particular points of space, atom-dependent domains, such as a union of spheres centered on the atoms, or the definition of Bader volumes. The latter option follows the Quantum Theory of Atoms in Molecules (QTAIM)¹⁵⁵ and associates a density region with a given atom through a density gradient path such that the boundaries of each volume are defined as the surfaces through which the charge density gradient has a zero flux. Note that some of these options allow for the user to specify overlapping regions. In this case, the overlap must be taken into account when analyzing the results and extra care is required when comparing results from different domains.

Let us exemplify the applicability of the local domain partitioning by decomposing the optical response for a coupled chromophore system. Similar to Frozen Density Embedding real-time TDDFT (FDE-rt-TDDFT),¹⁵⁶ we decompose the total optical spectrum as a sum of the local response within each domain

$$\alpha(\omega) = \sum_{i} \alpha_{i}(\omega), \quad (59)$$

where $\alpha$ is the dynamic polarizability. In addition, from each local dynamic polarizability tensor, we can compute for each fragment the

corresponding cross section. This is justified by the fact that the relevant operator to calculate $\alpha$ is the dipole operator, and as shown above, this is an additive operator that fulfills Eq. (56).

As an example, we performed a series of optical spectrum simulations using real-time TDDFT for a benzene-fulvene dimer with different $\pi$-stacking separation ranging from 4 Å to 8 Å. We use the standard PBE exchange-correlation functional¹⁵⁷,¹⁵⁸ and the optimized norm-conserving Vanderbilt PBE pseudopotential (sg15) set.¹⁵⁹ The real-space grid is defined as a parallelepiped box with length 16 Å, 17 Å, and 20 Å in the three Cartesian axes. The spacing between points is set to be 0.13 Å.

For each intermolecular separation, we carry out a single ground-state calculation and three time-propagations. For each time-propagation, a dipolar electric perturbation is applied along one of the Cartesian axes.¹²⁵ We let the perturbed Kohn-Sham states evolve for a propagation time of $T = 24 \, \text{ħ/eV}$ (15.8 fs) with a resolution of 0.26 eV ($2\pi/T$). The ground-state electron density is fragmented following the Bader atomic decomposition, and each molecular domain is defined as the sum of these atomic volumes. Then, the corresponding dipole operator is applied over each defined domain. Finally, by Fourier transform of the local time-dependent dipole moment, the local polarizability tensor is recovered.

Figure 13 shows the spectral decomposition of the benzene-fulvene system as a function of the intermolecular separation. The full system at large intermolecular separation (8 Å) shows two major peaks located at the same frequencies as for the isolated molecules. As the intermolecular distance becomes smaller, the global spectrum changes due to electrostatic effects induced by the electronic cloud

<figure><img src="image_23.png" /><figcaption>FIG. 13. (Top) Schematic representation of the benzene-fulvene dimer with $\pi$-stacking along the z axis and an intermolecular distance of 4 Å. From left to right, the electronic density for the global system for an isosurface of 0.006 Å⁻³ and the QTAIM local density for the fulvene and benzene molecules, respectively. Visualized with USCF Chimera software.¹⁶⁰ (Bottom) Photo-absorption cross section obtained from, from left to right, the full system density, the fulvene local density, and the benzene local density.</figcaption></figure>
$$E_{xc}^{\text{LDA}}[n] = \int n(\mathbf{r})e_{xc}^{\text{LDA}}(n(\mathbf{r}))d\mathbf{r} = \sum_{i}^{N} \int_{V_i} n_i(\mathbf{r})e_{xc}^{\text{LDA}}(n_i(\mathbf{r}))d\mathbf{r}, \quad (58)$$

where $e_{xc}^{\text{LDA}}$ is the exchange-correlation energy per unit particle.

Currently, several different options can be used in Octopus to define specific regions of the simulation box. These options include simple geometric shapes, such as spheres centered on particular points of space, atom-dependent domains, such as a union of spheres centered on the atoms, or the definition of Bader volumes. The latter option follows the Quantum Theory of Atoms in Molecules (QTAIM)¹⁵⁵ and associates a density region with a given atom through a density gradient path such that the boundaries of each volume are defined as the surfaces through which the charge density gradient has a zero flux. Note that some of these options allow for the user to specify overlapping regions. In this case, the overlap must be taken into account when analyzing the results and extra care is required when comparing results from different domains.

Let us exemplify the applicability of the local domain partitioning by decomposing the optical response for a coupled chromophore system. Similar to Frozen Density Embedding real-time TDDFT (FDE-rt-TDDFT),¹⁵⁶ we decompose the total optical spectrum as a sum of the local response within each domain

$$\alpha(\omega) = \sum_{i} \alpha_{i}(\omega), \quad (59)$$

where $\alpha$ is the dynamic polarizability. In addition, from each local dynamic polarizability tensor, we can compute for each fragment the

corresponding cross section. This is justified by the fact that the relevant operator to calculate $\alpha$ is the dipole operator, and as shown above, this is an additive operator that fulfills Eq. (56).

As an example, we performed a series of optical spectrum simulations using real-time TDDFT for a benzene-fulvene dimer with different $\pi$-stacking separation ranging from 4 Å to 8 Å. We use the standard PBE exchange-correlation functional¹⁵⁷,¹⁵⁸ and the optimized norm-conserving Vanderbilt PBE pseudopotential (sg15) set.¹⁵⁹ The real-space grid is defined as a parallelepiped box with length 16 Å, 17 Å, and 20 Å in the three Cartesian axes. The spacing between points is set to be 0.13 Å.

For each intermolecular separation, we carry out a single ground-state calculation and three time-propagations. For each time-propagation, a dipolar electric perturbation is applied along one of the Cartesian axes.¹²⁵ We let the perturbed Kohn-Sham states evolve for a propagation time of $T = 24 \text{ h/eV}$ (15.8 fs) with a resolution of 0.26 eV (2$\pi/T$). The ground-state electron density is fragmented following the Bader atomic decomposition, and each molecular domain is defined as the sum of these atomic volumes. Then, the corresponding dipole operator is applied over each defined domain. Finally, by Fourier transform of the local time-dependent dipole moment, the local polarizability tensor is recovered.

Figure 13 shows the spectral decomposition of the benzene-fulvene system as a function of the intermolecular separation. The full system at large intermolecular separation (8 Å) shows two major peaks located at the same frequencies as for the isolated molecules. As the intermolecular distance becomes smaller, the global spectrum changes due to electrostatic effects induced by the electronic cloud

<figure><img src="image_25.png" /><figcaption>FIG. 13. (Top) Schematic representation of the benzene-fulvene dimer with $\pi$-stacking along the z axis and an intermolecular distance of 4 Å. From left to right, the electronic density for the global system for an isosurface of 0.006 Å⁻³ and the QTAIM local density for the fulvene and benzene molecules, respectively. Visualized with USCF Chimera software.¹⁶⁰ (Bottom) Photo-absorption cross section obtained from, from left to right, the full system density, the fulvene local density, and the benzene local density.</figcaption></figure>
of the neighboring molecule. We can see that the peak located at around 5 eV suffers a red-shift, while the stronger peak between 6.5 eV and 7 eV reduces its intensity, giving rise to an oscillator strength transfer to excited states at higher energy. The local domain analysis reveals that the major change is caused by the reduction of the oscillator strength of the fulvene peak located at around 6.5 eV.

These results are in agreement with the FDE-rt-TDDFT calculations of Ref. 156, further validating our strategy. In addition, the local domain methodology does not require any previous fragmentation, selection, or localization of the basis set, allowing us also to treat very large systems, such as biological molecules like the major light harvesting complex LHC-II.161 Recently, we have shown that this technique also allows for exciton coupling calculations when combined with a new formulation of the transition densities from real-time TDDFT.162

### XIII. NEW PROPAGATORS FOR REAL-TIME TDDFT

At the core of most Octopus calculations lies the need to integrate the time-dependent Kohn–Sham (TDKS) equations. These are a set of non-linear equations, given the dependence of the Kohn–Sham Hamiltonian on the electronic density. Upon the discretization of the electronic Hilbert space, they take the generic form of a system of first-order ordinary differential equations (ODEs)

$$ \dot{\varphi} = f(\varphi(t), t), \quad (60) $$

$$ \varphi(t_0) = \varphi^0, \quad (61) $$

where $t_0$ is the initial time, $\varphi$ is an array containing all the Kohn–Sham orbitals, $\{\varphi_m^N\}$ ($\varphi_0$ is its initial value), and $f$ is a vector function ($f = (f_1, \dots, f_N)$) given by the action of the Kohn–Sham Hamiltonian $\hat{H}[n(t), t]$,

$$ f_i(\varphi(t), t) = -i\hat{H}[n(t)t]\varphi_i(t), \quad (i = 1, \dots, N). \quad (62) $$

Note that (i) if the nuclei are also to be propagated, the state should be supplemented with their position and momenta, and the equations with the corresponding nuclear equations of motion. (ii) Likewise, the formalism for solids includes a polarization vector field22 that must also be included in the system definition and propagation. (iii) Strictly speaking, the TDKS equations are not ordinary but “delay differential equations” (i.e., equations for which the derivative of the unknown function at a certain time depends on the function values at previous times), due to the dependence of the exact exchange and correlation potential functional on the past densities. This is ignored in the adiabatic approximation, which is almost always assumed.

A myriad of numerical propagators for ODEs are available, all of them theoretically applicable to the TDKS equations. Ideally, however, one should choose a propagator that respects all the mathematical properties of the equations that describe the problem at hand, for example, the preservation of the norm of the orbitals. In the case of the TDKS equations in the adiabatic approximation, another such property is symplecticity (a fact that is demonstrated in Ref. 163). This property is usually stated in terms of the conservation of the volume of the flow of the system of ODEs in the phase space, although the precise definition is as follows:

Any system of ODEs can be viewed as a “flow,” a differentiable map $g: \mathbb{R}^p \to \mathbb{R}^p$ that transforms the state $y(t_0)$ at some point in time $t_0$ into the state at time $t$, $y(t) = g(y(t_0))$. For systems described with complex variables such as ours, since it can be split into a real and an imaginary part, $p$ is even: $g: \mathbb{R}^{2N} \to \mathbb{R}^{2N}$. A map with an even number of variables such as this one is defined to be symplectic if and only if

$$ \frac{\partial g^T}{\partial x} J \frac{\partial g}{\partial x} = J, \quad J = \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix}, \quad (63) $$

where $I$ is the unit matrix of dimension $N$, and $x \in \mathbb{R}^{2N}$ ( conventionally, the first $N$ variables of $x$ are called the “coordinates” and the second half are the “momenta,” but no physical meaning should be assumed for them in this purely mathematical definition). A numerical propagator is also a differentiable map, which relates the solution $y(t)$ to $y(t + \Delta t)$. As such, it may respect or not the symplecticity—and other properties—of the original flow.

The symplecticity also means that the system has to be “Hamiltonian”: one may find a set of coordinates $q \in \mathbb{R}^N$, $p \in \mathbb{R}^N$ (in this case, the real and imaginary parts of the orbitals coefficients) and some scalar function $H(q, p)$ (the Hamiltonian expectation value), which permits us to rewrite the system into the well-known form of Hamilton’s equations

$$ \dot{q}_i = \frac{\partial H(p, q)}{\partial p_i}, \quad (64a) $$

$$ \dot{p}_i = - \frac{\partial H(p, q)}{\partial q_i}. \quad (64b) $$

The relevance of the symplectic numerical propagators stems from the fact that they present some features, such as a better conservation of the energy at long times, where the value of the energy ends up oscillating around the true value instead of diverging. For a more detailed discussion about symplecticity on numerical propagators, see Ref. 164.

The Octopus code has several propagator options, many of them already described in Ref. 165, such as the Crank–Nicolson, standard Runge–Kutta, exponential midpoint rule (EMR) and variations (e.g., the “enforced time-reversal symmetry” scheme), and split operator techniques. In a recent paper,163 some of the present authors have studied propagators of different families that had been scarcely (or not at all) tested for the TDKS equations: multistep, exponential Runge–Kutta, and commutator-free Magnus (CFM) expansions. After considering the accuracy, stability, and the performance of the propagators, the CFM techniques were identified as suitable schemes for TDDFT problems and implemented in Octopus. In this section, we make a brief description.

Developed in Ref. 166 for linear non-autonomous systems, the CFM expansion offers an alternative to the “standard” Magnus expansion, which requires expensive application of nested commutators of the Hamiltonian with itself at different times. In essence, CFM expansions $\Gamma(t + \Delta t, t)$ consist of substituting the propagator $\hat{U}(t + \Delta t, t)$ by a product of exponentials

$$ \Gamma(t + \Delta t, t) = \prod_{i=1}^{m} \exp(\hat{D}_i). \quad (65) $$

The $m$ linear operators $\hat{D}_i$ are either the Hamiltonian at different times within the interval $[t, t + \Delta t]$ or parts of it.
of the neighboring molecule. We can see that the peak located at around 5 eV suffers a red-shift, while the stronger peak between 6.5 eV and 7 eV reduces its intensity, giving rise to an oscillator strength transfer to excited states at higher energy. The local domain analysis reveals that the major change is caused by the reduction of the oscillator strength of the fulvene peak located at around 6.5 eV.

These results are in agreement with the FDE-rt-TDDFT calculations of Ref. 156, further validating our strategy. In addition, the local domain methodology does not require any previous fragmentation, selection, or localization of the basis set, allowing us also to treat very large systems, such as biological molecules like the major light harvesting complex LHC-II.161 Recently, we have shown that this technique also allows for exciton coupling calculations when combined with a new formulation of the transition densities from real-time TDDFT.162

### XIII. NEW PROPAGATORS FOR REAL-TIME TDDFT

At the core of most Octopus calculations lies the need to integrate the time-dependent Kohn-Sham (TDKS) equations. These are a set of non-linear equations, given the dependence of the Kohn-Sham Hamiltonian on the electronic density. Upon the discretization of the electronic Hilbert space, they take the generic form of a system of first-order ordinary differential equations (ODEs)

$$ \dot{\varphi} = f(\varphi(t), t), \quad (60) $$

$$ \varphi(t_0) = \varphi^0, \quad (61) $$

where $t_0$ is the initial time, $\varphi$ is an array containing all the Kohn-Sham orbitals, $\{\varphi_m^N\}$ ($\varphi_0$ is its initial value), and $f$ is a vector function ($f = (f_1, \dots, f_N)$) given by the action of the Kohn-Sham Hamiltonian $\hat{H}[n(t), t]$,

$$ f_i(\varphi(t), t) = -i\hat{H}[n(t)t]\varphi_i(t), \quad (i = 1, \dots, N). \quad (62) $$

Note that (i) if the nuclei are also to be propagated, the state should be supplemented with their position and momenta, and the equations with the corresponding nuclear equations of motion. (ii) Likewise, the formalism for solids includes a polarization vector field22 that must also be included in the system definition and propagation. (iii) Strictly speaking, the TDKS equations are not ordinary but “delay differential equations” (i.e., equations for which the derivative of the unknown function at a certain time depends on the function values at previous times), due to the dependence of the exact exchange and correlation potential functional on the past densities. This is ignored in the adiabatic approximation, which is almost always assumed.

A myriad of numerical propagators for ODEs are available, all of them theoretically applicable to the TDKS equations. Ideally, however, one should choose a propagator that respects all the mathematical properties of the equations that describe the problem at hand, for example, the preservation of the norm of the orbitals. In the case of the TDKS equations in the adiabatic approximation, another such property is symplecticity (a fact that is demonstrated in Ref. 163). This property is usually stated in terms of the conservation of the volume of the flow of the system of ODEs in the phase space, although the precise definition is as follows:

Any system of ODEs can be viewed as a “flow,” a differentiable map $g: \mathbb{R}^p \to \mathbb{R}^p$ that transforms the state $y(t_0)$ at some point in time $t_0$ into the state at time $t$, $y(t) = g(y(t_0))$. For systems described with complex variables such as ours, since it can be split into a real and an imaginary part, $p$ is even: $g: \mathbb{R}^{2N} \to \mathbb{R}^{2N}$. A map with an even number of variables such as this one is defined to be symplectic if and only if

$$ \frac{\partial g^T}{\partial x} J \frac{\partial g}{\partial x} = J, \quad J = \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix}, \quad (63) $$

where $I$ is the unit matrix of dimension $N$, and $x \in \mathbb{R}^{2N}$ (conventionally, the first $N$ variables of $x$ are called the “coordinates” and the second half are the “momenta,” but no physical meaning should be assumed for them in this purely mathematical definition). A numerical propagator is also a differentiable map, which relates the solution $y(t)$ to $y(t + \Delta t)$. As such, it may respect or not the symplecticity—and other properties—of the original flow.

The symplecticity also means that the system has to be “Hamiltonian”: one may find a set of coordinates $q \in \mathbb{R}^N$, $p \in \mathbb{R}^N$ (in this case, the real and imaginary parts of the orbitals coefficients) and some scalar function $H(q, p)$ (the Hamiltonian expectation value), which permits us to rewrite the system into the well-known form of Hamilton’s equations

$$ \dot{q}_i = \frac{\partial H(p, q)}{\partial p_i}, \quad (64a) $$

$$ \dot{p}_i = - \frac{\partial H(p, q)}{\partial q_i}. \quad (64b) $$

The relevance of the symplectic numerical propagators stems from the fact that they present some features, such as a better conservation of the energy at long times, where the value of the energy ends up oscillating around the true value instead of diverging. For a more detailed discussion about symplecticity on numerical propagators, see Ref. 164.

The Octopus code has several propagator options, many of them already described in Ref. 165, such as the Crank-Nicolson, standard Runge-Kutta, exponential midpoint rule (EMR) and variations (e.g., the “enforced time-reversal symmetry” scheme), and split operator techniques. In a recent paper,163 some of the present authors have studied propagators of different families that had been scarcely (or not at all) tested for the TDKS equations: multistep, exponential Runge-Kutta, and commutator-free Magnus (CFM) expansions. After considering the accuracy, stability, and the performance of the propagators, the CFM techniques were identified as suitable schemes for TDDFT problems and implemented in Octopus. In this section, we make a brief description.

Developed in Ref. 166 for linear non-autonomous systems, the CFM expansion offers an alternative to the “standard” Magnus expansion, which requires expensive application of nested commutators of the Hamiltonian with itself at different times. In essence, CFM expansions $\Gamma(t + \Delta t, t)$ consist of substituting the propagator $\hat{U}(t + \Delta t, t)$ by a product of exponentials

$$ \Gamma(t + \Delta t, t) = \prod_{i=1}^{m} \exp(\hat{D}_i). \quad (65) $$

The $m$ linear operators $\hat{D}_i$ are either the Hamiltonian at different times within the interval $[t, t + \Delta t]$ or parts of it.
We implemented in Octopus an order four ($q = 4$) version, given, for example, in Eq. (43) of Ref. 166, which we call hereafter CFM4. This propagator requires two exponentials ($m = 2$),

$$ \varphi(t + \Delta t) = \exp(-i\Delta t(\alpha_1 \hat{H}[t_1] + \alpha_2 \hat{H}[t_2])) \times \exp(-i\Delta t(\alpha_2 \hat{H}[t_1] + \alpha_1 \hat{H}[t_2]))\varphi(t), \quad (66) $$

where

$$ \alpha_1 = \frac{3 - 2\sqrt{3}}{12}, \qquad t_1 = t + \left(\frac{1}{2} - \frac{\sqrt{3}}{6}\right)\Delta t, \quad (67) $$

$$ \alpha_2 = \frac{3 + 2\sqrt{3}}{12}, \qquad t_2 = t + \left(\frac{1}{2} + \frac{\sqrt{3}}{6}\right)\Delta t. \quad (68) $$

$\hat{H}[t_1]$ and $\hat{H}[t_2]$ are the Hamiltonians at times $t_1$ and $t_2$, which are, in fact, unknown because they depend on the Kohn–Sham states through the density: we are dealing with a non-linear problem and the CFM expansions were, in fact, developed for linear systems. We have various options to extend them for our non-linear problem: for example, one could define $\hat{H}[t_1]$ and $\hat{H}[t_2]$ as interpolated Hamiltonians from $\hat{H}[t]$ and $\hat{H}[t + \Delta t]$ in which case we end up with an implicit equation for $\varphi(t + \Delta t)$ that we would have to solve at a substantial cost.

The alternative that we implemented, however, is to approximate $\hat{H}[t_1]$ via an extrapolation from the Hamiltonian at various previous time steps (in practice, it is the Hartree, exchange, and correlation parts that must be extrapolated). The resulting method is then explicit, i.e., no linear or non-linear algebraic equations need to be solved. The fourth order accuracy is preserved as long as the extrapolation is also done at order four.

As an example, Fig. 14 shows the performance of this CFM4 method against the well-known exponential midpoint rule (EMR). The benchmark system consists of a benzene molecule. It is subject to an instantaneous perturbation at the beginning of the propagation, and then, it evolves freely for some fixed interval of time. The propagations are performed at varying values of $\Delta t$. The plot displays the cost of the propagation vs the accuracy achieved in each

<figure><img src="image_24.png" /><figcaption>FIG. 14. Wall-time computational cost of the method (in seconds) against the error in the wave function [defined as $Error(T, \Delta t) = \sqrt{\sum_m ||\varphi_m(T) - \varphi_m^{exact}(T)||^2}$, where the $\varphi_m^{exact}$ are the values of the Kohn–Sham orbitals computed using the explicit fourth-order Runge–Kutta integrator with a very small time step] for the EMR and CFM4 propagators.</figcaption></figure>

case. As we can see, the CFM4 outperforms the EMR for all values of the error examined, although the differences become more obvious when higher accuracy is demanded. The key difference between the two methods is the fourth-order accuracy of the CFM4 scheme—at the cost of requiring two exponentials, instead of just one, while the EMR is only second-order accurate. This is reflected in the different slopes of the curves as the error becomes smaller, i.e., as $\Delta t \to 0$.

## XIV. CONJUGATE GRADIENT IMPLEMENTATION IN RDMFT

The RDMFT implementation in Octopus has been described in detail in a previous paper.¹⁹ In this section, we review briefly the existing implementation, providing some new insights, and introduce a recently implemented method to solve the RDMFT equations, which is better suited to the real-space grids used in the code.

The optimization of the natural orbitals in RDMFT is subject to an orthonormalization constraint for the orbitals. One minimizes the functional

$$ E_{\text{total}} - \sum_{j,k=1}^{M} \lambda_{jk} \left( \int d\mathbf{r} \psi_j^*(\mathbf{r}) \psi_k(\mathbf{r}) - \delta_{jk} \right), \quad (69) $$

where $E_{total}$ and $M$ denote the total energy and the number of natural orbitals $\psi_j$ included in the calculation, respectively. The parameters $\lambda_{jk}$ are the Lagrange multipliers, which ensure the orthonormality of the natural orbitals at the solution point. We note that, for an RDMFT calculation, the number of natural orbitals has to be larger than the number of electrons in the system (core electrons that are included in the pseudopotential do not count). The exact number $M$ is system-dependent and should be treated as an additional parameter with respect to which a convergence study should be carried out, just like it is done for the basis set or the grid parameters. Typically, the optimization is performed using the so-called Piris method,¹⁶⁷ which was the method previously implemented in Octopus. Within this approach, one uses the orthonormality constraint of the natural orbitals, which implies that a certain matrix constructed from the Lagrange multipliers $\lambda_{jk}$ is diagonal at the solution point. As an immediate consequence of the Piris method, the natural orbitals at the solution point are linear combinations of the orbitals used as a starting point for the minimization. In other words, the initial orbitals serve as a basis. Consequently, the necessary matrix elements for different energy contributions can be calculated for the basis functions (initial orbitals) before the iterative optimization of the natural orbitals and their occupation numbers is started. In addition, the optimization of the occupation numbers can be turned off completely, resulting in a Hartree–Fock calculation on the basis of the initial orbitals.

For the existing implementation in Octopus of the Piris method, the initial orbitals are taken to be the solutions obtained with a different level of theory, such as independent particles or DFT. In order to better understand the effect of the choice of basis, we tested the following choices: (i) independent particles, density functional theory within (ii) the local density approximation (LDA) or (iii) the exact exchange (EXX) approximation, and (iv) the Hartree–Fock approximation. In all cases, we have to ensure that the number
We implemented in Octopus an order four ($q = 4$) version, given, for example, in Eq. (43) of Ref. 166, which we call hereafter CFM4. This propagator requires two exponentials ($m = 2$),

$$ \varphi(t + \Delta t) = \exp(-i\Delta t(\alpha_1 \hat{H}[t_1] + \alpha_2 \hat{H}[t_2])) \times \exp(-i\Delta t(\alpha_2 \hat{H}[t_1] + \alpha_1 \hat{H}[t_2]))\varphi(t), \quad (66) $$

where

$$ \alpha_1 = \frac{3 - 2\sqrt{3}}{12}, \qquad t_1 = t + \left(\frac{1}{2} - \frac{\sqrt{3}}{6}\right)\Delta t, \quad (67) $$

$$ \alpha_2 = \frac{3 + 2\sqrt{3}}{12}, \qquad t_2 = t + \left(\frac{1}{2} + \frac{\sqrt{3}}{6}\right)\Delta t. \quad (68) $$

$\hat{H}[t_1]$ and $\hat{H}[t_2]$ are the Hamiltonians at times $t_1$ and $t_2$, which are, in fact, unknown because they depend on the Kohn–Sham states through the density: we are dealing with a non-linear problem and the CFM expansions were, in fact, developed for linear systems. We have various options to extend them for our non-linear problem: for example, one could define $\hat{H}[t_1]$ and $\hat{H}[t_2]$ as interpolated Hamiltonians from $\hat{H}[t]$ and $\hat{H}[t + \Delta t]$ in which case we end up with an implicit equation for $\varphi(t + \Delta t)$ that we would have to solve at a substantial cost.

The alternative that we implemented, however, is to approximate $\hat{H}[t_1]$ via an extrapolation from the Hamiltonian at various previous time steps (in practice, it is the Hartree, exchange, and correlation parts that must be extrapolated). The resulting method is then explicit, i.e., no linear or non-linear algebraic equations need to be solved. The fourth order accuracy is preserved as long as the extrapolation is also done at order four.

As an example, Fig. 14 shows the performance of this CFM4 method against the well-known exponential midpoint rule (EMR). The benchmark system consists of a benzene molecule. It is subject to an instantaneous perturbation at the beginning of the propagation, and then, it evolves freely for some fixed interval of time. The propagations are performed at varying values of $\Delta t$. The plot displays the cost of the propagation vs the accuracy achieved in each

<figure><img src="image_26.png" /><figcaption>FIG. 14. Wall-time computational cost of the method (in seconds) against the error in the wave function [defined as $Error(T, \Delta t) = \sqrt{\sum_m ||\varphi_m(T) - \varphi_m^{exact}(T)||^2}$, where the $\varphi_m^{exact}$ are the values of the Kohn–Sham orbitals computed using the explicit fourth-order Runge–Kutta integrator with a very small time step] for the EMR and CFM4 propagators.</figcaption></figure>

case. As we can see, the CFM4 outperforms the EMR for all values of the error examined, although the differences become more obvious when higher accuracy is demanded. The key difference between the two methods is the fourth-order accuracy of the CFM4 scheme—at the cost of requiring two exponentials, instead of just one, while the EMR is only second-order accurate. This is reflected in the different slopes of the curves as the error becomes smaller, i.e., as $\Delta t \to 0$.

## XIV. CONJUGATE GRADIENT IMPLEMENTATION IN RDMFT

The RDMFT implementation in Octopus has been described in detail in a previous paper.¹⁹ In this section, we review briefly the existing implementation, providing some new insights, and introduce a recently implemented method to solve the RDMFT equations, which is better suited to the real-space grids used in the code.

The optimization of the natural orbitals in RDMFT is subject to an orthonormalization constraint for the orbitals. One minimizes the functional

$$ E_{\text{total}} - \sum_{j,k=1}^{M} \lambda_{jk} \left( \int d\mathbf{r} \psi_j^*(\mathbf{r}) \psi_k(\mathbf{r}) - \delta_{jk} \right), \quad (69) $$

where $E_{\text{total}}$ and $M$ denote the total energy and the number of natural orbitals $\psi_j$ included in the calculation, respectively. The parameters $\lambda_{jk}$ are the Lagrange multipliers, which ensure the orthonormality of the natural orbitals at the solution point. We note that, for an RDMFT calculation, the number of natural orbitals has to be larger than the number of electrons in the system (core electrons that are included in the pseudopotential do not count). The exact number $M$ is system-dependent and should be treated as an additional parameter with respect to which a convergence study should be carried out, just like it is done for the basis set or the grid parameters. Typically, the optimization is performed using the so-called Piris method,¹⁶⁷ which was the method previously implemented in Octopus. Within this approach, one uses the orthonormality constraint of the natural orbitals, which implies that a certain matrix constructed from the Lagrange multipliers $\lambda_{jk}$ is diagonal at the solution point. As an immediate consequence of the Piris method, the natural orbitals at the solution point are linear combinations of the orbitals used as a starting point for the minimization. In other words, the initial orbitals serve as a basis. Consequently, the necessary matrix elements for different energy contributions can be calculated for the basis functions (initial orbitals) before the iterative optimization of the natural orbitals and their occupation numbers is started. In addition, the optimization of the occupation numbers can be turned off completely, resulting in a Hartree–Fock calculation on the basis of the initial orbitals.

For the existing implementation in Octopus of the Piris method, the initial orbitals are taken to be the solutions obtained with a different level of theory, such as independent particles or DFT. In order to better understand the effect of the choice of basis, we tested the following choices: (i) independent particles, density functional theory within (ii) the local density approximation (LDA) or (iii) the exact exchange (EXX) approximation, and (iv) the Hartree–Fock approximation. In all cases, we have to ensure that the number
of unoccupied states in the calculation is sufficient to cover all the natural orbitals that will obtain significant occupation in the following RDMFT calculation. The results for the convergence of the total energy of a one-dimensional (1D) hydrogen molecule using the Müller functional¹⁶⁸ are given in Fig. 15. The calculations were performed on a 1D grid extending from –12.0 bohr to 12.0 bohrs, with a grid spacing of 0.03 bohr. The nuclear potential for the 1D hydrogen molecule reads

$$v(x) = - \frac{1}{\sqrt{(x-d)^2 + 1}} - \frac{1}{\sqrt{(x+d)^2 + 1}}, \quad (70)$$

with $d = 1.628$ bohr, which corresponds to the equilibrium geometry. The electron-electron interaction in one dimension is described by the soft-Coulomb interaction

$$w(x, x') = \frac{1}{\sqrt{(x-x')^2 + 1}}. \quad (71)$$

Since Octopus performs all calculations on a finite grid, we typically obtain a finite number of bound states in the calculations for the basis set and any additional orbitals that extend over the whole grid, i.e., they are unbound and, therefore, delocalized. However, all natural orbitals with non-zero occupation, because they decay with the ionization potential of the system,¹⁷⁰ are localized on the system. Hence, the extended basis states will only contribute with very small coefficients, if at all, and their inclusion in the basis set does not lead to a significant improvement of the results. In the past, this problem was addressed by performing an additional step to localize the initial states before starting the RDMFT calculation. However, further investigations showed that this only improves the results for a small number of natural orbitals in the calculation. Testing the convergence with respect to the number of natural orbitals then showed that the additional localization step slows down the convergence

<figure><img src="image_25.png" /><figcaption>FIG. 15. Total energy for the RDMFT calculation for one-dimensional H₂¹⁶⁹ using the Piris method with different basis sets and using the conjugate gradient implementation. The inset shows a zoomed-in view of the area where convergence is reached. We employ the Müller functional¹⁶⁸ for all calculations.</figcaption></figure>

with respect to the number of basis functions. The fastest convergence and lowest total energies are obtained by using the results from an independent particle calculation, as those yield the largest number of localized orbitals from all the different basis sets that were tested, as shown in Fig. 15. As the additional localization step proved to be unnecessary and even hinders convergence, it has been removed from the new version of the code.

Since the natural way of representing quantities in Octopus is directly on the real-space grid, and to circumvent the limitations with the quality of the basis sets available for the Piris method, we decided to implement a conjugate gradient optimization of the natural orbitals. This implementation follows the procedure for DFT explained in Ref. 171, which we adapted for RDMFT. This procedure allows us to take advantage of the full flexibility of a real-space grid and provides a systematic way of improving the results by enlarging the grid and reducing the spacing between grid points. The conjugate gradient algorithm requires a set of initial orbitals to start the self-consistent calculation; however, at convergence, the results are independent of that starting point. Therefore, while the calculation using the Piris method requires a set of initial states that serve as the basis, the conjugate gradient algorithm can be used starting from an initial set of random states. In our tests of the conjugate gradient implementation, the quality of the initial states only had an influence on the number of iterations necessary for the convergence, but not on the final result. We suggest to use the orbitals obtained from an independent particle calculation as initial states since they can be obtained for a small numerical cost and simultaneously can serve as a basis set in the Piris implementation.

Since we are not solving an eigenvalue equation to obtain the natural orbitals, they are not automatically orthogonal. As mentioned above, the orthogonality is taken into account via a constraint. Compared to Ref. 171, the non-diagonal Lagrange multipliers $\lambda_{ik}$ lead to two modifications in the conjugate gradient procedure. First, the steepest-descent direction [cf. Eq. (5.10) of Ref. 171] reads here

$$\zeta_i^m = -\hat{H}\psi_i^m + \sum_{k=1}^{M} \lambda_{ik}\psi_k^m, \quad (72)$$

where

$$\lambda_{ik} = \langle \psi_i^m | \hat{H} \psi_k^m \rangle. \quad (73)$$

Second, the one parameter of the line-minimization [cf. Eq. (5.26) of Ref. 171] is changed to

$$\frac{\partial E_{\text{total}}}{\partial \Theta} = \langle \phi_i^{m'} | \hat{H} \psi_i^m \rangle + \langle \psi_i^m | \hat{H} \phi_i^{m'} \rangle \\ - \sum_k (\lambda_{ki} \langle \phi_i^{m'} | \psi_i^m \rangle + \lambda_{ik} \langle \psi_i^m | \phi_i^{m'} \rangle), \quad (74)$$

where $\Theta$ parameterizes the descend in the direction of the gradient and the single-particle Hamiltonian acting on a state $|\phi\rangle$ is defined as

$$\hat{H}|\phi\rangle = \frac{\partial E_{\text{total}}}{\partial \phi^*}. \quad (75)$$

For details of the notation, we refer the reader to Ref. 171.

The convergence study with respect to the number of natural orbitals for the conjugate gradient algorithm is also included in
of unoccupied states in the calculation is sufficient to cover all the natural orbitals that will obtain significant occupation in the following RDMFT calculation. The results for the convergence of the total energy of a one-dimensional (1D) hydrogen molecule using the Müller functional¹⁶⁸ are given in Fig. 15. The calculations were performed on a 1D grid extending from -12.0 bohr to 12.0 bohrs, with a grid spacing of 0.03 bohr. The nuclear potential for the 1D hydrogen molecule reads

$$v(x) = - \frac{1}{\sqrt{(x-d)^2 + 1}} - \frac{1}{\sqrt{(x+d)^2 + 1}}, \quad (70)$$

with d = 1.628 bohr, which corresponds to the equilibrium geometry. The electron-electron interaction in one dimension is described by the soft-Coulomb interaction

$$w(x, x') = \frac{1}{\sqrt{(x-x')^2 + 1}}. \quad (71)$$

Since Octopus performs all calculations on a finite grid, we typically obtain a finite number of bound states in the calculations for the basis set and any additional orbitals that extend over the whole grid, i.e., they are unbound and, therefore, delocalized. However, all natural orbitals with non-zero occupation, because they decay with the ionization potential of the system,¹⁷⁰ are localized on the system. Hence, the extended basis states will only contribute with very small coefficients, if at all, and their inclusion in the basis set does not lead to a significant improvement of the results. In the past, this problem was addressed by performing an additional step to localize the initial states before starting the RDMFT calculation. However, further investigations showed that this only improves the results for a small number of natural orbitals in the calculation. Testing the convergence with respect to the number of natural orbitals then showed that the additional localization step slows down the convergence

<figure><img src="image_27.png" /><figcaption>FIG. 15. Total energy for the RDMFT calculation for one-dimensional H₂¹⁶⁹ using the Piris method with different basis sets and using the conjugate gradient implementation. The inset shows a zoomed-in view of the area where convergence is reached. We employ the Müller functional¹⁶⁸ for all calculations.</figcaption></figure>

with respect to the number of basis functions. The fastest convergence and lowest total energies are obtained by using the results from an independent particle calculation, as those yield the largest number of localized orbitals from all the different basis sets that were tested, as shown in Fig. 15. As the additional localization step proved to be unnecessary and even hinders convergence, it has been removed from the new version of the code.

Since the natural way of representing quantities in Octopus is directly on the real-space grid, and to circumvent the limitations with the quality of the basis sets available for the Piris method, we decided to implement a conjugate gradient optimization of the natural orbitals. This implementation follows the procedure for DFT explained in Ref. 171, which we adapted for RDMFT. This procedure allows us to take advantage of the full flexibility of a real-space grid and provides a systematic way of improving the results by enlarging the grid and reducing the spacing between grid points. The conjugate gradient algorithm requires a set of initial orbitals to start the self-consistent calculation; however, at convergence, the results are independent of that starting point. Therefore, while the calculation using the Piris method requires a set of initial states that serve as the basis, the conjugate gradient algorithm can be used starting from an initial set of random states. In our tests of the conjugate gradient implementation, the quality of the initial states only had an influence on the number of iterations necessary for the convergence, but not on the final result. We suggest to use the orbitals obtained from an independent particle calculation as initial states since they can be obtained for a small numerical cost and simultaneously can serve as a basis set in the Piris implementation.

Since we are not solving an eigenvalue equation to obtain the natural orbitals, they are not automatically orthogonal. As mentioned above, the orthogonality is taken into account via a constraint. Compared to Ref. 171, the non-diagonal Lagrange multipliers $\lambda_{ik}$ lead to two modifications in the conjugate gradient procedure. First, the steepest-descent direction [cf. Eq. (5.10) of Ref. 171] reads here

$$\zeta_i^m = -\hat{H}\psi_i^m + \sum_{k=1}^{M} \lambda_{ik}\psi_k^m, \quad (72)$$

where

$$\lambda_{ik} = \langle \psi_i^m | \hat{H} \psi_k^m \rangle. \quad (73)$$

Second, the one parameter of the line-minimization [cf. Eq. (5.26) of Ref. 171] is changed to

$$\frac{\partial E_{\text{total}}}{\partial \Theta} = \langle \phi_i^{m'} | \hat{H} \psi_i^m \rangle + \langle \psi_i^m | \hat{H} \phi_i^{m'} \rangle \\ - \sum_k (\lambda_{ki} \langle \phi_i^{m'} | \psi_i^m \rangle + \lambda_{ik} \langle \psi_i^m | \phi_i^{m'} \rangle), \quad (74)$$

where $\Theta$ parameterizes the descend in the direction of the gradient and the single-particle Hamiltonian acting on a state $|\phi\rangle$ is defined as

$$\hat{H}|\phi\rangle = \frac{\partial E_{\text{total}}}{\partial \phi^*}. \quad (75)$$

For details of the notation, we refer the reader to Ref. 171.

The convergence study with respect to the number of natural orbitals for the conjugate gradient algorithm is also included in
Fig. 15. As one can see, a smaller number of natural orbitals needs to be included in the calculation than for the basis set implemen- tations with the Piris method. As the number of natural orbitals equals the number of basis functions in these calculations, this is mostly due to the fact that the available basis sets are of rather poor quality. In addition, the converged total energy for the con- jugate gradient algorithm is slightly lower than for all the basis sets using the Piris method (see the inset of Fig. 15), which shows that the conjugate gradient algorithm exploits the full flexibility of the grid implementation, allowing for contributions to the nat- ural orbitals that are not covered by any of the basis sets. We have also verified that the converged result of the conjugate gradi- ent method is indeed independent of the choice of initial state for the Müller functional¹⁶⁸ employed in our calculations. The Müller functional is known to be convex when all infinitely many natu- ral orbitals are included.¹⁷² This property is most likely not shared by all available RDMFT functionals. In addition, the number of natural orbitals in any practical calculation is always finite. Conse- quently, in practice, one needs to test the convergence for the dif- ferent starting points, as the appearance of local minima cannot be excluded.

XV. PERIODIC SYSTEMS AND SYMMETRIES

Electronic structure in periodic systems is usually described using plane waves, but real-space grids have been shown to be a viable alternative when performing DFT and TDDFT calcula- tions.²,¹⁷³ Unfortunately, the discretization introduced by the real- space grid often breaks the direct connection between the physical system and the basis set, as symmetries and translations are, in most of the cases, incompatible with the discretized grid. However, real- space grids offer many advantages, such as natural mixed periodic- boundary conditions for semi-periodic systems, and the calculation of the exchange and correlation term of DFT is straightforwardly obtained on the grid.

One of the main challenges for treating periodic systems in real
space is the generation of the real-space grid and the corresponding
weights for the finite differences. This becomes relevant when the
grid is generated along the primitive axes of the Wigner-Seitz cell
(primitive cell) of a solid, where the generating axes are usually non-
orthogonal.

In Octopus, the grid points are generated along the primitive
axes, and the calculation of the finite-difference weights follows the
implementation described in Ref. 174. Once the grid is generated,
and the weights for the finite differences (gradient and Laplacian) are
obtained, it is necessary to deal with the discretization of the recip-
rocral space. Indeed, for periodic systems, the full crystal is replaced
by the primitive cell of the crystal in real-space, thanks to the Bloch
theorem, but is then complemented by the Brillouin zone, which also
needs to be sampled. This grid, usually called k-point grid, is com-
mon to any type of basis sets, as long as one decides to reduce the
crystal to its primitive cell.

In order to reduce the numerical effort associated with the
description of periodic systems, we make use of the symmetries to
reduce the Brillouin zone to its irreducible zone, which can drasti-
cally reduce the number of k-points. The space group of the crys-
tal, as well as its symmetries, are obtained, thanks to the spglib
library.¹⁷⁵ In the present implementation, we restrict the symmetries

to the symmorphic symmetries (inversion, rotations, and mirror
planes), leaving for later the so-called non-symmorphic symmetries,
i.e., symmetries involving a fractional translation, as they are not
compatible with arbitrary real-space grids.

In order to assess the validity of our implementation, we used the so-called delta-factor test.¹⁷⁶ Using the Schlipf-Gygi ONCVPSP 2015 pseudopotential set, we obtained the value of 1.50 meV/atom, which is very close to the value obtained by other codes using the same pseudopotential set.

When investigating the electron dynamics driven by a laser
field or any type of symmetry-breaking perturbation (vector poten-
tials, kicks with a finite momentum, strain, etc.), some of the origi-
nal symmetries are lost. To deal with that, Octopus finds the small
group of symmetries corresponding to the original space-group of
the solid, retaining only the symmetries that leave invariant the per-
turbation direction. This defines the symmetries that are used for a
time-dependent calculation.

One important aspect when using symmetries is the sym-
metrization (in real-space) of the charge and current densities, as
well as other observables, such as kinetic energy density. In com-
plement to the reduction of the Brillouin zone to its irreducible
zone, we implemented a real-space symmetrization. We found that
performing the real-space symmetrization is very important to
achieve good and stable numerical results when taking into account
symmetries.

For a comprehensive description of periodic systems, the ion
dynamics has to be considered in addition to the electron dynam-
ics. For isolated systems, it is often described with the TDDFT-
Ehrenfest dynamics method,¹⁷,¹⁷⁷ where the ions obey the New-
ton equation with force fields computed by TDDFT. In contrast,
for periodic systems, the ion dynamics has to be described with
two sets of equations: one is the Newton equation for ions on
the reduced coordinates in the primitive cell and the other is the
equation of motion for the primitive cell itself, as the ionic coor-
dinates of periodic systems are described by the combination of
the reduced coordinates with the lattice vectors of the primitive
cell. The lattice dynamics is often treated with the Parrinello-
Rahman method,¹⁷⁸,¹⁷⁹ and the key ingredient of the equation of
motion is the stress tensor. Therefore, to realize the *ab initio* sim-
ulations for electron-ion-lattice dynamics based on TDDFT, the
calculation of the stress-tensor has been implemented based on
Ref. 180, and the implementation of the lattice dynamics is under
way.

XVI. ADDITIONAL TECHNICAL CODE IMPROVEMENTS

In order to make all the new developments and applications
described in the above sections possible, the code needs to be accu-
rate, efficient, and reliable. In addition, when a code reaches the
size and complexity of Octopus (currently more than 200 000 lines
of source code), the amount of time required for maintenance and
for adding new features becomes considerable. Therefore, continu-
ous efforts at optimizing, validating, and improving the source code
quality are needed. These efforts are essential, but often do not get
the attention they deserve. In this section, we present some note-
worthy developments to improve the code reliability and efficiency,
to make the project easier to maintain, and to tackle new computer
architecture challenges.
Fig. 15. As one can see, a smaller number of natural orbitals needs to be included in the calculation than for the basis set implemen- tations with the Piris method. As the number of natural orbitals equals the number of basis functions in these calculations, this is mostly due to the fact that the available basis sets are of rather poor quality. In addition, the converged total energy for the con- jugate gradient algorithm is slightly lower than for all the basis sets using the Piris method (see the inset of Fig. 15), which shows that the conjugate gradient algorithm exploits the full flexibility of the grid implementation, allowing for contributions to the nat- ural orbitals that are not covered by any of the basis sets. We have also verified that the converged result of the conjugate gradi- ent method is indeed independent of the choice of initial state for the Müller functional¹⁶⁸ employed in our calculations. The Müller functional is known to be convex when all infinitely many natu- ral orbitals are included.¹⁷² This property is most likely not shared by all available RDMFT functionals. In addition, the number of natural orbitals in any practical calculation is always finite. Conse- quently, in practice, one needs to test the convergence for the dif- ferent starting points, as the appearance of local minima cannot be excluded.

XV. PERIODIC SYSTEMS AND SYMMETRIES

Electronic structure in periodic systems is usually described using plane waves, but real-space grids have been shown to be a viable alternative when performing DFT and TDDFT calcula- tions.<sup>2,173</sup> Unfortunately, the discretization introduced by the real- space grid often breaks the direct connection between the physical system and the basis set, as symmetries and translations are, in most of the cases, incompatible with the discretized grid. However, real- space grids offer many advantages, such as natural mixed periodic- boundary conditions for semi-periodic systems, and the calculation of the exchange and correlation term of DFT is straightforwardly obtained on the grid.

One of the main challenges for treating periodic systems in real
space is the generation of the real-space grid and the corresponding
weights for the finite differences. This becomes relevant when the
grid is generated along the primitive axes of the Wigner-Seitz cell
(primitive cell) of a solid, where the generating axes are usually non-
orthogonal.

In Octopus, the grid points are generated along the primitive
axes, and the calculation of the finite-difference weights follows the
implementation described in Ref. 174. Once the grid is generated,
and the weights for the finite differences (gradient and Laplacian) are
obtained, it is necessary to deal with the discretization of the recip-
rocral space. Indeed, for periodic systems, the full crystal is replaced
by the primitive cell of the crystal in real-space, thanks to the Bloch
theorem, but is then complemented by the Brillouin zone, which also
needs to be sampled. This grid, usually called k-point grid, is com-
mon to any type of basis sets, as long as one decides to reduce the
crystal to its primitive cell.

In order to reduce the numerical effort associated with the
description of periodic systems, we make use of the symmetries to
reduce the Brillouin zone to its irreducible zone, which can drasti-
cally reduce the number of k-points. The space group of the crys-
tal, as well as its symmetries, are obtained, thanks to the spglib
library.¹⁷⁵ In the present implementation, we restrict the symmetries

to the symmorphic symmetries (inversion, rotations, and mirror
planes), leaving for later the so-called non-symmorphic symmetries,
i.e., symmetries involving a fractional translation, as they are not
compatible with arbitrary real-space grids.

In order to assess the validity of our implementation, we used the so-called delta-factor test.<sup>176</sup> Using the Schlipf-Gygi ONCVPSP 2015 pseudopotential set, we obtained the value of 1.50 meV/atom, which is very close to the value obtained by other codes using the same pseudopotential set.

When investigating the electron dynamics driven by a laser
field or any type of symmetry-breaking perturbation (vector poten-
tials, kicks with a finite momentum, strain, etc.), some of the origi-
nal symmetries are lost. To deal with that, Octopus finds the small
group of symmetries corresponding to the original space-group of
the solid, retaining only the symmetries that leave invariant the per-
turbation direction. This defines the symmetries that are used for a
time-dependent calculation.

One important aspect when using symmetries is the sym-
metrization (in real-space) of the charge and current densities, as
well as other observables, such as kinetic energy density. In com-
plement to the reduction of the Brillouin zone to its irreducible
zone, we implemented a real-space symmetrization. We found that
performing the real-space symmetrization is very important to
achieve good and stable numerical results when taking into account
symmetries.

For a comprehensive description of periodic systems, the ion
dynamics has to be considered in addition to the electron dynam-
ics. For isolated systems, it is often described with the TDDFT-
Ehrenfest dynamics method,<sup>17,177</sup> where the ions obey the New-
ton equation with force fields computed by TDDFT. In contrast,
for periodic systems, the ion dynamics has to be described with
two sets of equations: one is the Newton equation for ions on
the reduced coordinates in the primitive cell and the other is the
equation of motion for the primitive cell itself, as the ionic coor-
dinates of periodic systems are described by the combination of
the reduced coordinates with the lattice vectors of the primitive
cell. The lattice dynamics is often treated with the Parrinello-
Rahman method,<sup>17,179</sup> and the key ingredient of the equation of
motion is the stress tensor. Therefore, to realize the *ab initio* sim-
ulations for electron-ion-lattice dynamics based on TDDFT, the
calculation of the stress-tensor has been implemented based on
Ref. 180, and the implementation of the lattice dynamics is under
way.

XVI. ADDITIONAL TECHNICAL CODE IMPROVEMENTS

In order to make all the new developments and applications
described in the above sections possible, the code needs to be accu-
rate, efficient, and reliable. In addition, when a code reaches the
size and complexity of Octopus (currently more than 200 000 lines
of source code), the amount of time required for maintenance and
for adding new features becomes considerable. Therefore, continu-
ous efforts at optimizing, validating, and improving the source code
quality are needed. These efforts are essential, but often do not get
the attention they deserve. In this section, we present some note-
worthy developments to improve the code reliability and efficiency,
to make the project easier to maintain, and to tackle new computer
architecture challenges.
## A. Web application for analyzing regression tests

In a complex code such as Octopus, every new development may have unintended side effects, possibly introducing errors to already existing features. To avoid such regressions, every change to the code is required to pass a suite of tests before being accepted. This suite is executed automatically using Buildbot¹⁸¹ and is integrated with the continuous integration framework of the gitlab platform where Octopus is hosted for several years now. The testsuite is executed with 30 different toolchains spanning a variety of architectures, compilers, and Message-Passing Interface (MPI) libraries. It contains currently about 180 tests that execute Octopus roughly 740 times and that contain about 12 000 comparisons to reference values. The overall coverage of the tests is determined using the codecov.io service and lies at about 71% at the moment.

Although the testsuite is efficient in avoiding regressions, analyzing why a test failed has been difficult so far due to the large amount of data generated. Thus, we implemented an interactive analysis and visualization of the testsuite results as a web application available at https://octopus-code.org/testsuite/. This application makes it easier to understand why a test failed: Is it a problem with just one toolchain, e.g., a particular architecture, or is it simply a larger numerical variation of the results? Moreover, it facilitates updating the tests and improving the testsuite itself.

The application was implemented using the web framework Django coupled to a Postgres database. The testsuite results are automatically uploaded by the Buildbot service as soon as they are available. The application allows us to analyze single toolchain runs and single comparison matches and to compare all toolchain runs for a single commit in git. It provides histograms to judge if there are outliers or if it is a broad distribution; moreover, this allows us to determine, e.g., if there is a difference between MPI and non-MPI toolchains. This application has already helped us to identify bugs causing regressions and will continue to be useful to understand failed tests, to update them, and to improve the testsuite.

## B. Improving ground-state calculations

To improve the reliability of ground-state calculations with Octopus, the default eigensolver used inside the self-consistent field (SCF) cycle was improved and different real-space preconditioners for the eigensolvers were evaluated.

The default eigensolver, a conjugate-gradient algorithm, has been improved over the previous implementation by now following closely Ref. 171, which greatly improves the reliability of the solver. The updated implementation differs from Ref. 171 in only one point: by default, the current band is not orthogonalized against all bands, instead it is only orthogonalized against previously computed bands with lower energies. According to our tests, in Octopus, this leads to a faster convergence for most cases.

Preconditioners for eigensolvers are an integral part of SCF calculations because they greatly accelerate convergence. They achieve this by applying an approximate inverse of the Hamiltonian in each iteration that leaves the solution invariant and brings the system closer to the real solution. As Octopus uses a real-space grid (as opposed to plane waves like in Ref. 171), different preconditioners are needed. A range of preconditioners has been compared for a wide variety of systems, comprising molecules, semiconductors,

metallic systems, and surfaces to test convergence for very different cases.

The default preconditioner in Octopus is a low-pass filter obtained by adding the neighboring values in each dimension to the current value at a grid point, weighted by a certain factor $\alpha$,

$$ \psi'_{i,j,k} = \alpha \psi_{i,j,k} + (1 - \alpha) (\psi_{i-1,j,k} + \psi_{i+1,j,k} + \psi_{i,j-1,k} \\ + \psi_{i,j+1,k} + \psi_{i,j,k-1} + \psi_{i,j,k+1}). \quad (76) $$

This preconditioner was first described by Saad et al.¹⁸² with $\alpha = 0.5$. It was proven to be the most effective preconditioner in our comparison because it is quite cheap to apply and it, nevertheless, decreases the number of SCF iterations noticeably. It can also be understood as two weighted Jacobi iterations (up to a prefactor) to solve for the inverse of the kinetic term of the Hamiltonian (i.e., $0.5\Delta\psi' = \psi$). From the convergence radius of the Jacobi iterations, we can conclude that the allowed values for $\alpha$ are between 0.5 and 1. From theoretical considerations of damping of different spatial wavelengths during Jacobi iterations (see Secs. 4.1.3 and 4.2 of Ref. 183), the ideal value should be 0.75 to most effectively damp high spatial frequencies. In practice, we find that this depends on the system. Moreover, increasing the number of Jacobi iterations does not make the preconditioner more effective. In addition, using a single Jacobi iteration (i.e., dividing by the diagonal of the Laplacian) does not speed up the convergence significantly.

A multigrid method similar to the one implemented in GPAW¹⁸⁴ also uses Jacobi iterations to solve for the kinetic Hamiltonian but employs different grids to reach a faster decrease in the error. Although this method reduces the number of SCF iterations, it is computationally more expensive because the Laplacian is applied several times for each iteration and is, thus, less effective than the filter preconditioner.

A preconditioner specifically targeting real-space methods was proposed by Seitsonen et al.¹⁸⁵ It uses the ratio between the difference of the energy and the potential to the kinetic energy [their Eq. (3)] with a preconditioning function [their Eq. (4)] originally from Ref. 171. However, this preconditioner did not speed up convergence significantly when used in Octopus.

Another way to get an approximate inverse of the Hamiltonian is to solve the Poisson equation associated with the kinetic part of the Hamiltonian. This reduces the number of iterations needed to converge the SCF calculation, but it increases the total calculation time, as each iteration is much more costly.

In summary, we find that for our real-space code, the filter preconditioner is most effective in reducing the total time needed to compute ground states because it reduces the total number of iterations without being computationally too expensive.

Nonetheless, sometimes certain states cannot be fully converged using a preconditioner, especially for calculations with many unoccupied states. In this case, restarting without preconditioner is needed to obtain a complete SCF convergence. The reasons for this are still under investigation.

## C. Novel multi-system framework

In the way Octopus was originally designed, the entire code was structured around the idea that there was only one system (albeit of arbitrary size and complexity) with an associated Hamiltonian.
## A. Web application for analyzing regression tests

In a complex code such as Octopus, every new development may have unintended side effects, possibly introducing errors to already existing features. To avoid such regressions, every change to the code is required to pass a suite of tests before being accepted. This suite is executed automatically using Buildbot¹⁸¹ and is integrated with the continuous integration framework of the gitlab platform where Octopus is hosted for several years now. The testsuite is executed with 30 different toolchains spanning a variety of architectures, compilers, and Message-Passing Interface (MPI) libraries. It contains currently about 180 tests that execute Octopus roughly 740 times and that contain about 12 000 comparisons to reference values. The overall coverage of the tests is determined using the codecov.io service and lies at about 71% at the moment.

Although the testsuite is efficient in avoiding regressions, analyzing why a test failed has been difficult so far due to the large amount of data generated. Thus, we implemented an interactive analysis and visualization of the testsuite results as a web application available at https://octopus-code.org/testsuite/. This application makes it easier to understand why a test failed: Is it a problem with just one toolchain, e.g., a particular architecture, or is it simply a larger numerical variation of the results? Moreover, it facilitates updating the tests and improving the testsuite itself.

The application was implemented using the web framework Django coupled to a Postgres database. The testsuite results are automatically uploaded by the Buildbot service as soon as they are available. The application allows us to analyze single toolchain runs and single comparison matches and to compare all toolchain runs for a single commit in git. It provides histograms to judge if there are outliers or if it is a broad distribution; moreover, this allows us to determine, e.g., if there is a difference between MPI and non-MPI toolchains. This application has already helped us to identify bugs causing regressions and will continue to be useful to understand failed tests, to update them, and to improve the testsuite.

## B. Improving ground-state calculations

To improve the reliability of ground-state calculations with Octopus, the default eigensolver used inside the self-consistent field (SCF) cycle was improved and different real-space preconditioners for the eigensolvers were evaluated.

The default eigensolver, a conjugate-gradient algorithm, has been improved over the previous implementation by now following closely Ref. 171, which greatly improves the reliability of the solver. The updated implementation differs from Ref. 171 in only one point: by default, the current band is not orthogonalized against all bands, instead it is only orthogonalized against previously computed bands with lower energies. According to our tests, in Octopus, this leads to a faster convergence for most cases.

Preconditioners for eigensolvers are an integral part of SCF calculations because they greatly accelerate convergence. They achieve this by applying an approximate inverse of the Hamiltonian in each iteration that leaves the solution invariant and brings the system closer to the real solution. As Octopus uses a real-space grid (as opposed to plane waves like in Ref. 171), different preconditioners are needed. A range of preconditioners has been compared for a wide variety of systems, comprising molecules, semiconductors,

metallic systems, and surfaces to test convergence for very different cases.

The default preconditioner in Octopus is a low-pass filter obtained by adding the neighboring values in each dimension to the current value at a grid point, weighted by a certain factor $\alpha$,

$$ \psi'_{i,j,k} = \alpha \psi_{i,j,k} + (1 - \alpha) (\psi_{i-1,j,k} + \psi_{i+1,j,k} + \psi_{i,j-1,k} \\ + \psi_{i,j+1,k} + \psi_{i,j,k-1} + \psi_{i,j,k+1}). \quad (76) $$

This preconditioner was first described by Saad et al.¹⁸² with $\alpha = 0.5$. It was proven to be the most effective preconditioner in our comparison because it is quite cheap to apply and it, nevertheless, decreases the number of SCF iterations noticeably. It can also be understood as two weighted Jacobi iterations (up to a prefactor) to solve for the inverse of the kinetic term of the Hamiltonian (i.e., $0.5\Delta\psi' = \psi$). From the convergence radius of the Jacobi iterations, we can conclude that the allowed values for $\alpha$ are between 0.5 and 1. From theoretical considerations of damping of different spatial wavelengths during Jacobi iterations (see Secs. 4.1.3 and 4.2 of Ref. 183), the ideal value should be 0.75 to most effectively damp high spatial frequencies. In practice, we find that this depends on the system. Moreover, increasing the number of Jacobi iterations does not make the preconditioner more effective. In addition, using a single Jacobi iteration (i.e., dividing by the diagonal of the Laplacian) does not speed up the convergence significantly.

A multigrid method similar to the one implemented in GPAW¹⁸⁴ also uses Jacobi iterations to solve for the kinetic Hamiltonian but employs different grids to reach a faster decrease in the error. Although this method reduces the number of SCF iterations, it is computationally more expensive because the Laplacian is applied several times for each iteration and is, thus, less effective than the filter preconditioner.

A preconditioner specifically targeting real-space methods was proposed by Seitsonen et al.¹⁸⁵ It uses the ratio between the difference of the energy and the potential to the kinetic energy [their Eq. (3)] with a preconditioning function [their Eq. (4)] originally from Ref. 171. However, this preconditioner did not speed up convergence significantly when used in Octopus.

Another way to get an approximate inverse of the Hamiltonian is to solve the Poisson equation associated with the kinetic part of the Hamiltonian. This reduces the number of iterations needed to converge the SCF calculation, but it increases the total calculation time, as each iteration is much more costly.

In summary, we find that for our real-space code, the filter preconditioner is most effective in reducing the total time needed to compute ground states because it reduces the total number of iterations without being computationally too expensive.

Nonetheless, sometimes certain states cannot be fully converged using a preconditioner, especially for calculations with many unoccupied states. In this case, restarting without preconditioner is needed to obtain a complete SCF convergence. The reasons for this are still under investigation.

## C. Novel multi-system framework

In the way Octopus was originally designed, the entire code was structured around the idea that there was only one system (albeit of arbitrary size and complexity) with an associated Hamiltonian.
This was typically a combined system of electrons and nuclei, where the former were described using DFT and the latter were treated classically. All possible sorts of interactions of this system with some external source (e.g., external electromagnetic fields, solvents described with the PCM, etc.) were then added in an *ad hoc* fashion. Although this approach worked well for many features and applications, its limitations became evident with several recent developments, in particular, the coupled Maxwell-Kohn-Sham equations described in Sec. II, where the Kohn-Sham orbitals need to be time-propagated alongside the Maxwell fields.

This has prompted us to start a major refactoring of the code, where the full physical system is treated as several subsystems interacting with each other. Such subsystems can be electrons, nuclei, and Maxwell fields. While this framework was mainly motivated by the Maxwell-TDDFT coupling, it has been implemented in a very general fashion so that all existing features can be converted to it and that, in the future, many other developments that require the coupling of several subsystems can be based on it.

## D. Memory layout

A new memory layout for storing the orbitals was introduced in Ref. 186, where all states are stored in a number of smaller batches with the innermost index being the state index instead of the real-space grid index. Thus, for all states in a batch, exactly the same operations can be executed when looping over the grid while accessing memory contiguously. This allows efficient parallelization on graphical processing units (GPUs), where all threads in a warp need to execute the same instructions, as well as vectorization on central processing units (CPUs), where one instruction can operate on several data points. To fully utilize these instructions, the kernel for computing finite differences (e.g., the Laplacian) has been specialized to explicitly use SSE, AVX, or AVX512 instructions. Now, more parts of the code have been ported to use this new layout to increase their performance, and whenever possible, the code will use this layout by default.

## E. GPUs

The first GPU implementation was described in Ref. 186. This was based on OpenCL and was limited to a small selection of numerical algorithms and calculation modes. These included some of the most commonly used features of the code or algorithms that were particularly suited for GPU porting, such as the residual minimization method-direct inversion in the iterative subspace (RMM-DIIS) eigensolver for ground-state calculations or the enforced time-reversal-symmetry propagator for time-dependent calculations. Since then, the GPU implementation has

been expanded in several different ways. Most notably, it now supports CUDA through an additional compatibility layer. The number of supported features and algorithms has also been increasing steadily such that most time-dependent calculations can now be run efficiently on GPUs. The implementation was also expanded to support multiple devices per host. Using the packed storage format described above, Octopus is able to store the states fully in the GPU memory, provided the memory is large enough, thus reducing memory transfers to a minimum. Recently, this was improved further by removing frequent allocations and deallocations of temporary variables on the GPU, now using a custom memory management for those. This proved very effective in improving the scaling to several GPUs and also to several nodes with GPUs.

We show two examples of time-dependent runs in Fig. 16. For the first example (left panel), β-cyclodextrine was used as an input system and the simulation was run on the GPU island of the COBRA supercomputer at the Max Planck Computing and Data Facility. Each node in this island consists of two 20-core Intel Xeon 6148 Gold sockets (Skylake architecture) with two Nvidia Volta 100 GPUs (PCIe); the nodes are connected with an Omnipath fabric (100 Gbit/s). When executed on a full node with GPUs, the average time step for the example is a factor of 4.8 faster than on one full CPU node, i.e., the time to solution is reduced by a factor of 4.8. This also means that the GPU version consumes less energy: although one GPU node draws more power than a CPU node (~950 W vs 450 W), the faster execution time reduces the overall energy to solution by a factor of 2.3 when running on a GPU node. When scaling to 2, 4, and 8 nodes, the speed-up is 1.26, 2.14, and 3.37, respectively, and the corresponding parallel efficiency is 63%, 54%, and 42%. Although the scaling is not yet perfect, it has improved considerably and also hints at the need for improving inter-node communication. For the second example (Fig. 16, right panel), a time-dependent simulation of a part of a chlorophyll complex was executed on a Supermicro server with two 8-core Intel Xeon 6134 Gold CPUs (Skylake architecture) and eight Nvidia Volta 100 GPUs interconnected with NVLink. Here, the speed-up of the average time step of the simulation when using 2, 4, and 8 GPUs is 1.95, 3.70, and 5.97, which corresponds to parallel efficiencies of 97%, 93%, and 75%, respectively. These efficiencies are much better than for the multi-node example, probably because data are only communicated within one node. For both examples, all states were stored in the GPU memory, and the parallelization was achieved by distributing the states only. Distributing the grid (i.e., parallelization over domain) is not as effective, because it leads to more frequent communication, and thus has a larger overhead.

<figure><img src="image_26.png" /><figcaption>FIG. 16. Scaling plots for the GPU implementation of Octopus. The left panel shows speed-up for a machine with two 20-core Intel Xeon 6148 Gold sockets (Skylake architecture) and two V100 GPUs (PCIe) per node. The right panel shows speed-up for a Supermicro server with 8 V100 GPUs (NVLink).</figcaption></figure>
This was typically a combined system of electrons and nuclei, where the former were described using DFT and the latter were treated classically. All possible sorts of interactions of this system with some external source (e.g., external electromagnetic fields, solvents described with the PCM, etc.) were then added in an *ad hoc* fashion. Although this approach worked well for many features and applications, its limitations became evident with several recent developments, in particular, the coupled Maxwell-Kohn-Sham equations described in Sec. II, where the Kohn-Sham orbitals need to be time-propagated alongside the Maxwell fields.

This has prompted us to start a major refactoring of the code, where the full physical system is treated as several subsystems interacting with each other. Such subsystems can be electrons, nuclei, and Maxwell fields. While this framework was mainly motivated by the Maxwell-TDDFT coupling, it has been implemented in a very general fashion so that all existing features can be converted to it and that, in the future, many other developments that require the coupling of several subsystems can be based on it.

## D. Memory layout

A new memory layout for storing the orbitals was introduced in Ref. 186, where all states are stored in a number of smaller batches with the innermost index being the state index instead of the real-space grid index. Thus, for all states in a batch, exactly the same operations can be executed when looping over the grid while accessing memory contiguously. This allows efficient parallelization on graphical processing units (GPUs), where all threads in a warp need to execute the same instructions, as well as vectorization on central processing units (CPUs), where one instruction can operate on several data points. To fully utilize these instructions, the kernel for computing finite differences (e.g., the Laplacian) has been specialized to explicitly use SSE, AVX, or AVX512 instructions. Now, more parts of the code have been ported to use this new layout to increase their performance, and whenever possible, the code will use this layout by default.

## E. GPUs

The first GPU implementation was described in Ref. 186. This was based on OpenCL and was limited to a small selection of numerical algorithms and calculation modes. These included some of the most commonly used features of the code or algorithms that were particularly suited for GPU porting, such as the residual minimization method-direct inversion in the iterative subspace (RMM-DIIS) eigensolver for ground-state calculations or the enforced time-reversal-symmetry propagator for time-dependent calculations. Since then, the GPU implementation has

been expanded in several different ways. Most notably, it now supports CUDA through an additional compatibility layer. The number of supported features and algorithms has also been increasing steadily such that most time-dependent calculations can now be run efficiently on GPUs. The implementation was also expanded to support multiple devices per host. Using the packed storage format described above, Octopus is able to store the states fully in the GPU memory, provided the memory is large enough, thus reducing memory transfers to a minimum. Recently, this was improved further by removing frequent allocations and deallocations of temporary variables on the GPU, now using a custom memory management for those. This proved very effective in improving the scaling to several GPUs and also to several nodes with GPUs.

We show two examples of time-dependent runs in Fig. 16. For the first example (left panel), β-cyclodextrine was used as an input system and the simulation was run on the GPU island of the COBRA supercomputer at the Max Planck Computing and Data Facility. Each node in this island consists of two 20-core Intel Xeon 6148 Gold sockets (Skylake architecture) with two Nvidia Volta 100 GPUs (PCIe); the nodes are connected with an Omnipath fabric (100 Gbit/s). When executed on a full node with GPUs, the average time step for the example is a factor of 4.8 faster than on one full CPU node, i.e., the time to solution is reduced by a factor of 4.8. This also means that the GPU version consumes less energy: although one GPU node draws more power than a CPU node (~950 W vs 450 W), the faster execution time reduces the overall energy to solution by a factor of 2.3 when running on a GPU node. When scaling to 2, 4, and 8 nodes, the speed-up is 1.26, 2.14, and 3.37, respectively, and the corresponding parallel efficiency is 63%, 54%, and 42%. Although the scaling is not yet perfect, it has improved considerably and also hints at the need for improving inter-node communication. For the second example (Fig. 16, right panel), a time-dependent simulation of a part of a chlorophyll complex was executed on a Supermicro server with two 8-core Intel Xeon 6134 Gold CPUs (Skylake architecture) and eight Nvidia Volta 100 GPUs interconnected with NVLink. Here, the speed-up of the average time step of the simulation when using 2, 4, and 8 GPUs is 1.95, 3.70, and 5.97, which corresponds to parallel efficiencies of 97%, 93%, and 75%, respectively. These efficiencies are much better than for the multi-node example, probably because data are only communicated within one node. For both examples, all states were stored in the GPU memory, and the parallelization was achieved by distributing the states only. Distributing the grid (i.e., parallelization over domain) is not as effective, because it leads to more frequent communication, and thus has a larger overhead.

<figure><img src="image_28.png" /><figcaption>FIG. 16. Scaling plots for the GPU implementation of Octopus. The left panel shows speed-up for a machine with two 20-core Intel Xeon 6148 Gold sockets (Skylake architecture) and two V100 GPUs (PCIe) per node. The right panel shows speed-up for a Supermicro server with 8 V100 GPUs (NVLink).</figcaption></figure>
Our current and future efforts regarding the GPU implementation are focused on two aspects. First, we are planning to enhance the existing implementation by improving the communication, possibly by overlapping computation and communication, and by optimizing the existing kernels. Second, we want to port more algorithms that are currently only implemented for the CPU version, especially for spin-dependent calculations and ground-state runs. This is a long-term effort, and the ultimate goal is to have as many features and algorithms as possible running efficiently on GPUs.

## XVII. CONCLUSIONS

It has been almost 20 years since the development of Octopus started. During this period, the code has matured and expanded to cover an ever-growing range of methods, theories, applications, and systems. It has also adapted continuously to the available computer architectures and computing paradigms, allowing researchers to tackle increasingly challenging problems in electronic structure theory.

Although many of the code capabilities have been routinely used by many groups around the world, mainly to study electronic excited state properties and dynamics, which still remain, to a large extent, the core of Octopus, we believe the main reason for the code's success lies elsewhere. From the beginning, the code has been designed to take full advantage of the flexibility and versatility offered by the use of real-space grids and provides developers with a framework to easily implement and test new ideas and methods that can be later on adopted by other codes (as it has already been the case with quite a few features previously developed within Octopus). This is demonstrated by the number of new theoretical methodologies and frameworks presented in this paper to deal with non-equilibrium phenomena of complex systems and to address the combined dynamics of electrons, phonons, and photons that go beyond what can be found in other electronic structure codes. This includes several novel approaches to treat the coupling of the electronic systems to the photons, such as the coupled Maxwell-Kohn-Sham equations, the OEP approach to the electron-photon coupling, and the dressed RDMFT. Other examples include the description of magnons in real-time and the orbital magneto-optical response in solids and molecules using the Sternheimer approach. In the case of magnons, supercells need to be employed, and the scalability of real-space grid methods make this approach very promising for future applications in more complex correlated magnetic systems. We expect that many of these methods and approaches will be integrated into other electronic-structure codes and will become standard tools in the near future.

The flexibility of the code is also demonstrated by the variety of systems that it can efficiently treat, as the applications described in this paper include molecules, nanoparticles, model systems, solvents, solids, and monolayers. Particularly noteworthy is the efficient treatment of periodic systems, which traditionally have been described using plane-wave basis-sets, and this has been achieved without loss of accuracy, as shown by our results for the delta-factor test. In the end, it is our hope that combining this flexibility with a growing range of methodologies will provide researchers with the necessary tools to study new challenging phenomena, such as novel correlated materials, out of equilibrium physics, or coupled electron-boson systems.

We also discussed recent improvements in performance and scalability, with particular emphasis on the GPU support. These developments are crucial in the view of the new challenges that electronic-structure applications are facing with the upcoming exaflop supercomputers.

Finally, several other implementations are in the pipeline, such as Floquet and cavity QED materials engineering, multitrajec-tory methods to deal with the nonadiabatic electron-ion dynamics, treatment of open quantum dissipative systems, and spectroscopies with entangled photons. All of these should be available to the users in the next years, so we invite you to stay tuned to the Octopus webpage at https://octopus-code.org/.

## SUPPLEMENTARY MATERIAL

The supplementary material contains the following: (i) Movie SM1—The nitrobenzene molecule is surrounded by the solvation cavity on which the PCM charges rest. PCM charges are modified in an iterative scheme coupled to the KS problem to get the molecular ground state in equilibrium with the solvent (here, water). The size and color of the spheres reflect the PCM charge value at each point of the surface. The animation is the 3D version of Fig. 7 in the main text. (ii) Movie SM2—Time dependent evolution of the molecular dipole of the nitrobenzene molecule in water subject to an electri-cal dipole perturbation in the x direction applied at the initial time. (Left) The PCM surface charges (represented as in movie SM1) and the molecular dipole (red arrow) are animated according to the equi-librium method (see main text for details). The white arrow is the molecular dipole *in vacuo*. (Right) The molecular dipole magnitude (top panel) and *xy* in-plane angle (bottom panel) calculated with the different PCM methods implemented in Octopus. Observe how the presence of the solvent enhances the intrinsic molecular dipole magnitude and produces a noticeable dephasing with respect to the *in vacuo* case already after few fs. This dephasing reflects into the solvatochromic shift reported in Fig. 8 of the main text.

## ACKNOWLEDGMENTS

The authors would like to thank all the people who have contributed to the development of Octopus over the last two decades. They would also like to thank Lin Lin for useful and interesting discussions and acknowledge the open discussions about real space methods with the group of Professor Chelikowsky. This work was supported by the European Research Council (Grant No. ERC-2015-AdG694097), the Cluster of Excellence “Advanced Imaging of Matter” (AIM), Grupos Consolidados (IT1249-19), and SFB925. The Flatiron Institute is a division of the Simons Foundation. X.A., A.W., and A.C. acknowledge that part of this work was performed under the auspices of the U.S. Department of Energy at Lawrence Livermore National Laboratory under Contract No. DE-AC52-07A27344. J.J.-S. gratefully acknowledges the funding from the European Union Horizon 2020 Research and Innovation Program under the Marie Sklodowska-Curie Grant Agreement No. 795246-StrongLights. J.F. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG Forschungsstipendium FL 997/1-1). D.A.S. acknowledges University of California, Merced start-up funding.
Our current and future efforts regarding the GPU implementation are focused on two aspects. First, we are planning to enhance the existing implementation by improving the communication, possibly by overlapping computation and communication, and by optimizing the existing kernels. Second, we want to port more algorithms that are currently only implemented for the CPU version, especially for spin-dependent calculations and ground-state runs. This is a long-term effort, and the ultimate goal is to have as many features and algorithms as possible running efficiently on GPUs.

## XVII. CONCLUSIONS

It has been almost 20 years since the development of Octopus started. During this period, the code has matured and expanded to cover an ever-growing range of methods, theories, applications, and systems. It has also adapted continuously to the available computer architectures and computing paradigms, allowing researchers to tackle increasingly challenging problems in electronic structure theory.

Although many of the code capabilities have been routinely used by many groups around the world, mainly to study electronic excited state properties and dynamics, which still remain, to a large extent, the core of Octopus, we believe the main reason for the code's success lies elsewhere. From the beginning, the code has been designed to take full advantage of the flexibility and versatility offered by the use of real-space grids and provides developers with a framework to easily implement and test new ideas and methods that can be later on adopted by other codes (as it has already been the case with quite a few features previously developed within Octopus). This is demonstrated by the number of new theoretical methodologies and frameworks presented in this paper to deal with non-equilibrium phenomena of complex systems and to address the combined dynamics of electrons, phonons, and photons that go beyond what can be found in other electronic structure codes. This includes several novel approaches to treat the coupling of the electronic systems to the photons, such as the coupled Maxwell-Kohn-Sham equations, the OEP approach to the electron-photon coupling, and the dressed RDMFT. Other examples include the description of magnons in real-time and the orbital magneto-optical response in solids and molecules using the Sternheimer approach. In the case of magnons, supercells need to be employed, and the scalability of real-space grid methods make this approach very promising for future applications in more complex correlated magnetic systems. We expect that many of these methods and approaches will be integrated into other electronic-structure codes and will become standard tools in the near future.

The flexibility of the code is also demonstrated by the variety of systems that it can efficiently treat, as the applications described in this paper include molecules, nanoparticles, model systems, solvents, solids, and monolayers. Particularly noteworthy is the efficient treatment of periodic systems, which traditionally have been described using plane-wave basis-sets, and this has been achieved without loss of accuracy, as shown by our results for the delta-factor test. In the end, it is our hope that combining this flexibility with a growing range of methodologies will provide researchers with the necessary tools to study new challenging phenomena, such as novel correlated materials, out of equilibrium physics, or coupled electron-boson systems.

We also discussed recent improvements in performance and scalability, with particular emphasis on the GPU support. These developments are crucial in the view of the new challenges that electronic-structure applications are facing with the upcoming exaflop supercomputers.

Finally, several other implementations are in the pipeline, such as Floquet and cavity QED materials engineering, multitrajec-tory methods to deal with the nonadiabatic electron-ion dynamics, treatment of open quantum dissipative systems, and spectroscopies with entangled photons. All of these should be available to the users in the next years, so we invite you to stay tuned to the Octopus webpage at https://octopus-code.org/.

## SUPPLEMENTARY MATERIAL

The supplementary material contains the following: (i) Movie SM1—The nitrobenzene molecule is surrounded by the solvation cavity on which the PCM charges rest. PCM charges are modified in an iterative scheme coupled to the KS problem to get the molecular ground state in equilibrium with the solvent (here, water). The size and color of the spheres reflect the PCM charge value at each point of the surface. The animation is the 3D version of Fig. 7 in the main text. (ii) Movie SM2—Time dependent evolution of the molecular dipole of the nitrobenzene molecule in water subject to an electric-cal dipole perturbation in the x direction applied at the initial time. (Left) The PCM surface charges (represented as in movie SM1) and the molecular dipole (red arrow) are animated according to the equi-librium method (see main text for details). The white arrow is the molecular dipole *in vacuo*. (Right) The molecular dipole magnitude (top panel) and *xy* in-plane angle (bottom panel) calculated with the different PCM methods implemented in Octopus. Observe how the presence of the solvent enhances the intrinsic molecular dipole magnitude and produces a noticeable dephasing with respect to the *in vacuo* case already after few fs. This dephasing reflects into the solvatochromic shift reported in Fig. 8 of the main text.

## ACKNOWLEDGMENTS

The authors would like to thank all the people who have contributed to the development of Octopus over the last two decades. They would also like to thank Lin Lin for useful and interesting discussions and acknowledge the open discussions about real space methods with the group of Professor Chelikowsky. This work was supported by the European Research Council (Grant No. ERC-2015-AdG694097), the Cluster of Excellence “Advanced Imaging of Matter” (AIM), Grupos Consolidados (IT1249-19), and SFB925. The Flatiron Institute is a division of the Simons Foundation. X.A., A.W., and A.C. acknowledge that part of this work was performed under the auspices of the U.S. Department of Energy at Lawrence Livermore National Laboratory under Contract No. DE-AC52-07A27344. J.J.-S. gratefully acknowledges the funding from the European Union Horizon 2020 Research and Innovation Program under the Marie Sklodowska-Curie Grant Agreement No. 795246-StrongLights. J.F. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG Forschungsstipendium FL 997/1-1). D.A.S. acknowledges University of California, Merced start-up funding.
REFERENCES

¹M. A. L. Marques, A. Rubio, E. K. Gross, K. Burke, F. Nogueira, and C. A. Ullrich, *Time-Dependent Density Functional Theory* (Springer Science & Business Media, 2006), Vol. 706.

²M. A. L. Marques, N. T. Maitra, F. M. Nogueira, E. K. Gross, and A. Rubio, *Fundamentals of Time-Dependent Density Functional Theory* (Springer Science & Business Media, 2012), Vol. 837.

³T. W. Ebbesen, “Hybrid light–matter states in a molecular and material science perspective,” *Acc. Chem. Res.* **49**, 2403–2412 (2016).

⁴M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly, and A. Rubio, “Quantum-electrodynamical density-functional theory: Bridging quantum optics and electronic-structure theory,” *Phys. Rev. A* **90**, 012508 (2014).

⁵J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, “Kohn-Sham approach to quantum electrodynamical density-functional theory: Exact time-dependent effective potentials in real space,” *Proc. Natl. Acad. Sci. U. S. A.* **112**, 15285–15290 (2015).

⁶J. Flick, N. Rivera, and P. Narang, “Strong light-matter coupling in quantum chemistry and quantum photonics,” *Nanophotonics* **7**, 1479–1501 (2018).

⁷J. Feist, J. Galego, and F. J. Garcia-Vidal, “Polaritonic chemistry with organic molecules,” *ACS Photonics* **5**, 205–216 (2017).

⁸R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, “Polariton chemistry: Controlling molecular dynamics with optical cavities,” *Chem. Sci.* **9**, 6325–6339 (2018).

⁹A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, “Ultrastrong coupling between light and matter,” *Nat. Rev. Phys.* **1**, 19–40 (2019).

¹⁰I. V. Tokatly, “Time-dependent density functional theory for many-electron systems interacting with cavity photons,” *Phys. Rev. Lett.* **110**, 233001 (2013).

¹¹J. Flick, “Exact nonadiabatic many-body dynamics: Electron-phonon coupling in photoelectron spectroscopy and light-matter interactions in quantum electrodynamical density-functional theory,” Ph.D. thesis, Humboldt-Universität zu Berlin, Berlin, 2016.

¹²M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, “From a quantum-electrodynamical light-matter description to novel spectroscopies,” *Nat. Rev. Chem.* **2**, 0118 (2018).

¹³J. Flick, C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, “Ab initio optimized effective potentials for real molecules in optical cavities: Photon contributions to the molecular ground state,” *ACS Photonics* **5**, 992–1005 (2018).

¹⁴C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, “Modification of excitation and charge transfer in cavity quantum-electrodynamical chemistry,” *Proc. Natl. Acad. Sci. U. S. A.* **116**, 4883–4892 (2019).

¹⁵A. Castro, “A first-principles time-dependent density functional theory scheme for the computation of the electromagnetic response of nanostructures,” Ph.D. thesis, University of Valladolid, 2004.

¹⁶X. Andrade, “Linear and non-linear response phenomena of molecular systems within time-dependent density functional theory,” Ph.D. thesis, Universidad del País Vasco, 2010.

¹⁷M. A. Marques, A. Castro, G. F. Bertsch, and A. Rubio, “Octopus: A first-principles tool for excited electron–ion dynamics,” *Comput. Phys. Commun.* **151**, 60–78 (2003).

¹⁸A. Castro, H. Appel, M. Oliveira, C. A. Rozzi, X. Andrade, F. Lorenzen, M. A. L. Marques, E. K. U. Gross, and A. Rubio, “Octopus: A tool for the application of time-dependent density functional theory,” *Phys. Status Solidi B* **243**, 2465–2488 (2006).

¹⁹X. Andrade, D. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. J. Verstraete, L. Stella, F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and A. Rubio, “Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems,” *Phys. Chem. Chem. Phys.* **17**, 31371–31396 (2015).

²⁰See https://octopus-code.org/ for the Octopus code, tutorials, and examples can be found on its website.

²¹R. Jestädt, M. Ruggenthaler, M. J. T. Oliveira, A. Rubio, and H. Appel, “Light-matter interactions within the Ehrenfest-Maxwell-Pauli-Kohn-Sham framework:

Fundamentals, implementation, and nano-optical applications,” *Adv. Phys.* **68**, 225–333 (2019); arXiv:1812.05049.

²²G. F. Bertsch, J.-I. Iwata, A. Rubio, and K. Yabana, “Real-space, real-time method for the dielectric function,” *Phys. Rev. B* **62**, 7998–8002 (2000).

²³A. Taflove and S. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, Artech House Antennas and Propagation Library (Artech House, 2005).

²⁴I. Bialynicki-Birula, “On the wave function of the photon,” *Acta Phys. Pol., Ser. A* **86**, 97–116 (1994).

²⁵S. Ludwig, “Elektromagnetische grundgleichungen in bivektorieller behandlung,” *Ann. Phys.* **327**, 579–586 (1907).

²⁶R. Loudon, *The Quantum Theory of Light* (Oxford Science Publications, 1988).

²⁷D. P. Craig and T. Thirunamachandran, *Molecular Quantum Electrodynamics: An Introduction to Radiation-Molecule Interactions* (Courier Corporation, 1998).

²⁸J. J. Baumberg, J. Aizpurua, M. H. Mikkelsen, and D. R. Smith, “Extreme nanophotonics from ultrathin metallic gaps,” *Nat. Mater.* **18**, 668–678 (2019).

²⁹R. van Leeuwen, “The Sham-Schlüter equation in time-dependent density-functional theory,” *Phys. Rev. Lett.* **76**, 3610–3613 (1996).

³⁰C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel, and A. Rubio, “Optimized effective potential for quantum electrodynamical time-dependent density functional theory,” *Phys. Rev. Lett.* **115**, 093001 (2015).

³¹S. Kümmel and L. Kronik, “Orbital-dependent density functionals: Theory and applications,” *Rev. Mod. Phys.* **80**, 3–60 (2008).

³²M. Mundt and S. Kümmel, “Optimized effective potential in real time: Problems and prospects in time-dependent density-functional theory,” *Phys. Rev. A* **74**, 022511 (2006).

³³H. O. Wijewardane and C. A. Ullrich, “Real-time electron dynamics with exact-exchange time-dependent density-functional theory,” *Phys. Rev. Lett.* **100**, 056404 (2008).

³⁴S.-L. Liao, T.-S. Ho, H. Rabitz, and S.-I. Chu, “Time-local equation for the exact optimized effective potential in time-dependent density functional theory,” *Phys. Rev. Lett.* **118**, 243001 (2017).

³⁵J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, “Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry,” *Proc. Natl. Acad. Sci. U. S. A.* **114**, 3026–3034 (2017).

³⁶S. E. B. Nielsen, C. Schäfer, M. Ruggenthaler, and A. Rubio, “Dressed-orbital approach to cavity quantum electrodynamics and beyond,” arXiv:1812.00388 (2018).

³⁷S. Kümmel and J. P. Perdew, “Simple iterative construction of the optimized effective potential for orbital functionals, including exact exchange,” *Phys. Rev. Lett.* **90**, 043004 (2003).

³⁸T. W. Hollins, S. J. Clark, K. Refson, and N. I. Gidopoulos, “Optimized effective potential using the hylleraas variational method,” *Phys. Rev. B* **85**, 235126 (2012).

³⁹J. Krieger, Y. Li, and G. Iafrate, “Derivation and application of an accurate Kohn-Sham potential with integer discontinuity,” *Phys. Lett. A* **146**, 256–260 (1990).

⁴⁰J. B. Krieger, Y. Li, and G. J. Iafrate, “Construction and application of an accurate local spin-polarized Kohn-Sham potential with integer discontinuity: Exchange-only theory,” *Phys. Rev. A* **45**, 101–126 (1992).

⁴¹J. B. Krieger, Y. Li, and G. J. Iafrate, “Systematic approximations to the optimized effective potential: Application to orbital-density-functional theory,” *Phys. Rev. A* **46**, 5453–5458 (1992).

⁴²C. Schäfer, M. Ruggenthaler, and A. Rubio, “Ab initio nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling,” *Phys. Rev. A* **98**, 043801 (2018).

⁴³V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, “Light-matter interaction in the long-wavelength limit: No ground-state without dipole self-energy,” *J. Phys. B: At., Mol. Opt. Phys.* **51**, 034005 (2018).

⁴⁴J. Flick and P. Narang, “Cavity-correlated electron-nuclear dynamics from first principles,” *Phys. Rev. Lett.* **121**, 113002 (2018).

⁴⁵J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel, and A. Rubio, “Light-matter response in nonrelativistic quantum electrodynamics,” *ACS Photonics* **6**, 2757–2778 (2019).
REFERENCES

¹M. A. L. Marques, A. Rubio, E. K. Gross, K. Burke, F. Nogueira, and C. A. Ullrich, *Time-Dependent Density Functional Theory* (Springer Science & Business Media, 2006), Vol. 706.

²M. A. L. Marques, N. T. Maitra, F. M. Nogueira, E. K. Gross, and A. Rubio, *Fundamentals of Time-Dependent Density Functional Theory* (Springer Science & Business Media, 2012), Vol. 837.

³T. W. Ebbesen, “Hybrid light–matter states in a molecular and material science perspective,” *Acc. Chem. Res.* **49**, 2403–2412 (2016).

⁴M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly, and A. Rubio, “Quantum-electrodynamical density-functional theory: Bridging quantum optics and electronic-structure theory,” *Phys. Rev. A* **90**, 012508 (2014).

⁵J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, “Kohn-Sham approach to quantum electrodynamical density-functional theory: Exact time-dependent effective potentials in real space,” *Proc. Natl. Acad. Sci. U. S. A.* **112**, 15285–15290 (2015).

⁶J. Flick, N. Rivera, and P. Narang, “Strong light-matter coupling in quantum chemistry and quantum photonics,” *Nanophotonics* **7**, 1479–1501 (2018).

⁷J. Feist, J. Galego, and F. J. Garcia-Vidal, “Polaritonic chemistry with organic molecules,” *ACS Photonics* **5**, 205–216 (2017).

⁸R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo, and J. Yuen-Zhou, “Polariton chemistry: Controlling molecular dynamics with optical cavities,” *Chem. Sci.* **9**, 6325–6339 (2018).

⁹A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, “Ultrastrong coupling between light and matter,” *Nat. Rev. Phys.* **1**, 19–40 (2019).

¹⁰I. V. Tokatly, “Time-dependent density functional theory for many-electron systems interacting with cavity photons,” *Phys. Rev. Lett.* **110**, 233001 (2013).

¹¹J. Flick, “Exact nonadiabatic many-body dynamics: Electron-phonon coupling in photoelectron spectroscopy and light-matter interactions in quantum electrodynamical density-functional theory,” Ph.D. thesis, Humboldt-Universität zu Berlin, Berlin, 2016.

¹²M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, “From a quantum-electrodynamical light-matter description to novel spectroscopies,” *Nat. Rev. Chem.* **2**, 0118 (2018).

¹³J. Flick, C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, “Ab initio optimized effective potentials for real molecules in optical cavities: Photon contributions to the molecular ground state,” *ACS Photonics* **5**, 992–1005 (2018).

¹⁴C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, “Modification of excitation and charge transfer in cavity quantum-electrodynamical chemistry,” *Proc. Natl. Acad. Sci. U. S. A.* **116**, 4883–4892 (2019).

¹⁵A. Castro, “A first-principles time-dependent density functional theory scheme for the computation of the electromagnetic response of nanostructures,” Ph.D. thesis, University of Valladolid, 2004.

¹⁶X. Andrade, “Linear and non-linear response phenomena of molecular systems within time-dependent density functional theory,” Ph.D. thesis, Universidad del País Vasco, 2010.

¹⁷M. A. Marques, A. Castro, G. F. Bertsch, and A. Rubio, “Octopus: A first-principles tool for excited electron–ion dynamics,” *Comput. Phys. Commun.* **151**, 60–78 (2003).

¹⁸A. Castro, H. Appel, M. Oliveira, C. A. Rozzi, X. Andrade, F. Lorenzen, M. A. L. Marques, E. K. U. Gross, and A. Rubio, “Octopus: A tool for the application of time-dependent density functional theory,” *Phys. Status Solidi B* **243**, 2465–2488 (2006).

¹⁹X. Andrade, D. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. J. Verstraete, L. Stella, F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and A. Rubio, “Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems,” *Phys. Chem. Chem. Phys.* **17**, 31371–31396 (2015).

²⁰See https://octopus-code.org/ for the Octopus code, tutorials, and examples can be found on its website.

²¹R. Jestädt, M. Ruggenthaler, M. J. T. Oliveira, A. Rubio, and H. Appel, “Light-matter interactions within the Ehrenfest-Maxwell-Pauli-Kohn-Sham framework:

Fundamentals, implementation, and nano-optical applications,” *Adv. Phys.* **68**, 225–333 (2019); arXiv:1812.05049.

²²G. F. Bertsch, J.-I. Iwata, A. Rubio, and K. Yabana, “Real-space, real-time method for the dielectric function,” *Phys. Rev. B* **62**, 7998–8002 (2000).

²³A. Taflove and S. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, Artech House Antennas and Propagation Library (Artech House, 2005).

²⁴I. Bialynicki-Birula, “On the wave function of the photon,” *Acta Phys. Pol., Ser. A* **86**, 97–116 (1994).

²⁵S. Ludwig, “Elektromagnetische grundgleichungen in bivektorieller behandlung,” *Ann. Phys.* **327**, 579–586 (1907).

²⁶R. Loudon, *The Quantum Theory of Light* (Oxford Science Publications, 1988).

²⁷D. P. Craig and T. Thirunamachandran, *Molecular Quantum Electrodynamics: An Introduction to Radiation-Molecule Interactions* (Courier Corporation, 1998).

²⁸J. J. Baumberg, J. Aizpurua, M. H. Mikkelsen, and D. R. Smith, “Extreme nanophotonics from ultrathin metallic gaps,” *Nat. Mater.* **18**, 668–678 (2019).

²⁹R. van Leeuwen, “The Sham-Schlüter equation in time-dependent density-functional theory,” *Phys. Rev. Lett.* **76**, 3610–3613 (1996).

³⁰C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel, and A. Rubio, “Optimized effective potential for quantum electrodynamical time-dependent density functional theory,” *Phys. Rev. Lett.* **115**, 093001 (2015).

³¹S. Kümmel and L. Kronik, “Orbital-dependent density functionals: Theory and applications,” *Rev. Mod. Phys.* **80**, 3–60 (2008).

³²M. Mundt and S. Kümmel, “Optimized effective potential in real time: Problems and prospects in time-dependent density-functional theory,” *Phys. Rev. A* **74**, 022511 (2006).

³³H. O. Wijewardane and C. A. Ullrich, “Real-time electron dynamics with exact-exchange time-dependent density-functional theory,” *Phys. Rev. Lett.* **100**, 056404 (2008).

³⁴S.-L. Liao, T.-S. Ho, H. Rabitz, and S.-I. Chu, “Time-local equation for the exact optimized effective potential in time-dependent density functional theory,” *Phys. Rev. Lett.* **118**, 243001 (2017).

³⁵J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, “Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry,” *Proc. Natl. Acad. Sci. U. S. A.* **114**, 3026–3034 (2017).

³⁶S. E. B. Nielsen, C. Schäfer, M. Ruggenthaler, and A. Rubio, “Dressed-orbital approach to cavity quantum electrodynamics and beyond,” arXiv:1812.00388 (2018).

³⁷S. Kümmel and J. P. Perdew, “Simple iterative construction of the optimized effective potential for orbital functionals, including exact exchange,” *Phys. Rev. Lett.* **90**, 043004 (2003).

³⁸T. W. Hollins, S. J. Clark, K. Refson, and N. I. Gidopoulos, “Optimized effective potential using the hylleraas variational method,” *Phys. Rev. B* **85**, 235126 (2012).

³⁹J. Krieger, Y. Li, and G. Iafrate, “Derivation and application of an accurate Kohn-Sham potential with integer discontinuity,” *Phys. Lett. A* **146**, 256–260 (1990).

⁴⁰J. B. Krieger, Y. Li, and G. J. Iafrate, “Construction and application of an accurate local spin-polarized Kohn-Sham potential with integer discontinuity: Exchange-only theory,” *Phys. Rev. A* **45**, 101–126 (1992).

⁴¹J. B. Krieger, Y. Li, and G. J. Iafrate, “Systematic approximations to the optimized effective potential: Application to orbital-density-functional theory,” *Phys. Rev. A* **46**, 5453–5458 (1992).

⁴²C. Schäfer, M. Ruggenthaler, and A. Rubio, “Ab initio nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling,” *Phys. Rev. A* **98**, 043801 (2018).

⁴³V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, “Light-matter interaction in the long-wavelength limit: No ground-state without dipole self-energy,” *J. Phys. B: At., Mol. Opt. Phys.* **51**, 034005 (2018).

⁴⁴J. Flick and P. Narang, “Cavity-correlated electron-nuclear dynamics from first principles,” *Phys. Rev. Lett.* **121**, 113002 (2018).

⁴⁵J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel, and A. Rubio, “Light-matter response in nonrelativistic quantum electrodynamics,” *ACS Photonics* **6**, 2757–2778 (2019).
<footnote n="46">⁴⁶N. M. Hoffmann, C. Schäfer, A. Rubio, A. Kelly, and H. Appel, "Capturing vacuum fluctuations and photon correlations in cavity quantum electrodynamics with multitrajectory ehrenfest dynamics," Phys. Rev. A 99, 063819 (2019).</footnote>

<footnote n="47">⁴⁷J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, "Cavity Casimir-Polder forces and their effects in ground state chemical reactivity," Phys. Rev. X 9, 021057 (2019).</footnote>

<footnote n="48">⁴⁸Note that Hartree-Fock theory is also included within RDMFT by fixing the orbital occupations to 0 and 1.</footnote>

<footnote n="49">⁴⁹F. Buchholz, I. Theophilou, S. E. B. Nielsen, M. Ruggenthaler, and A. Rubio, "Reduced density-matrix approach to strong matter-photon interaction," ACS Photonics 6, 2694-2711 (2019).</footnote>

<footnote n="50">⁵⁰At the current state of the code, we have only implemented the so-called Müller functional,¹⁷¹ but one could potentially use any RDMFT functional to approximate two-body expression.</footnote>

<footnote n="51">⁵¹A. J. Coleman, "Structure of fermion density matrices," Rev. Mod. Phys. 35, 668-686 (1963).</footnote>

<footnote n="52">⁵²The respective results for the equilibrium distance *b* = 1.61 bohr are shown in Ref. 50, Sec. 7.</footnote>

<footnote n="53">⁵³Performed with the many-body routine of Octopus, see Ref. 19, Sec. 13 for details.</footnote>

<footnote n="54">⁵⁴V. I. Anisimov, J. Zaanen, and O. K. Andersen, "Band theory and Mott insulators: Hubbard U instead of Stoner I," Phys. Rev. B 44, 943 (1991).</footnote>

<footnote n="55">⁵⁵V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, "Density-functional theory and NiO photoemission spectra," Phys. Rev. B 48, 16929 (1993).</footnote>

<footnote n="56">⁵⁶A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, "Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators," Phys. Rev. B 52, R5467 (1995).</footnote>

<footnote n="57">⁵⁷V. I. Anisimov, F. Aryasetiawan, and A. Lichtenstein, "First-principles calculations of the electronic structure and spectra of strongly correlated systems: The LDA+U method," J. Phys.: Condens. Matter 9, 767 (1997).</footnote>

<footnote n="58">⁵⁸K. Haule, "Exact double counting in combining the dynamical mean field theory and the density functional theory," Phys. Rev. Lett. 115, 196403 (2015).</footnote>

<footnote n="59">⁵⁹A. G. Petukhov, I. I. Mazin, L. Chioncel, and A. I. Lichtenstein, "Correlated metals and the LDA+U method," Phys. Rev. B 67, 153106 (2003).</footnote>

<footnote n="60">⁶⁰S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, "Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study," Phys. Rev. B 57, 1505-1509 (1998).</footnote>

<footnote n="61">⁶¹M. T. Czyżyk and G. A. Sawatzky, "Local-density functional and on-site correlations: The electronic structure of La₂CuO₄ and LaCuO₃," Phys. Rev. B 49, 14211-14228 (1994).</footnote>

<footnote n="62">⁶²N. Tancogne-Dejean, M. J. T. Oliveira, and A. Rubio, "Self-consistent DFT+U method for real-space time-dependent density functional theory calculations," Phys. Rev. B 96, 245133 (2017).</footnote>

<footnote n="63">⁶³L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, "Multiflat bands and strong correlations in twisted bilayer boron nitride: Doping-induced correlated insulator and superconductor," Nano Lett. 19, 4934-4940 (2019).</footnote>

<footnote n="64">⁶⁴N. Tancogne-Dejean and A. Rubio, "Parameter-free hybrid functional based on an extended Hubbard model: DFT+U+V," arXiv:1911.10813 (2019).</footnote>

<footnote n="65">⁶⁵L. A. Agapito, S. Curtarolo, and M. Buongiorno Nardelli, "Reformulation of DFT+U as a pseudohybrid hubbard density functional for accelerated materials discovery," Phys. Rev. X 5, 011006 (2015).</footnote>

<footnote n="66">⁶⁶G. E. Topp, N. Tancogne-Dejean, A. F. Kemper, A. Rubio, and M. A. Sentef, "All-optical nonequilibrium pathway to stabilising magnetic Weyl semimetals in pyrochlore iridates," Nature Commun. 9, 4452 (2018).</footnote>

<footnote n="67">⁶⁷N. Tancogne-Dejean, M. A. Sentef, and A. Rubio, "Ultrafast modification of hubbard *U* in a strongly correlated material: Ab initio high-harmonic generation in NiO," Phys. Rev. Lett. 121, 097402 (2018).</footnote>

<footnote n="68">⁶⁸K. Berland, V. R. Cooper, K. Lee, E. Schröder, T. Thonhauser, P. Hyldgaard, and B. I. Lundqvist, "van der Waals forces in density functional theory: A review of the vdW-DF method," Rep. Prog. Phys. 78, 066501 (2015).</footnote>

<footnote n="69">⁶⁹S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, "A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu," J. Chem. Phys. 132, 154104 (2010).</footnote>

<footnote n="70">⁷⁰A. Tkatchenko and M. Scheffler, "Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data," Phys. Rev. Lett. 102, 073005 (2009).</footnote>

<footnote n="71">⁷¹A. H. Larsen, M. Kuisma, J. Löfgren, Y. Pouillon, P. Erhart, and P. Hyldgaard, "libvdwxc: A library for exchange-correlation functionals in the vdW-DF family," Modell. Simul. Mater. Sci. Eng. 25, 065004 (2017).</footnote>

<footnote n="72">⁷²J. P. Perdew and Y. Wang, "Accurate and simple analytic representation of the electron-gas correlation energy," Phys. Rev. B 45, 13244-13249 (1992).</footnote>

<footnote n="73">⁷³G. Román-Pérez and J. M. Soler, "Efficient implementation of a van der Waals density functional: Application to double-wall carbon nanotubes," Phys. Rev. Lett. 103, 096102 (2009).</footnote>

<footnote n="74">⁷⁴M. Frigo and S. G. Johnson, "The design and implementation of FFTW3," Proc. IEEE 93, 216-231 (2005), part of special issue: Program generation, optimization, and platform adaptation.</footnote>

<footnote n="75">⁷⁵M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, "van der Waals density functional for general geometries," Phys. Rev. Lett. 92, 246401 (2004).</footnote>

<footnote n="76">⁷⁶M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, "Erratum: "van der Waals density functional for general geometries" [Phys. Rev. Lett. 92, 246401 (2004)]," Phys. Rev. Lett. 95, 109902 (2005).</footnote>

<footnote n="77">⁷⁷K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, "Higher-accuracy van der Waals density functional," Phys. Rev. B 82, 081101 (2010).</footnote>

<footnote n="78">⁷⁸K. Berland and P. Hyldgaard, "Exchange functional that tests the robustness of the plasmon description of the van der Waals density functional," Phys. Rev. B 89, 035412 (2014).</footnote>

<footnote n="79">⁷⁹J. Klimeš, D. R. Bowler, and A. Michaelides, "Chemical accuracy for the van der Waals density functional," J. Phys.: Condens. Matter 22, 022201 (2010).</footnote>

<footnote n="80">⁸⁰V. R. Cooper, "van der Waals density functional: An appropriate exchange functional," Phys. Rev. B 81, 161104 (2010).</footnote>

<footnote n="81">⁸¹J. Tomasi, B. Mennucci, and R. Cammi, "Quantum mechanical continuum solvation models," Chem. Rev. 105, 2999-3094 (2005).</footnote>

<footnote n="82">⁸²J. Tomasi and M. Persico, "Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent," Chem. Rev. 94, 2027-2094 (1994).</footnote>

<footnote n="83">⁸³J. Tomasi, "The physical model," in Continuum Solvation Models in Chemical Physics: From Theory to Applications, edited by B. Mennucci and R. Cammi (Wiley & Sons, Chichester, UK, 2007), Chap. 1.1, p. 16.</footnote>

<footnote n="84">⁸⁴A. Delgado, S. Corni, S. Pittalis, and C. A. Rozzi, "Modeling solvation effects in real-space and real-time within density functional approaches," J. Chem. Phys. 143, 144111 (2015).</footnote>

<footnote n="85">⁸⁵E. Cancés, B. Mennucci, and J. Tomasi, "A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics," J. Chem. Phys. 107, 3032 (1997).</footnote>

<footnote n="86">⁸⁶J.-L. Pascual-ahuir, E. Silla, and I. Tunon, "GEPOL: An improved description of molecular surfaces. III. A new algorithm for the computation of a solvent-excluding surface," J. Comput. Chem. 15, 1127-1138 (1994).</footnote>

<footnote n="87">⁸⁷M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su et al., "General atomic and molecular electronic structure system," J. Comput. Chem. 14, 1347-1363 (1993).</footnote>

<footnote n="88">⁸⁸J. Perdew, K. Burke, and M. Ernzerhof, "Perdew, Burke, and Ernzerhof reply," Phys. Rev. Lett. 80, 891 (1998).</footnote>

<footnote n="89">⁸⁹M. Caricato, F. Ingrosso, B. Mennucci, and J. Tomasi, "A time-dependent polarizable continuum model: Theory and application," J. Chem. Phys. 122, 154501 (2005).</footnote>

<footnote n="90">⁹⁰R. Cammi and J. Tomasi, "Nonequilibrium solvation theory for the polarizable continuum model: A new formulation at the SCF level with application to the case of the frequency-dependent linear electric response function," Int. J. Quantum Chem. 56, 465-474 (1995).</footnote>

<footnote n="91">⁹¹S. Corni, S. Pipolo, and R. Cammi, "Equation of motion for the solvent polarization apparent charges in the polarizable continuum model: Application to real-time TDDFT," J. Phys. Chem. A 119, 5405-5416 (2014).</footnote>
<footnote n="46">⁴⁶N. M. Hoffmann, C. Schäfer, A. Rubio, A. Kelly, and H. Appel, "Capturing vacuum fluctuations and photon correlations in cavity quantum electrodynamics with multitrajectory ehrenfest dynamics," Phys. Rev. A 99, 063819 (2019).</footnote>

<footnote n="47">⁴⁷J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, "Cavity Casimir-Polder forces and their effects in ground state chemical reactivity," Phys. Rev. X 9, 021057 (2019).</footnote>

<footnote n="48">⁴⁸Note that Hartree-Fock theory is also included within RDMFT by fixing the orbital occupations to 0 and 1.</footnote>

<footnote n="49">⁴⁹F. Buchholz, I. Theophilou, S. E. B. Nielsen, M. Ruggenthaler, and A. Rubio, "Reduced density-matrix approach to strong matter-photon interaction," ACS Photonics 6, 2694-2711 (2019).</footnote>

<footnote n="50">⁵⁰At the current state of the code, we have only implemented the so-called Müller functional,¹⁷¹ but one could potentially use any RDMFT functional to approximate two-body expression.</footnote>

<footnote n="51">⁵¹A. J. Coleman, "Structure of fermion density matrices," Rev. Mod. Phys. 35, 668-686 (1963).</footnote>

<footnote n="52">⁵²The respective results for the equilibrium distance *b* = 1.61 bohr are shown in Ref. 50, Sec. 7.</footnote>

<footnote n="53">⁵³Performed with the many-body routine of Octopus, see Ref. 19, Sec. 13 for details.</footnote>

<footnote n="54">⁵⁴V. I. Anisimov, J. Zaanen, and O. K. Andersen, "Band theory and Mott insulators: Hubbard U instead of Stoner I," Phys. Rev. B 44, 943 (1991).</footnote>

<footnote n="55">⁵⁵V. I. Anisimov, I. V. Solovyev, M. A. Korotin, M. T. Czyżyk, and G. A. Sawatzky, "Density-functional theory and NiO photoemission spectra," Phys. Rev. B 48, 16929 (1993).</footnote>

<footnote n="56">⁵⁶A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, "Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators," Phys. Rev. B 52, R5467 (1995).</footnote>

<footnote n="57">⁵⁷V. I. Anisimov, F. Aryasetiawan, and A. Lichtenstein, "First-principles calculations of the electronic structure and spectra of strongly correlated systems: The LDA+U method," J. Phys.: Condens. Matter 9, 767 (1997).</footnote>

<footnote n="58">⁵⁸K. Haule, "Exact double counting in combining the dynamical mean field theory and the density functional theory," Phys. Rev. Lett. 115, 196403 (2015).</footnote>

<footnote n="59">⁵⁹A. G. Petukhov, I. I. Mazin, L. Chioncel, and A. I. Lichtenstein, "Correlated metals and the LDA+U method," Phys. Rev. B 67, 153106 (2003).</footnote>

<footnote n="60">⁶⁰S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, "Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study," Phys. Rev. B 57, 1505-1509 (1998).</footnote>

<footnote n="61">⁶¹M. T. Czyżyk and G. A. Sawatzky, "Local-density functional and on-site correlations: The electronic structure of La₂CuO₄ and LaCuO₃," Phys. Rev. B 49, 14211-14228 (1994).</footnote>

<footnote n="62">⁶²N. Tancogne-Dejean, M. J. T. Oliveira, and A. Rubio, "Self-consistent DFT+U method for real-space time-dependent density functional theory calculations," Phys. Rev. B 96, 245133 (2017).</footnote>

<footnote n="63">⁶³L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, "Multiflat bands and strong correlations in twisted bilayer boron nitride: Doping-induced correlated insulator and superconductor," Nano Lett. 19, 4934-4940 (2019).</footnote>

<footnote n="64">⁶⁴N. Tancogne-Dejean and A. Rubio, "Parameter-free hybrid functional based on an extended Hubbard model: DFT+U+V," arXiv:1911.10813 (2019).</footnote>

<footnote n="65">⁶⁵L. A. Agapito, S. Curtarolo, and M. Buongiorno Nardelli, "Reformulation of DFT+U as a pseudohybrid hubbard density functional for accelerated materials discovery," Phys. Rev. X 5, 011006 (2015).</footnote>

<footnote n="66">⁶⁶G. E. Topp, N. Tancogne-Dejean, A. F. Kemper, A. Rubio, and M. A. Sentef, "All-optical nonequilibrium pathway to stabilising magnetic Weyl semimetals in pyrochlore iridates," Nature Commun. 9, 4452 (2018).</footnote>

<footnote n="67">⁶⁷N. Tancogne-Dejean, M. A. Sentef, and A. Rubio, "Ultrafast modification of hubbard *U* in a strongly correlated material: Ab initio high-harmonic generation in NiO," Phys. Rev. Lett. 121, 097402 (2018).</footnote>

<footnote n="68">⁶⁸K. Berland, V. R. Cooper, K. Lee, E. Schröder, T. Thonhauser, P. Hyldgaard, and B. I. Lundqvist, "van der Waals forces in density functional theory: A review of the vdW-DF method," Rep. Prog. Phys. 78, 066501 (2015).</footnote>

<footnote n="69">⁶⁹S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, "A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu," J. Chem. Phys. 132, 154104 (2010).</footnote>

<footnote n="70">⁷⁰A. Tkatchenko and M. Scheffler, "Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data," Phys. Rev. Lett. 102, 073005 (2009).</footnote>

<footnote n="71">⁷¹A. H. Larsen, M. Kuisma, J. Löfgren, Y. Pouillon, P. Erhart, and P. Hyldgaard, "libvdwxc: A library for exchange-correlation functionals in the vdW-DF family," Modell. Simul. Mater. Sci. Eng. 25, 065004 (2017).</footnote>

<footnote n="72">⁷²J. P. Perdew and Y. Wang, "Accurate and simple analytic representation of the electron-gas correlation energy," Phys. Rev. B 45, 13244-13249 (1992).</footnote>

<footnote n="73">⁷³G. Román-Pérez and J. M. Soler, "Efficient implementation of a van der Waals density functional: Application to double-wall carbon nanotubes," Phys. Rev. Lett. 103, 096102 (2009).</footnote>

<footnote n="74">⁷⁴M. Frigo and S. G. Johnson, "The design and implementation of FFTW3," Proc. IEEE 93, 216-231 (2005), part of special issue: Program generation, optimization, and platform adaptation.</footnote>

<footnote n="75">⁷⁵M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, "van der Waals density functional for general geometries," Phys. Rev. Lett. 92, 246401 (2004).</footnote>

<footnote n="76">⁷⁶M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, and B. I. Lundqvist, "Erratum: "van der Waals density functional for general geometries" [Phys. Rev. Lett. 92, 246401 (2004)]," Phys. Rev. Lett. 95, 109902 (2005).</footnote>

<footnote n="77">⁷⁷K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, "Higher-accuracy van der Waals density functional," Phys. Rev. B 82, 081101 (2010).</footnote>

<footnote n="78">⁷⁸K. Berland and P. Hyldgaard, "Exchange functional that tests the robustness of the plasmon description of the van der Waals density functional," Phys. Rev. B 89, 035412 (2014).</footnote>

<footnote n="79">⁷⁹J. Klimeš, D. R. Bowler, and A. Michaelides, "Chemical accuracy for the van der Waals density functional," J. Phys.: Condens. Matter 22, 022201 (2010).</footnote>

<footnote n="80">⁸⁰V. R. Cooper, "van der Waals density functional: An appropriate exchange functional," Phys. Rev. B 81, 161104 (2010).</footnote>

<footnote n="81">⁸¹J. Tomasi, B. Mennucci, and R. Cammi, "Quantum mechanical continuum solvation models," Chem. Rev. 105, 2999-3094 (2005).</footnote>

<footnote n="82">⁸²J. Tomasi and M. Persico, "Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent," Chem. Rev. 94, 2027-2094 (1994).</footnote>

<footnote n="83">⁸³J. Tomasi, "The physical model," in Continuum Solvation Models in Chemical Physics: From Theory to Applications, edited by B. Mennucci and R. Cammi (Wiley & Sons, Chichester, UK, 2007), Chap. 1.1, p. 16.</footnote>

<footnote n="84">⁸⁴A. Delgado, S. Corni, S. Pittalis, and C. A. Rozzi, "Modeling solvation effects in real-space and real-time within density functional approaches," J. Chem. Phys. 143, 144111 (2015).</footnote>

<footnote n="85">⁸⁵E. Cancés, B. Mennucci, and J. Tomasi, "A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics," J. Chem. Phys. 107, 3032 (1997).</footnote>

<footnote n="86">⁸⁶J.-L. Pascual-ahuir, E. Silla, and I. Tunon, "GEPOL: An improved description of molecular surfaces. III. A new algorithm for the computation of a solvent-excluding surface," J. Comput. Chem. 15, 1127-1138 (1994).</footnote>

<footnote n="87">⁸⁷M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su et al., "General atomic and molecular electronic structure system," J. Comput. Chem. 14, 1347-1363 (1993).</footnote>

<footnote n="88">⁸⁸J. Perdew, K. Burke, and M. Ernzerhof, "Perdew, Burke, and Ernzerhof reply," Phys. Rev. Lett. 80, 891 (1998).</footnote>

<footnote n="89">⁸⁹M. Caricato, F. Ingrosso, B. Mennucci, and J. Tomasi, "A time-dependent polarizable continuum model: Theory and application," J. Chem. Phys. 122, 154501 (2005).</footnote>

<footnote n="90">⁹⁰R. Cammi and J. Tomasi, "Nonequilibrium solvation theory for the polarizable continuum model: A new formulation at the SCF level with application to the case of the frequency-dependent linear electric response function," Int. J. Quantum Chem. 56, 465-474 (1995).</footnote>

<footnote n="91">⁹¹S. Corni, S. Pipolo, and R. Cammi, "Equation of motion for the solvent polarization apparent charges in the polarizable continuum model: Application to real-time TDDFT," J. Phys. Chem. A 119, 5405-5416 (2014).</footnote>
<footnote n="92">⁹²G. Gil, S. Pipolo, A. Delgado, C. A. Rozzi, and S. Corni, "Nonequilibrium solvent polarization effects in real-time electronic dynamics of solute molecules subject to time-dependent electric fields: A new feature of the polarizable continuum model," *J. Chem. Theory Comput.* **15**, 2306–2319 (2019).</footnote>

<footnote n="93">⁹³L. Onsager, "Electric moments of molecules in liquids," *J. Am. Chem. Soc.* **58**, 1486–1493 (1936).</footnote>

<footnote n="94">⁹⁴S. Corni, R. Cammi, B. Mennucci, and J. Tomasi, "Electronic excitation energies of molecules in solution within continuum solvation models: Investigating the discrepancy between state-specific and linear-response methods," *J. Chem. Phys.* **123**, 134512 (2005).</footnote>

<footnote n="95">⁹⁵P. Buczek, A. Ernst, and L. M. Sandratskii, "Different dimensionality trends in the Landau damping of magnons in iron, cobalt, and nickel: Time-dependent density functional study," *Phys. Rev. B* **84**, 174418 (2011).</footnote>

<footnote n="96">⁹⁶M. Niesert, "Ab initio calculations of spin-wave excitation spectra from time-dependent density-functional theory," Ph.D. dissertation (RWTH University, Aachen, 2011), record converted from VDB: 12.11.2012.</footnote>

<footnote n="97">⁹⁷B. Rousseau, A. Eiguren, and A. Bergara, "Efficient computation of magnon dispersions within time-dependent density functional theory using maximally localized wannier functions," *Phys. Rev. B* **85**, 054305 (2012).</footnote>

<footnote n="98">⁹⁸T. Gorni, I. Timrov, and S. Baroni, "Spin dynamics from time-dependent density functional perturbation theory," *Eur. Phys. J. B* **91**, 249 (2018).</footnote>

<footnote n="99">⁹⁹K. Cao, H. Lambert, P. G. Radaelli, and F. Giustino, "Ab initio calculation of spin fluctuation spectra using time-dependent density functional perturbation theory, plane waves, and pseudopotentials," *Phys. Rev. B* **97**, 024420 (2018).</footnote>

<footnote n="100">¹⁰⁰N. Singh, P. Elliott, T. Nautiyal, J. K. Dewhurst, and S. Sharma, "Adiabatic generalized gradient approximation kernel in time-dependent density functional theory," *Phys. Rev. B* **99**, 035151 (2019).</footnote>

<footnote n="101">¹⁰¹S. Y. Savrasov, "Linear response calculations of spin fluctuations," *Phys. Rev. Lett.* **81**, 2570–2573 (1998).</footnote>

<footnote n="102">¹⁰²N. Tancogne-Dejean, F. G. Eich, and A. Rubio, "Time-dependent magnons from first principles," *J. Chem. Theory Comput.* **16**, 1007–1017 (2020).</footnote>

<footnote n="103">¹⁰³L. M. Sandratskii, "Energy band structure calculations for crystals with spiral magnetic structure," *Phys. Status Solidi B* **136**, 167–180 (1986).</footnote>

<footnote n="104">¹⁰⁴H. A. Mook and D. M. Paul, "Neutron-scattering measurement of the spin-wave spectra for nickel," *Phys. Rev. Lett.* **54**, 227–229 (1985).</footnote>

<footnote n="105">¹⁰⁵U. De Giovannini, G. Brunetto, A. Castro, J. Walkenhorst, and A. Rubio, "Simulating pump-probe photoelectron and absorption spectroscopy on the attosecond timescale with time-dependent density functional theory," *ChemPhysChem* **14**, 1363–1376 (2013).</footnote>

<footnote n="106">¹⁰⁶J. Walkenhorst, U. De Giovannini, A. Castro, and A. Rubio, "Tailored pump-probe transient spectroscopy with time-dependent density-functional theory: Controlling absorption spectra," *Eur. Phys. J. B* **89**, 128 (2016).</footnote>

<footnote n="107">¹⁰⁷U. De Giovannini, "Pump-probe photoelectron spectra," in *Handbook of Materials Modeling Methods: Theory Modeling* (Springer, 2018), pp. 1–19.</footnote>

<footnote n="108">¹⁰⁸S. Sato, H. Hübener, U. De Giovannini, and A. Rubio, "Ab initio simulation of attosecond transient absorption spectroscopy in two-dimensional materials," *Appl. Sci.* **8**, 1777 (2018).</footnote>

<footnote n="109">¹⁰⁹L. D. Barron, *Molecular Light Scattering and Optical Activity*, 2nd ed. (Cambridge University Press, Cambridge, 2004).</footnote>

<footnote n="110">¹¹⁰F. R. Keßler and J. Metzdorf, "Landau level spectroscopy: Interband effects and Faraday rotation," in *Modern Problems in Condensed Matter Sciences*, edited by V. M. Agranovich and A. A. Maradudin (Elsevier Science Publishers, Amsterdam, 1991), Vol. 27.1, Chap. 11.</footnote>

<footnote n="111">¹¹¹H. Solheim, K. Ruud, S. Coriani, and P. Norman, "Complex polarization propagator calculations of magnetic circular dichroism spectra," *J. Chem. Phys.* **128**, 094103 (2008).</footnote>

<footnote n="112">¹¹²H. Solheim, K. Ruud, S. Coriani, and P. Norman, "The A and B terms of magnetic circular dichroism revisited," *J. Phys. Chem. A* **112**, 9615–9618 (2008).</footnote>

<footnote n="113">¹¹³K.-M. Lee, K. Yabana, and G. F. Bertsch, "Magnetic circular dichroism in real-time time-dependent density functional theory," *J. Chem. Phys.* **134**, 144106 (2011).</footnote>

<footnote n="114">¹¹⁴M. Seth, M. Krykunov, T. Ziegler, J. Autschbach, and A. Banerjee, "Application of magnetically perturbed time-dependent density functional theory to magnetic circular dichroism: Calculation of B terms," *J. Chem. Phys.* **128**, 144105 (2008).</footnote>

<footnote n="115">¹¹⁵M. Seth, M. Krykunov, T. Ziegler, and J. Autschbach, "Application of magnetically perturbed time-dependent density functional theory to magnetic circular dichroism. II. Calculation of A terms," *J. Chem. Phys.* **128**, 234102 (2008).</footnote>

<footnote n="116">¹¹⁶I. V. Lebedeva, D. A. Strubbe, I. V. Tokatly, and A. Rubio, "Orbital magnet optical response of periodic insulators from first principles," *npj Comput. Mater.* **5**, 32 (2019).</footnote>

<footnote n="117">¹¹⁷R. D. King-Smith and D. Vanderbilt, "Theory of polarization of crystalline solids," *Phys. Rev. B* **47**, 1651–1654 (1993).</footnote>

<footnote n="118">¹¹⁸D. Vanderbilt and R. D. King-Smith, "Electric polarization as a bulk quantity and its relation to surface charge," *Phys. Rev. B* **48**, 4442–4455 (1993).</footnote>

<footnote n="119">¹¹⁹R. Resta, "Macroscopic polarization in crystalline dielectrics: The geometric phase approach," *Rev. Mod. Phys.* **66**, 899–915 (1994).</footnote>

<footnote n="120">¹²⁰A. M. Essin, A. M. Turner, J. E. Moore, and D. Vanderbilt, "Orbital magnetoelectric coupling in band insulators," *Phys. Rev. B* **81**, 205104 (2010).</footnote>

<footnote n="121">¹²¹K.-T. Chen and P. A. Lee, "Unified formalism for calculating polarization, magnetization, and more in a periodic insulator," *Phys. Rev. B* **84**, 205137 (2011).</footnote>

<footnote n="122">¹²²X. Gonze and J. W. Zwanziger, "Density-operator theory of orbital magnetic susceptibility in periodic insulators," *Phys. Rev. B* **84**, 064445 (2011).</footnote>

<footnote n="123">¹²³X. Andrade, S. Botti, M. A. L. Marques, and A. Rubio, "Time-dependent density functional theory scheme for efficient calculations of dynamic (hyper)polarizabilities," *J. Chem. Phys.* **126**, 184106 (2007).</footnote>

<footnote n="124">¹²⁴D. A. Strubbe, "Optical and transport properties of organic molecules: Methods and applications," Ph.D. thesis, University of California, Berkeley, USA, 2012.</footnote>

<footnote n="125">¹²⁵D. A. Strubbe, L. Lehtovaara, A. Rubio, M. A. L. Marques, and S. G. Louie, "Response functions in TDDFT: Concepts and implementation," in *Fundamentals of Time-Dependent Density Functional Theory* (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012), pp. 139–166.</footnote>

<footnote n="126">¹²⁶M. Lazzeri and F. Mauri, "High-order density-matrix perturbation theory," *Phys. Rev. B* **68**, 161101 (2003).</footnote>

<footnote n="127">¹²⁷X. Gonze and J.-P. Vigneron, "Density-functional approach to nonlinear-response coefficients of solids," *Phys. Rev. B* **39**, 13120–13128 (1989).</footnote>

<footnote n="128">¹²⁸A. D. Corso, F. Mauri, and A. Rubio, "Density-functional theory of the nonlinear optical susceptibility: Application to cubic semiconductors," *Phys. Rev. B* **53**, 15638–15642 (1996).</footnote>

<footnote n="129">¹²⁹J. A. Berger, "Fully parameter-free calculation of optical spectra for insulators, semiconductors, and metals from a simple polarization functional," *Phys. Rev. Lett.* **115**, 137402 (2015).</footnote>

<footnote n="130">¹³⁰S. Albrecht, L. Reining, R. Del Sole, and G. Onida, "Ab initio calculation of excitonic effects in the optical spectra of semiconductors," *Phys. Rev. Lett.* **80**, 4510–4513 (1998).</footnote>

<footnote n="131">¹³¹S. Botti, F. Sottile, N. Vast, V. Olevano, L. Reining, H.-C. Weissker, A. Rubio, G. Onida, R. Del Sole, and R. W. Godby, "Long-range contribution to the exchange-correlation kernel of time-dependent density functional theory," *Phys. Rev. B* **69**, 155112 (2004).</footnote>

<footnote n="132">¹³²P. Lautenschlager, M. Garriga, L. Viñà, and M. Cardona, "Temperature dependence of the dielectric function and interband critical points in silicon," *Phys. Rev. B* **36**, 4821–4830 (1987).</footnote>

<footnote n="133">¹³³U. De Giovannini and A. Castro, "Real-time and real-space time-dependent density-functional theory approach to attosecond dynamics," in *Attosecond Molecular Dynamics* (Royal Society of Chemistry, Cambridge, 2018), pp. 424–461.</footnote>

<footnote n="134">¹³⁴U. De Giovannini, H. Hübener, and A. Rubio, "A first-principles time-dependent density functional theory framework for spin and time-resolved angular-resolved photoelectron spectroscopy in periodic systems," *J. Chem. Theory Comput.* **13**, 265–273 (2017).</footnote>

<footnote n="135">¹³⁵L. Tao and A. Scrinzi, "Photo-electron momentum spectra from minimal volumes: The time-dependent surface flux method," *New J. Phys.* **14**, 013021 (2012).</footnote>

<footnote n="136">¹³⁶P. Wopperer, U. De Giovannini, and A. Rubio, "Efficient and accurate modeling of electron photoemission in nanostructures with TDDFT," *Eur. Phys. J. B* **90**, 1307 (2017).</footnote>
<footnote n="92">⁹²G. Gil, S. Pipolo, A. Delgado, C. A. Rozzi, and S. Corni, "Nonequilibrium solvent polarization effects in real-time electronic dynamics of solute molecules subject to time-dependent electric fields: A new feature of the polarizable continuum model," *J. Chem. Theory Comput.* **15**, 2306–2319 (2019).</footnote>

<footnote n="93">⁹³L. Onsager, "Electric moments of molecules in liquids," *J. Am. Chem. Soc.* **58**, 1486–1493 (1936).</footnote>

<footnote n="94">⁹⁴S. Corni, R. Cammi, B. Mennucci, and J. Tomasi, "Electronic excitation energies of molecules in solution within continuum solvation models: Investigating the discrepancy between state-specific and linear-response methods," *J. Chem. Phys.* **123**, 134512 (2005).</footnote>

<footnote n="95">⁹⁵P. Buczek, A. Ernst, and L. M. Sandratskii, "Different dimensionality trends in the Landau damping of magnons in iron, cobalt, and nickel: Time-dependent density functional study," *Phys. Rev. B* **84**, 174418 (2011).</footnote>

<footnote n="96">⁹⁶M. Niesert, "Ab initio calculations of spin-wave excitation spectra from time-dependent density-functional theory," Ph.D. dissertation (RWTH University, Aachen, 2011), record converted from VDB: 12.11.2012.</footnote>

<footnote n="97">⁹⁷B. Rousseau, A. Eiguren, and A. Bergara, "Efficient computation of magnon dispersions within time-dependent density functional theory using maximally localized wannier functions," *Phys. Rev. B* **85**, 054305 (2012).</footnote>

<footnote n="98">⁹⁸T. Gorni, I. Timrov, and S. Baroni, "Spin dynamics from time-dependent density functional perturbation theory," *Eur. Phys. J. B* **91**, 249 (2018).</footnote>

<footnote n="99">⁹⁹K. Cao, H. Lambert, P. G. Radaelli, and F. Giustino, "Ab initio calculation of spin fluctuation spectra using time-dependent density functional perturbation theory, plane waves, and pseudopotentials," *Phys. Rev. B* **97**, 024420 (2018).</footnote>

<footnote n="100">¹⁰⁰N. Singh, P. Elliott, T. Nautiyal, J. K. Dewhurst, and S. Sharma, "Adiabatic generalized gradient approximation kernel in time-dependent density functional theory," *Phys. Rev. B* **99**, 035151 (2019).</footnote>

<footnote n="101">¹⁰¹S. Y. Savrasov, "Linear response calculations of spin fluctuations," *Phys. Rev. Lett.* **81**, 2570–2573 (1998).</footnote>

<footnote n="102">¹⁰²N. Tancogne-Dejean, F. G. Eich, and A. Rubio, "Time-dependent magnons from first principles," *J. Chem. Theory Comput.* **16**, 1007–1017 (2020).</footnote>

<footnote n="103">¹⁰³L. M. Sandratskii, "Energy band structure calculations for crystals with spiral magnetic structure," *Phys. Status Solidi B* **136**, 167–180 (1986).</footnote>

<footnote n="104">¹⁰⁴H. A. Mook and D. M. Paul, "Neutron-scattering measurement of the spin-wave spectra for nickel," *Phys. Rev. Lett.* **54**, 227–229 (1985).</footnote>

<footnote n="105">¹⁰⁵U. De Giovannini, G. Brunetto, A. Castro, J. Walkenhorst, and A. Rubio, "Simulating pump-probe photoelectron and absorption spectroscopy on the attosecond timescale with time-dependent density functional theory," *ChemPhysChem* **14**, 1363–1376 (2013).</footnote>

<footnote n="106">¹⁰⁶J. Walkenhorst, U. De Giovannini, A. Castro, and A. Rubio, "Tailored pump-probe transient spectroscopy with time-dependent density-functional theory: Controlling absorption spectra," *Eur. Phys. J. B* **89**, 128 (2016).</footnote>

<footnote n="107">¹⁰⁷U. De Giovannini, "Pump-probe photoelectron spectra," in *Handbook of Materials Modeling Methods: Theory Modeling* (Springer, 2018), pp. 1–19.</footnote>

<footnote n="108">¹⁰⁸S. Sato, H. Hübener, U. De Giovannini, and A. Rubio, "Ab initio simulation of attosecond transient absorption spectroscopy in two-dimensional materials," *Appl. Sci.* **8**, 1777 (2018).</footnote>

<footnote n="109">¹⁰⁹L. D. Barron, *Molecular Light Scattering and Optical Activity*, 2nd ed. (Cambridge University Press, Cambridge, 2004).</footnote>

<footnote n="110">¹¹⁰F. R. Keßler and J. Metzdorf, "Landau level spectroscopy: Interband effects and Faraday rotation," in *Modern Problems in Condensed Matter Sciences*, edited by V. M. Agranovich and A. A. Maradudin (Elsevier Science Publishers, Amsterdam, 1991), Vol. 27.1, Chap. 11.</footnote>

<footnote n="111">¹¹¹H. Solheim, K. Ruud, S. Coriani, and P. Norman, "Complex polarization propagator calculations of magnetic circular dichroism spectra," *J. Chem. Phys.* **128**, 094103 (2008).</footnote>

<footnote n="112">¹¹²H. Solheim, K. Ruud, S. Coriani, and P. Norman, "The A and B terms of magnetic circular dichroism revisited," *J. Phys. Chem. A* **112**, 9615–9618 (2008).</footnote>

<footnote n="113">¹¹³K.-M. Lee, K. Yabana, and G. F. Bertsch, "Magnetic circular dichroism in real-time time-dependent density functional theory," *J. Chem. Phys.* **134**, 144106 (2011).</footnote>

<footnote n="114">¹¹⁴M. Seth, M. Krykunov, T. Ziegler, J. Autschbach, and A. Banerjee, "Application of magnetically perturbed time-dependent density functional theory to magnetic circular dichroism: Calculation of B terms," *J. Chem. Phys.* **128**, 144105 (2008).</footnote>

<footnote n="115">¹¹⁵M. Seth, M. Krykunov, T. Ziegler, and J. Autschbach, "Application of magnetically perturbed time-dependent density functional theory to magnetic circular dichroism. II. Calculation of A terms," *J. Chem. Phys.* **128**, 234102 (2008).</footnote>

<footnote n="116">¹¹⁶I. V. Lebedeva, D. A. Strubbe, I. V. Tokatly, and A. Rubio, "Orbital magneto-optical response of periodic insulators from first principles," *npj Comput. Mater.* **5**, 32 (2019).</footnote>

<footnote n="117">¹¹⁷R. D. King-Smith and D. Vanderbilt, "Theory of polarization of crystalline solids," *Phys. Rev. B* **47**, 1651–1654 (1993).</footnote>

<footnote n="118">¹¹⁸D. Vanderbilt and R. D. King-Smith, "Electric polarization as a bulk quantity and its relation to surface charge," *Phys. Rev. B* **48**, 4442–4455 (1993).</footnote>

<footnote n="119">¹¹⁹R. Resta, "Macroscopic polarization in crystalline dielectrics: The geometric phase approach," *Rev. Mod. Phys.* **66**, 899–915 (1994).</footnote>

<footnote n="120">¹²⁰A. M. Essin, A. M. Turner, J. E. Moore, and D. Vanderbilt, "Orbital magneto-electric coupling in band insulators," *Phys. Rev. B* **81**, 205104 (2010).</footnote>

<footnote n="121">¹²¹K.-T. Chen and P. A. Lee, "Unified formalism for calculating polarization, magnetization, and more in a periodic insulator," *Phys. Rev. B* **84**, 205137 (2011).</footnote>

<footnote n="122">¹²²X. Gonze and J. W. Zwanziger, "Density-operator theory of orbital magnetic susceptibility in periodic insulators," *Phys. Rev. B* **84**, 064445 (2011).</footnote>

<footnote n="123">¹²³X. Andrade, S. Botti, M. A. L. Marques, and A. Rubio, "Time-dependent density functional theory scheme for efficient calculations of dynamic (hyper)polarizabilities," *J. Chem. Phys.* **126**, 184106 (2007).</footnote>

<footnote n="124">¹²⁴D. A. Strubbe, "Optical and transport properties of organic molecules: Methods and applications," Ph.D. thesis, University of California, Berkeley, USA, 2012.</footnote>

<footnote n="125">¹²⁵D. A. Strubbe, L. Lehtovaara, A. Rubio, M. A. L. Marques, and S. G. Louie, "Response functions in TDDFT: Concepts and implementation," in *Fundamentals of Time-Dependent Density Functional Theory* (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012), pp. 139–166.</footnote>

<footnote n="126">¹²⁶M. Lazzeri and F. Mauri, "High-order density-matrix perturbation theory," *Phys. Rev. B* **68**, 161101 (2003).</footnote>

<footnote n="127">¹²⁷X. Gonze and J.-P. Vigneron, "Density-functional approach to nonlinear-response coefficients of solids," *Phys. Rev. B* **39**, 13120–13128 (1989).</footnote>

<footnote n="128">¹²⁸A. D. Corso, F. Mauri, and A. Rubio, "Density-functional theory of the nonlinear optical susceptibility: Application to cubic semiconductors," *Phys. Rev. B* **53**, 15638–15642 (1996).</footnote>

<footnote n="129">¹²⁹J. A. Berger, "Fully parameter-free calculation of optical spectra for insulators, semiconductors, and metals from a simple polarization functional," *Phys. Rev. Lett.* **115**, 137402 (2015).</footnote>

<footnote n="130">¹³⁰S. Albrecht, L. Reining, R. Del Sole, and G. Onida, "Ab initio calculation of excitonic effects in the optical spectra of semiconductors," *Phys. Rev. Lett.* **80**, 4510–4513 (1998).</footnote>

<footnote n="131">¹³¹S. Botti, F. Sottile, N. Vast, V. Olevano, L. Reining, H.-C. Weissker, A. Rubio, G. Onida, R. Del Sole, and R. W. Godby, "Long-range contribution to the exchange-correlation kernel of time-dependent density functional theory," *Phys. Rev. B* **69**, 155112 (2004).</footnote>

<footnote n="132">¹³²P. Lautenschlager, M. Garriga, L. Viñà, and M. Cardona, "Temperature dependence of the dielectric function and interband critical points in silicon," *Phys. Rev. B* **36**, 4821–4830 (1987).</footnote>

<footnote n="133">¹³³U. De Giovannini and A. Castro, "Real-time and real-space time-dependent density-functional theory approach to attosecond dynamics," in *Attosecond Molecular Dynamics* (Royal Society of Chemistry, Cambridge, 2018), pp. 424–461.</footnote>

<footnote n="134">¹³⁴U. De Giovannini, H. Hübener, and A. Rubio, "A first-principles time-dependent density functional theory framework for spin and time-resolved angular-resolved photoelectron spectroscopy in periodic systems," *J. Chem. Theory Comput.* **13**, 265–273 (2017).</footnote>

<footnote n="135">¹³⁵L. Tao and A. Scrinzi, "Photo-electron momentum spectra from minimal volumes: The time-dependent surface flux method," *New J. Phys.* **14**, 013021 (2012).</footnote>

<footnote n="136">¹³⁶P. Wopperer, U. De Giovannini, and A. Rubio, "Efficient and accurate modeling of electron photoemission in nanostructures with TDDFT," *Eur. Phys. J. B* **90**, 1307 (2017).</footnote>