Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems

It is a general challenge in the electronic structure community to develop accurate and efficient methods for modeling materials of ever increasing complexity in order to predict their properties. In this respect, time-dependent density functional theory (TDDFT) and related methods have become a natural choice for modeling systems in and out of their equilibrium.

During the last few years, novel directions of research emerged in the fields of chemistry, physics, and materials science that required the development of novel simulation tools needed to face the new challenges posed by those experimental advances and emerging phenomena. One can cite some examples among these new fields of research: the strong coupling between light and matter (including materials embedded in cavities), new states of matter (hidden phases and topological solids and molecules), and strong field dynamics in periodic systems. Whereas the strong-field dynamics of atoms and molecules is now well understood, strong-field dynamics in solids is an active field of research. Real-time (rt) TDDFT^1,2 represents a natural tool to study highly non-linear phenomena in solids and low-dimensional materials, and the development of efficient numerical methods to perform real-time TDDFT in such periodic or semi-periodic systems is crucial to explore this new phenomenon (including the quantum nature of light and phonons). Indeed, it allows us to describe highly nonlinear processes without having to resort to perturbation theory on top of equilibrium DFT calculations.

Recent years have seen tremendous experimental progress in the field of strong light–matter interactions,^3–5 where the strong coupling of light to chemical systems, quantum materials, or nanoplasmonic systems, among others, has been demonstrated. In this regime, light and matter meet on the same footing and the electron–photon interaction has to be explicitly considered.^4–9 A novel theoretical approach accelerating this field is quantum-electrodynamical density-functional theory (QEDFT),^4,5,10–14 which complements TDDFT with the photonic degrees of freedom and provides reliable and predictive simulations in this emerging field of research.

In this paper, we explore the recent advances in the Octopus project.^15–19 We focus our attention on the recently added features and, particularly, on the ones that have not been described in the previous papers.^17–19 These new features include the implementation of new levels of self-consistent and microscopic couplings of light and matter, the treatment of solvent effects through the polarizable continuum model (PCM), the implementation of various methods to treat van der Waals (vdW) interactions, new methods to calculate magnons, conductivities, and photoelectron spectroscopy from real-time TDDFT, and the calculation of orbital magneto-optical responses. Advances in numerical algorithms and methods, such as new propagators and the use of iterative eigensolvers in the context of reduced density matrix (RDM) functional theory, are also discussed. Finally, recent improvements in the treatment of periodic systems, as well as more technical code improvements, are also presented. For a more detailed explanation about how to use these newly introduced features in practice, we refer the readers to the Octopus webpage,²⁰ as new tutorials and examples are regularly being added to it.

This paper is organized as follows: First, we present in Secs. II–XII the new implementations of physical theories and algorithms that allow us to deal with non-equilibrium phenomena in materials and nanostructures. This is followed by Secs. XIII–XVI dealing with technical developments that are fundamental in improving the code performance and the stability of the algorithms. Finally, we draw our conclusions in Sec. XVII. Unless otherwise stated, atomic units are used throughout the paper.

In most cases when light–matter interactions are considered, a decoupling of light and matter is performed at the outset. Either the electromagnetic fields are prescribed and then the properties of the matter subsystem are determined, as frequently found in, e.g., quantum chemistry or solid-state physics, or the properties of matter are prescribed and then the properties of the photon subsystem are determined, as done in, e.g., quantum optics or photonics.^12,21

The microscopic interaction of light and matter in Octopus has followed this decoupling strategy and has been treated so far only in the forward-coupling direction. In this approximation, an external classical laser pulse or kick is prescribed and the response of the system is computed by evolving the Kohn–Sham orbitals.²² The back-action of the matter subsystem on the laser pulse, the subsequent effect of this modified pulse on the matter, and so on are ignored. This forward-coupling approximation is highly accurate when the total current generated in the system is comparably small, such as in atoms or small molecules. Small refers here to currents that induce transverse electric fields that are about an order of magnitude smaller than the longitudinal electric field that originates from the Hartree potential of the system. This has been exploited in Octopus and other codes, and many different types of spectroscopies were computed successfully in the past using this approach.

In contrast, in classical electromagnetic modeling, the opposite view is taken. Here, the material properties are prescribed and the resulting electromagnetic fields are computed. In practice, the material properties are routinely approximated by a local continuum model or by the dielectric function of the system (such as the Debye, Lorentz, or Drude models) and then Maxwell’s equations are solved for linear dielectric media arranged in appropriate geometries.²³ 

It is clear that both perspectives, a focus on matter dynamics alone or a focus on electromagnetic field dynamics alone, break down when the total induced currents become large and when electromagnetic near-field effects on the scale of the material system are not negligible anymore. Prime examples for such cases are nanoplasmonic systems, surface plasmon-polaritons, or tip-enhanced spectroscopies. In these cases, the back-action of the material response on the system itself has to be taken into account, leading to screening and retardation effects. The proper theoretical framework, which encompasses all these effects, is quantum electrodynamics. Starting with a generalized Pauli-Fierz field theory for the combined system of electrons, nuclei, and photons, we have recently derived different levels of self-consistent and microscopic couplings of light and matter. In the classical limit, this results in a coupled set of Ehrenfest–Maxwell–Pauli–Kohn–Sham equations.²¹ To implement these equations, we added a Maxwell solver to the Octopus code, which we couple self-consistently to the dynamics of the electrons and nuclei. In the following, we briefly summarize the basic ingredients for this implementation and show an example of self-consistent light–matter interactions for a nanoplasmonic system. Further details of the implementation and nano-optical applications can be found in Ref. 21.

Since over the years, Octopus has been optimized heavily to solve time-dependent Schrödinger and Kohn–Sham equations, we have exploited the fact that Maxwell’s equations can be formulated in Schrödinger form²⁴ to benefit from the efficient time-evolution in the code. This reformulation is based on the Riemann–Silberstein vector,²⁵ which is a combination of the electric E(r, t) and magnetic field B(r, t),

The sign of the imaginary part of the Riemann–Silberstein vector corresponds to different helicities. The reformulation of Maxwell’s equations in Schrödinger form is purely algebraic and starts out with the microscopic Maxwell’s equations

where E and B are the classical electric and magnetic fields, respectively, ρ and J are the charge and current densities, ϵ₀ and μ₀ are the vacuum permittivity and permeability, respectively, and $c=ϵ0μ0−1/2$ is the speed of light. Using the Riemann–Silberstein vector, the electric and magnetic Gauss laws may now be combined,

and likewise, the Faraday and Ampere laws can be combined into one evolution equation for the Riemann–Silberstein vector,

Here, S = (S_x, S_y, S_z) denotes a vector of spin one matrices

which are analogous to the Pauli matrices and show the spin-one character of the photon. Having cast Maxwell’s equations as an inhomogeneous Schrödinger equation, it is now straightforward to use the time-evolution algorithms in Octopus to time-evolve the Riemann–Silberstein vector. The only difference to the matter propagation is that we are now dealing with the “Maxwell Hamiltonian”

which acts upon six orbitals of the Riemann–Silberstein vector corresponding to the three components of the electric and magnetic field vectors. As in the matter case in Octopus, the discretization of the gradient in the Maxwell Hamiltonian is performed with finite-difference stencils, and the domain parallelization of Octopus can be used seamlessly for the Maxwell case as well. The difference to finite-difference time-domain (FDTD) codes based on the Yee algorithm is that we do not employ two shifted grids for the electric and magnetic fields,²³ but rather a single grid for the Riemann–Silberstein vector. This simplifies the coupling to matter and allows us to use higher-order finite-difference discretizations for the gradient. Since the spatial discretization is connected to the temporal discretization through the Courant condition, this, in turn, allows us to use larger time steps. Furthermore, from our experience, a unified grid also improves the stability compared to FDTD.

Instead of using the constitutive relations, we couple Maxwell’s equations directly to the microscopic current density of the matter subsystem, consisting of the usual paramagnetic current term, the diamagnetic current term, and the magnetization current term. For the coupling of the electromagnetic fields to the matter subsystem, we are relying on the Power–Zienau–Woolley transformation,^26,27 which leads to a multipole expansion. We have implemented the first two orders of this expansion: the dipole approximation in lowest order and electrical-quadrupole and magnetic-dipole coupling in the next order. In addition, we are currently working on implementing the full minimal coupling with a full position dependent vector potential.

The time evolution of the Kohn–Sham orbitals and the Maxwell fields is performed side-by-side and the two subsystems are coupled self-consistently in each time step, as illustrated in Fig. 1. To propagate different subsystems with different Hamiltonians and different sets of orbitals simultaneously, we implemented a multi-system framework in Octopus (for more details, refer Subsection XVI C).

Similar to the matter propagation, also in the Maxwell case, outgoing waves that reach the boundary of the simulation box have to be absorbed to avoid artificial reflections and backscattering. Our first attempt was to use also the mask functions that are used for the matter propagation in Octopus. However, in the electromagnetic case, the mask absorption of outgoing waves turned out to be not efficient enough so that we implemented a perfectly matched layer (PML)²³ for the Maxwell propagation.

When considering incoming electromagnetic fields with optical wavelengths, the coupling to atomistic or nano-scale systems leads to a multi-scale problem. The optical wavelength of the radiation is in this case much larger than the de Broglie wavelength of the matter. Likewise, the electromagnetic waves are traveling with the speed of light, which requires sub-attosecond time steps. We have, therefore, implemented different multi-scale couplings in space and time. For example, the Maxwell simulation box can be on the same scale as the matter box. In this case, the electromagnetic waves are represented as incoming analytical time-dependent boundary conditions and are propagated numerically inside the simulation box. Alternatively, the Maxwell simulation box can be much larger than the matter box to fully encompass laser pulses with optical wavelengths. In this case, prolongations and interpolations have to be used similarly to multi-grid methods. Since the electronic and nuclear motion is much slower than the time-evolution of the electromagnetic waves, we also implemented a multi-scale approach for the real-time propagation. The Riemann–Silberstein vector is propagated with frozen electronic current from the last point of interaction for many intermediate time steps before a coupling to the matter subsystem takes place. The number of intermediate steps is a convergence parameter and depends on the physical situation at hand.

Since we now include the description of classical electromagnetic fields explicitly in our real-time simulations, we have directly access to the outgoing electromagnetic radiation. This allows us to define electromagnetic detectors at the box boundaries, which accumulate the outgoing electromagnetic waves. We implemented such electromagnetic detectors in Octopus, and this allows us to run simulations in close analogy with experiments and to directly observe the outgoing radiation. For example, it is no longer needed to Fourier transform the time signal of the matter dipole to get optical spectra, but we rather have access to the spectrum directly on the Maxwell grid.

As an example of a coupled Ehrenfest–Maxwell–Pauli–Kohn–Sham propagation with Octopus, we selected a nano-optical application. We consider in this example two almost spherical sodium nanoparticles with 297 sodium atoms each, which are arranged in a dimer configuration. This system interacts with an incoming laser pulse that excites either the internal dipole or quadrupole plasmon motion of the dimer. In Fig. 2, we show the resulting electromagnetic field enhancements for different levels of light–matter coupling. In panels (a) and (b), we show the temporal profile of the incoming laser pulse and the resulting current density at the center point between the two nanoparticles. The field enhancement can be seen in panels (c)–(e). Including a self-consistent back reaction in the light–matter coupling, the field enhancement is reduced at the center point between the two nanoparticles, as shown in panel (c), while far away from the dimer, as shown in panels (d) and (e), the field enhancement is larger than in the forward coupled case. Furthermore, frequency shifts can also be observed, which are more pronounced in the near-field. We found that the field enhancement is also sensitive to the coupling terms of the multipole expansion,which are included, and that the quantitative difference of switching from local-density approximation (LDA) to the Perdew-Burke-Ernzerhof (PBE) functional for the exchange–correlation functional is in this case smaller than including the back reaction in the light–matter coupling (cf. Ref. 21).

To conclude, with our new efficient implementation for coupled Ehrenfest–Maxwell–Pauli–Kohn–Sham equations in Octopus, we have now a very versatile tool that allows us to compute fully self-consistent forward–backward light–matter coupling in real-time and real-space for a vast set of applications and in close analogy to experiments. As the coupling to Maxwell’s equations of motion (EOM) represents the classical limit of the light–matter interaction, this development leads to the classical limit of QEDFT.

The nascent field of strong light–matter interaction expanded over the past decades from small atomic structures to chemistry³ and materials science.²⁸ This development necessitates predictive first-principles methods capable to describe light and matter on the same footing. We introduced for the first time, a time-dependent density functional theory for quantum electrodynamics (QEDFT)^4,10 to treat ab initio weak and strong light–matter interactions and its applications to chemistry and materials, which provides a unique framework to explore, predict, and control new states of matter out of equilibrium. This generalization of the time-dependent density-functional method allows us, for the first time, to explore the effects of dressing electronic states with photons while retaining the electronic properties of real materials.

A general non-relativistic Hamiltonian for light–matter systems treating N interacting electrons coupled to N_p photonic modes in the case of the so-called length-gauge, and after employing the long-wavelength (dipole) approximation,¹³ reads as follows:

Here, the first two terms on the right-hand side correspond to the usual electronic many-body Hamiltonian, while the last term describes the photon modes, which is characterized for each photon mode α by its elongation q_α, frequency ω_α, and electron–photon coupling strength vector λ_α that includes the polarization of the photon mode and introduces the coupling to the total dipole $∑k=1Nrk$ of the electronic system. For a more detailed discussion of the different regimes of weak/strong/ultra-strong light–matter coupling, we refer to Ref. 6. The external variable for the photon system is the time-dependent current $jext(α)(t)$⁠.

QEDFT^4,10 is structurally similar to time-dependent density-functional theory in that it is based on a one-to-one correspondence between internal and external variables. If now photons are also considered, the set of internal variables has to be expanded. In the framework of Eq. (8), the internal variables become the density n(r, t) and the mode-resolved contributions to the electric displacement field q_α(t). By introducing and exploiting the bijective mapping of these internal and the external variables [$v$_ext(r, t) and j_ext(t)], the auxiliary Kohn–Sham system is described by the electronic Kohn–Sham equations as well as Maxwell’s equations,⁵ leading to no exchange–correlation contribution in the photon subsystem [$jxc(α)(t)=0$]. This reformulation subsumes the “quantumness” of the light–matter interaction solely into the local exchange–correlation potential that now features a component due to the electron–photon interactions, in addition to the part due to electron–electron interaction, i.e., $vxcσ(r,t)=vxcσee(r,t)+vxcσep(r,t)$⁠.^4,10,12 Extending the coupled Maxwell–Kohn–Sham equations to consider quantum photons is, thus, condensed into the calculation of the local exchange–correlation potential.

In practice, QEDFT requires the construction of an additional exchange–correlation potential that describes the electron–photon interaction. First attempts to construct this potential were based on many-body perturbation theory via the Sham–Schlüter equation²⁹ within the exact exchange (EXX) approximation^13,30 by utilizing the optimized-effective potential (OEP) method.³¹ 

The framework of the optimized effective potential (OEP) equation for electron–photon systems has been introduced in Ref. 30 for the static and time-dependent regimes. Solving the time-dependent OEP equation is unfortunately computationally challenging. Applications to electronic systems^32–34 have been realized so far only for small and low-dimensional systems while their utilization for interacting electron–photon systems remained restricted to models.^30,35,36 This section provides a brief introduction to the implementation that allows us to solve the computationally more tractable static OEP equation.

As introduced in Ref. 30, the OEP photon energy depends on both the occupied and unoccupied electronic orbitals. Alternatively, the OEP photon energy can also be formulated using occupied orbitals and orbital shifts only.¹³ The full energy expression is given by

where $Exc(ee)$ describes the electronic exchange–correlation energy and $Ex(α)$ describes the exchange energy due to the interaction of the electrons with the photon mode α. Avoiding unoccupied orbitals is computationally much more favorable for larger systems and the photonic induced exchange energy can be correspondingly expressed as a sum over the N_σ occupied orbitals of spin channel σ,

where ω_α describes the αth mode of the electromagnetic field and $d^α=λα⋅r$ describes the dipole operator and the electron–photon coupling strength. We can now reformulate the problem in terms of two electron–photon orbital shifts. The Kohn–Sham orbitals φ_iσ contribute to both electron–photon orbital shifts $Φiσ,α(1)$ and $Φiσ,α(2)$ that can be calculated using the Sternheimer equations.¹³ The first electron–photon orbital shift can be obtained explicitly by the solution of a linear Sternheimer equation

with the matrix element $dijσ(α)=φiσ|d^α|φjσ$⁠. In addition, the second electron–photon orbital shift $Φiσ,α(2)(r)$ can be defined explicitly as follows:

From Eq. (10), we can now deduce the potential using

After these reformulations, we find for the final OEP equation including electron–electron effects as well as electron–photon effects

where the inhomogeneity Λ_iσ(r) is given by

In Eq. (14), we defined a third orbital shift, the exchange–correlation orbital shift, which will be used to obtain the corresponding exchange–correlation potential and is defined along the lines of the orbital shift usually used in OEP calculations.³¹ We can also obtain $ψiσ*(r)$ using a Sternheimer equation

where $Miσ*(r)$ now consist of the electron–photon orbital shifts and the Kohn–Sham orbitals, as described in Ref. 13. Accordingly, we define this quantity as

where the effects of the electron–electron interaction are included in the quantity u_xiσ(r). For instance, in exchange-only calculations, this quantity is defined as $uxiσ(r)=1φiσ*(r)δEx(ee)[{φjτ}]δφiσ(r)$⁠, where $Ex(ee)$ is the usual Fock exchange energy.

Equation (15) has to be solved self-consistently with Eq. (11). By this reformulation, we replaced the problem of calculating the OEP equation using all unoccupied states by a problem of solving N_p+1 Sternheimer equations that only involve the occupied orbitals. In this way, the formulation of the problem becomes similar to the one of Ref. 37 for electrons only.

For practical implementation, we reformulate the OEP equation in the following form, as it is commonly done to construct the electronic OEP:³¹ 

and update the potential with

The quantity S_σ(r) becomes a measure for convergence, since it is vanishing at the solution point [compare Eqs. (17) and (14)]. For the function c(r), we have different possibilities, such as using a constant or using the inverse of the electron density, as used in Ref. 37. Other methods are also possible, such as the Barzilai–Borwein (BB) method.³⁸ We found stable algorithms when using a constant and the Barzilai–Borwein method.

While we will show the computational feasibility of this approach in the following, often a simplified solution is beneficial as a starting point for the self-consistency procedure. Such a simplified approximation can be deduced by reformulating Eq. (15) into the following equivalent form:

and subsequently assuming $(ĥsσ(r)−ϵiσ)ψiσ*(r)=0$ to start the iterative process. In situations where Λ_iσ(r) = 0, such as pure electronic exact exchange, this approximation is exact for a single electron and is referred to as the Krieger–Li–Iafrate (KLI) approximation.^{13,37,39–41} By multiplying Eq. (19) by |φ_jσ(r)|² and integrating over the spatial coordinates, we arrive at a linear equation that, in turn, can be solved for the approximate $v$_xσ. This approximation scheme has proven to often deliver sufficiently accurate results for electronic structure calculations with significantly lower computational effort and reliable stability. As this approximation can be seen as a diagonal approximation to the response function, it unavoidably fails in accurately describing polarization features. In the context of light–matter correlated ground-states, this leads to a slight unbalance when approximating components including photonic excitations [introduced by Φ⁽¹⁾] and self-polarization interaction [introduced by Φ⁽²⁾].^13,30,42,43 This results in a violation of translational invariance and introduces an artificial dependence on the permanent dipole. When performing KLI calculations including light–matter interaction, we thus suggest moving the set of coordinates into the electronic center-of-charge instead of the center-of-mass frame. To reduce the effect of this dependence during a self-consistent calculation, the optional input parameter KLIpt_coc has been introduced in the code. When activated, this option defines the dipole operator $d^α$ with respect to the electronic center-of-charge, thus improving the stability of the algorithm.

Finally, in Fig. 3, we show the capabilities of the new implementation, where we calculate two sodium dimers in the weak coupling regime under light–matter coupling. Figure 3(a) shows the electron density, and Fig. 3(b) shows the comparison of OEP and KLI results. As shown in Ref. 13, in the weak-coupling regime, the KLI is close to the OEP result. In Fig. 3(c), we show the convergence behavior when using a constant in Eq. (18) and when using the Barzilai–Borwein method.³⁸ 

For a more detailed analysis on how the OEP equation performs also in the strong and the ultra-strong coupling limit, we refer to Refs. 13 and 30, where we find accurate results for few-level systems and asymmetric external potentials up to high electron–photon coupling strength.

Extensions of QEDFT to the regime of vibrational strong coupling,⁴⁴ the linear-response regime,⁴⁵ and multitrajectory methods that capture quantum fluctuations⁴⁶ are currently work in progress and will further strengthen the capabilities of the Octopus code for the real-space description of strongly coupled light–matter systems. To describe the effects in the ultra-strong coupling regime, one can use an alternative method that is presented in Sec. IV.

The accurate description of the (ultra-)strong coupling regime of light–matter systems is a formidable task. In many cases, the known functionals for QEDFT (see Sec. III) are inaccurate, and for complex electronic systems, typical few-level approximations become unreliable.^42,47 In this section, we present the Octopus implementation of an alternative real-space ab initio method for coupled light–matter systems. Dressed Reduced Density Matrix Functional Theory (dressed RDMFT)⁴⁸ extends standard electronic RDMFT to coupled light–matter systems similarly to how QEDFT extends DFT. First tests on simple model systems suggest that the dressed RDMFT remains accurate from the weak to the ultra-strong coupling regime. A proper introduction of the theory, examples, and convergence studies can be found in Ref. 49.

This approach allows for the description of an interacting N-electron system coupled to one photonic mode within the dipole approximation. The respective Hamiltonian is given in Eq. (8) of Sec. III. Note that we set j_ext = 0 throughout this section and that the ground state Ψ of Eq. (8) depends on 4N + 1 coordinates, i.e., Ψ = Ψ(r₁, σ₁, …, r_N, σ_N, p), where {σ_i} denotes the spin degrees of freedom and p is the elongation of the photon mode.

Within the dressed RDMFT, the original Hamiltonian (8) of N electrons in d dimensions and one mode is replaced by an extended auxiliary Hamiltonian of N dressed fermions in d + 1 dimensions with coordinates $z=(r,q)∈Rd+1$⁠. This auxiliary Hamiltonian reads

and gives access to the same physics (see Refs. 49, Sec. 4). For a d = 1 matter subsystem, the operators introduced in Eq. (20) read as follows: the dressed Laplacian $Δ′=∂2∂x2+∂2∂q2$⁠, the dressed local potential

and the dressed interaction kernel

The ground state Ψ′(z₁, σ₁, …, z_N, σ_N) = Ψ(x₁, σ₁, …, x_N, σ_n) ⊗ χ(p₂, …, p_N) of Ĥ′ is a product of the original physical ground state Ψ and the ground state of χ, which, in turn, is the product of N − 1 harmonic oscillator ground states. The auxiliary Hamiltonian (20) contains only one-body and two-body terms in terms of the dressed coordinates, which allows us, in principle, to use any standard electronic structure method to solve it (see also Refs. 49, Secs. 4 and 5). We use this construction to develop the dressed RDMFT and dressed Hartree–Fock (HF). For that, we define the dressed (spin-summed) first order reduced density matrix (1RDM)

To apply RDMFT theory on the auxiliary system, we have to replace the total energy functional of electronic RDMFT, given in Ref. 19, with the newly introduced quantities of the dressed system, i.e., the auxiliary Hamiltonian of Eq. (20) (approximately) evaluated by the dressed 1RDM γ′ of Eq. (23). By that, common approximations for the two-body energy expression in terms of the 1RDM can be directly transferred from electronic theory to the dressed system.⁵⁰ The minimization is performed like in the electronic case and is based on the RDMFT implementation of Octopus, though the convergence of the dressed system requires a more complicated protocol that can be found in the supplementary material of Ref. 49. The current implementation in Octopus approximates the conditions under which the dressed 1RDM corresponds to a wave function by ensuring the fermionic ensemble N-representability conditions.⁵¹ However, the auxiliary wave function also exhibits another exchange symmetry with respect to the auxiliary coordinates, which is currently neglected. For the practical validity of this approach, the reader is referred to Sec. 5 and the supplementary material of Ref. 49.

As an example, we consider the one-dimensional (1D) H₂ molecule in a soft-Coulomb potential with a slightly stretched bond-length of b = 2.0 bohrs⁵² that is modeled by the local potential

and the soft Coulomb interaction

In Fig. 4, we show the total energy and the total photon number of the dressed RDMFT, the dressed HF, and the exact many-body calculations⁵³ for different coupling strengths. We see that for small couplings, both observables are captured well by the dressed RDMFT and dressed HF. With an increase in the coupling strength, both approximations fail to capture the strongly increasing photon number. For the total energy instead, the dressed RDMFT remains very close to the exact result, whereas the deviations to the dressed HF increase with increasing coupling strength. This shows the potential of the dressed RDMFT to describe correlated electron systems that are strongly coupled to a cavity mode. In the future, we plan to investigate better approximations to the polaritonic N-representability conditions that also account for the symmetry of the many-body wave function with respect to the exchange of photon coordinates.

It is well known that the standard local and semilocal functionals of DFT tend to over-delocalize the electrons, as usual approximations are based on the homogeneous electron gas. This leads to several failures of DFT for materials in which the localization of electrons plays a critical role in dictating the system’s properties. This is, for instance, the case for transition metal oxides. The DFT+U method was originally proposed to compensate for some of the failures of the LDA for such materials.^54–57 In essence, the DFT+U method aims at a better description of the local electron–electron interaction, which is achieved by adding the mean-field Hubbard model on a chosen localized subspace to the DFT total energy. The double counting of electron interaction in this localized subspace is then removed. The DFT+U total energy functional reads

where E_(ee) is the usual electron–electron interaction energy and E_dc accounts for the double counting of the electron–electron interaction already present in E_DFT. The analytical expression of this double-counting term is not known in the general case, a general problem to all +U methods. Several approximated forms have been proposed along the years.^58,59 In the Octopus code, we implemented the most commonly used double-counting terms: the fully localized limit (FLL)⁶⁰ and the around-mean field (AMF) double-counting terms.⁶¹ They read as⁵⁹ 

and

respectively, where $Nσ=∑mnmmσ$ and N = N_↑ + N_↓. The E_ee and E_dc energies depend on the density matrix of a localized orbital basis set {ϕ_I,m}, which are the localized orbitals attached to the atom I. In the following, we refer to the elements of the density matrix of the localized basis as occupation matrices and we denote them as ${nmm′I,σ}$⁠. In the rotational-invariant form of DFT+U proposed by Dudarev et al.,⁶⁰ and for the FLL double-counting term, we obtain the E_U energy to be added to the DFT total energy, which only depends on an effective Hubbard U parameter U^eff = U − J,

where I is an atom index, σ is the spin index, and n, l, and m refer to the principal, azimuthal, and angular quantum numbers, respectively. In the case of a periodic system, the occupation matrices $nmm′I,n,l,σ$ are given by

where $w$_k is the k-point weight and $fnkσ$ is the occupation of the Bloch state $|ψn,kσ$⁠. Here, |ϕ_I,n,l,m⟩ are the localized orbitals that form the basis used to describe the electron localization. Details of the implementation can be found in Ref. 62. We recently extended our original implementation to be able to construct a localized subspace from localized states, such as Wannier states,⁶³ and to treat the intersite interaction.⁶⁴ 

In its usual formulation, the DFT+U method is an empirical method in which the effective U is a parameter of the calculation. However, it recently became possible to evaluate U and J fully ab initio and self-consistently using the ACBN0 functional.⁶⁵ We also implemented this method in Octopus and extended it to the time-dependent case⁶² to investigate strongly correlated materials out-of equilibrium. We showed that the absorption spectra of transition metal oxides, such as NiO or MnO, are well reproduced by our TDDFT+U simulations.⁶² 

Figure 5 shows the calculated time profile of the effective Hubbard U_eff = U − J for the 3d orbitals of Ni for light-driven NiO. The top panel exhibits the time profile of the driving vector potential. This shows that strongly driven correlated materials cannot be described by simply assuming that the effective electronic parameters (here, the effective Hubbard U) remain constant out of equilibrium. Moreover, the possibility to tune these effective electronic parameters offers opportunities for light-driven phase transitions, such as light-induced magnetic Weyl semimetals.⁶⁶ 

We finally note that TDDFT+U is not the only method that can be used to describe correlated materials such as transition metal oxides. An alternative approach to this problem would be to use hybrid functionals. However, such types of calculations in complex systems, potentially involving strong spin–orbit coupling interactions, are usually numerically very expensive in the context of real-time real-space TDDFT.

The van der Waals (vdW) interactions arise from correlations between electrons and are, in principle, described by the exact energy functional through the correlation energy functional E_c[n(r)]. However, the vdW interactions are inherently non-local and long-range and, by construction, cannot be described by usual local and semi-local functionals.⁶⁸ Therefore, much work has been devoted to finding consistent ways to enhance available functionals to correctly describe the necessary correlations yielding vdW forces.

Within a pure DFT approach, the inclusion of vdW interactions should be done through the correlation functional. To this end, a family of non-local density functionals has been proposed.⁶⁸ The so-called vdW density functionals (vdW-DF) are derived ab initio from the screened response of the homogeneous electron gas.⁶⁸ Another route to include vdW interactions in DFT calculations is by adding explicit corrections to the energy and forces based on atomic parameterizations. The most well-known method of this type is probably the vdW-D3 scheme of dispersion corrections from Grimme.⁶⁹ Although computationally more favorable than the vdW-DF method, the major disadvantage of this type of approach, particularly from the perspective of time-dependent simulations, is that it fails to correctly describe the effects that depend on the electron dynamics. A trivial example would be a lone electron traveling through space. Another scheme based on explicit corrections to the energy and forces was proposed by Tkatchenko and Scheffler (vdW-TS).⁷⁰ This scheme has the advantage of retaining much of the low computational cost of the vdW-D3 scheme while making the atomic parameterization dependent on the electronic density. Recently, these three schemes (vdW-DF, vdW-TS, and vdW-D3) were implemented in the Octopus code to deal with vdW interactions in isolated and periodic systems.

In the context of real-time TDDFT, it is important to note that these various levels of the description of vdW interactions have different implications. Indeed, in the vdW-D3 method, the energy and forces do not depend on the electronic density and will not change in time, unless atoms are moved, as already mentioned above. This should be a good approximation only when the time-dependent density remains close enough to the ground-state density. On the contrary, the expressions for the vdW-TS and the vdW-DF methods have a direct dependence on the electronic density and should thus yield a better description of vdW energies and forces out of equilibrium. However, we note that the use of the vdW-TS and vdW-DF methods within real-time TDDFT must be done using the adiabatic approximation and will, therefore, behave poorly when this approximation breaks down.

Octopus supports vdW-DF functionals⁶⁸ through the libvdwxc library.⁷¹ The vdW-DF functionals are expressed as a sum,

of the LDA correlation energy⁷² and a fully non-local correlation term. The latter is the integral over a kernel function ϕ(q₀, q₀′, r),

where q₀(r) depends on the local density and its gradient. An explicit evaluation of this six-dimensional integral is very expensive and scales as the volume squared, $O(N2)$⁠. Roman-Pérez and Soler proposed an efficient method, which approximates it as a sum of three-dimensional integrals.⁷³ The method works by expressing the integrand as a convolution using a limited set of helper functions and then applying the Fourier convolution theorem for a scaling of $O(N⁡log⁡N)$ and a much lower prefactor. This has since become the standard method for evaluating vdW-DF functionals.

As implemented in Octopus, the density is redistributed from its normal uniform grid of arbitrary shape onto a uniform 3D grid forming a cube or parallelepiped, which is suitable for 3D Fourier transforms. The library libvdwxc then evaluates the non-local energy and contributions to the potential, relying on the FFTW library⁷⁴ for efficient parallel Fourier transforms. After calculating the energy and potential, the potential is redistributed back to the original form.

Octopus supports the standard functionals vdW-DF1,^75,76 vdW-DF2,⁷⁷ and vdW-DF-cx,⁷⁸ as well as other common forms that differ by combining the correlation with a different exchange functional. Some common supported variations are vdW-DF-optPBE, vdW-DF-optB88,⁷⁹ and vdW-DF-C09.⁸⁰ 

Since the vdW-TS approach depends on the density, the effect of the van der Waals interaction can be observed in properties other than the forces. In particular, we expect to observe an effect in the excited states of systems that interact through vdW forces. We use the hydrogen fluoride dimer as a simple model system for a proof-of-concept application of the modular implementation of the vdW-TS functional correction on TDDFT calculations. The dimer geometry is setup as shown in Fig. 6. The hydrogen fluoride monomers are placed in an anti-parallel fashion, each one with its main symmetry axis oriented along the y-axis. The hydrogen fluoride bond in each monomer is 0.92 Å long, and the molecules are separated by 2.8 Å along the z-axis. At this distance, the van der Waals interaction between the monomers is the strongest according to the vdW-TS model. We choose this dimer model because it is a conveniently small system in which the effects of including the dispersion correction can be shown at an affordable computational cost.

To calculate absorption, we excite the system with an infinitesimal electric-field pulse and then propagate the time-dependent Kohn–Sham (TDKS) equations for 30.385 35 ℏ/eV. The singlet dipole spectrum is evaluated from the time-dependent dipole moment. The strength of the perturbation is set to 0.01 Å⁻¹, and the polarization is in the z-axis. The time evolution is carried out using the enforced time-reversal symmetry propagator with (default) time steps of 0.033 52 ℏ/eV.

The results (Fig. 6) show a small van der Waals-induced bathochromic-like (red) shift in the optical spectrum of the hydrogen fluoride dimer calculated with the LDA. This example opens the door for a new series of applications in supra-molecular chemistry, structural biology, and polymer science, which incorporate van der Waals effects on real-time electron dynamics.

Octopus also supports the DFT-D3 van der Waals correction.⁶⁹ This correction depends only on the atomic positions and does not depend on the atomic density. It is implemented in Octopus by linking to the DFT-D3 library provided by the authors. As we have validated Octopus results with the reference data provided with the library, the results are guaranteed to be consistent with other codes that implement this correction.

The only modification we did to the library was to move the very large set of coefficients from a source file to a stand alone text file that is parsed only when necessary. This speeds up compilation and reduces the size of the binaries. The modified library is distributed with Octopus, so users need not compile it separately.

The Polarizable Continuum Model (PCM)⁸¹ comprises a family of implicit-solvent approaches to tackle quantum-mechanical calculations of molecules in solution. The PCM assumes that (i) the solvent is a continuum and infinite dielectric medium characterized by a frequency-dependent dielectric function and (ii) a void cavity of appropriate shape and size encapsulates the solute molecule, separating it from the solvent by a sharp interface. The numerical implementation of the PCM relies on the Apparent Surface Charge (ASC) approach and the Boundary Element Method (BEM).⁸² In this framework, the solvent polarization response induced by the molecule’s charge density is modeled by a reaction potential defined by a set of point charges q = {q₁, q₂, …, q_T} that spread over a tessellated cavity surface consisting of T finite surface elements or tesserae.⁸³ The reaction potential is the object through which the complex dielectric environment is accounted for.⁸⁴ The implementation of the PCM in Octopus rests on the Integral Equation Formalism (IEF).⁸⁵ The key IEF-PCM equation for computing the polarization charges is

where q and V are column vectors of size T storing the induced polarization charges and the molecule’s electrostatic potential at the tesserae representative points, respectively. Q is the T × T PCM response matrix, which depends on the geometry of the cavity and the dielectric function of the solvent.⁸¹ 

A realistic description of the molecular cavity is key to capture accurate electronic and optical properties of molecules in solution. In principle, the cavity should (i) exclude the solvent, (ii) comprise most of the solute electronic density, and (iii) conform with the molecular shape. Octopus uses the GEPOL algorithm,⁸⁶ which builds up the van der Waals cavity from the union of interlocking spheres with element-specific radii centered at each atom position (by default, no spheres are built around hydrogen atoms). Within GEPOL, the BEM tessellation of the solute–solvent interface is done starting from a 60- or 240-face circumscribed polyhedron per sphere, selecting only exposed tesserae and properly reshaping those that are partially exposed. The BEM tessellation is a surface grid constructed independently from the real-space three-dimensional grid used in Octopus to represent both the Kohn–Sham electronic Hamiltonian and molecular orbitals. The mismatch between the two discrete representations might cause numerical problems arising from the Coulomb singularities whenever tesserae and grid points are close to each other. The implementation of the PCM in Octopus regularizes such singularities by using normalized spherical Gaussian functions to smooth the discretized polarization charges q.⁸⁴ 

Having briefly described the implementation, we now look at specific and relevant cases for chemical reactions, namely, solvation energies and ground-state stabilization. The most evident effect of the presence of a solvent is to change the total energy of a system compared to its value in vacuum. In the framework of DFT, the Kohn–Sham Hamiltonian of the solvated molecule contains the solvent reaction potential, which is a functional of the electronic and nuclear densities of the solute molecule.⁸⁴ The latter implies that the Kohn–Sham and IEF-PCM equations become coupled and the polarization charges q_e = QV_Hartree[ρ^e], induced by the molecule’s electronic density, have to be computed at each Kohn–Sham iteration until convergence is reached. The resulting electronic density can be used to compute the electrostatic contribution to the solvation free energy $ΔGel=G[ρsolve]−E[ρvace]$⁠, where G and E denote the free and total energy functionals of the molecule in solvent and in vacuum, respectively. The current PCM implementation in Octopus⁸⁴ has been thoroughly tested on a benchmark set of organic molecules and compared with analogous calculations performed with the quantum chemistry package GAMESS.⁸⁷ As an example of the solvent stabilization effect in the ground state, we studied the nitrobenzene molecule in water (static dielectric constant set to ϵ_s = 80).

In Fig. 7, we show the convergence of the solvation free energy as a function of the Kohn–Sham iteration. The solvation free energy is stabilized already after 10 Kohn–Sham iterations out of the 19 required to optimize the molecular orbitals of the solvated molecule. We have also verified that the numerical error inherent to the discretization of the solute cavity surface is very small. For example, the total polarization charge $QPCM=∑i=1Tqie$ at each Kohn–Sham iteration deviates by only 0.03% from the actual number of valence electrons in the nitrobenzene molecule (the relation between q and the solute charge is determined by Gauss’s theorem⁸¹). However, the magnitude of this error will depend, in general, on the size and the geometry of the molecule.

Now we move forward and describe the extension of the previous PCM method to the time domain using real-time PCM and non-equilibrium solvation. The ground state PCM is unable to capture the complex dynamical interactions between the solvent and the solute when the latter is in an excited state. Fortunately, the PCM admits a generalization to account for solvation dynamics.⁸⁹ The time-dependent PCM (TD-PCM) implementation in Octopus comes in three flavors of increasing complexity, all coupled with real-time TDDFT calculations of the solute molecules. The first one, called equilibrium TD-PCM,⁸⁴ assumes that the solvent is fast enough to instantaneously equilibrate the solute charge density fluctuations and to polarize accordingly. This approach is physically sound for weakly polar solvents, having similar values for the static and dynamic dielectric constants. The second is a nonequilibrium approach called inertial TD-PCM and amounts to partitioning the solvent response in a fast (dynamic) and a slow (inertial) part.⁹⁰ Faster degrees of freedom respond instantaneously to the changes of the applied potential (either of molecular or external in origin), whereas slower degrees of freedom remain “frozen” and in equilibrium with the initial value of the field. This approximation works well when the solvent relaxation times are large enough with respect to the electronic excitations in the solute molecule (e.g., within the picosecond scale). The third TD-PCM approach, called equation of motion (EOM) TD-PCM, considers the full history-dependent evolution of the solvent polarization through a set of equations of motion for the polarization charges.⁹¹ Nonequilibrium polarization effects of this sort originate from the frequency-dependent dielectric response of the solvent, encoding the fact that it takes a different non-negligible time to adjust to fast or slow electrostatic perturbations. Solvation dynamics affect strongly and non-trivially the absorption spectrum of molecules, especially for fast and polar solvents, by inducing solvatochromic shifts of the peaks in the UV–visible absorption spectrum and modifying their relative amplitudes. The details about the implementation of all of these schemes and a detailed discussion of their effects can be found in Refs. 84 and 92.

Here, we show the differences among the TD-PCM flavors described above for the case of the nitrobenzene molecule in aqueous solution, excited with an electrical dipole perturbation in the x direction. We take the dynamic dielectric constant of water, entering in the inertial and EOM TD-PCMs, as ϵ_d = 1.786. Figure 8 shows the photo-absorption spectrum of the gas-phase and solvatedmolecule. We can see that, with respect to the absorption in $v$acuo, all of the TD-PCM methods produce shifts of the features (in this case, toward lower excitation energies). The shifts are not rigid overall but depend on the excited state and on the specific solvation scheme (equilibrium, inertial, and EOM). Within Debye’s model, the equilibrium and inertial TD-PCMs are limiting cases for the dynamics when the relaxation time is zero and infinite, respectively. The EOM TD-PCM, with a finite relaxation time, interpolates between these limits. Water has a large relaxation time of 3.37 ps, and therefore, the absorption spectrum is almost coincident for the EOM and the inertial TD-PCMs, as shown in Fig. 8. Whenever the solvent is as slow as water, there is not much gain in selecting the EOM over inertial TD-PCM and the latter is the method of choice in terms of computational performance. Instead, for faster solvents, the EOM results depart from those of the inertial algorithm, approaching those of the equilibrium TD-PCM. In Fig. 8, we considered an effective solvent having the same static and dynamic dielectric constants as water, but with a 1 fs relaxation time. The EOM TD-PCM for such a model solvent produces almost the same excitation energy shift equilibrium TD-PCM and peak intensities in between the inertial and the equilibrium TD-PCMs, as expected.

A real-time representation of the molecular dipole coupled with the PCM surface charges is found in the supplementary material, SM2, which also highlights how the different TD-PCM approaches affect the time-evolution of the dipole.

When studying the evolution of a solvated molecule under the effect of an external time-dependent electromagnetic field, a separate treatment is required to take into account that there are two solvent polarization contributions interacting with the solute molecule, namely, the reaction and the cavity fields.⁹³ The reaction field comes from the polarization induced by the solute molecule itself, while the cavity field arises as a polarization induced by the external electric field. Both reaction and cavity field effects can be accounted for in static and real-time quantum-mechanical calculations of molecules within the PCM and TD-PCM simulations in Octopus.⁹² 

All TD-PCM calculations shown here were performed using both reaction- and cavity-field effects, although the redshift of the peaks in each TD-PCM scheme is mainly a feature of the reaction-field, as ascertained in test calculations (not shown) and as expected.⁹⁴ Cavity field effects impact directly on the peak intensities by making the absorption more favorable depending on how large the effective local field acting on the molecule is allowed to be by solvent dielectric properties and the geometry of the cavity (normally, reassembling the molecular shape). Still, cavity field effects can have a non-trivial influence in the absorption spectrum shape when considering non-equilibrium solvation dynamics (EOM TD-PCM) by changing the relative peak intensities. This effect can be seen in the inset of Fig. 8, where the normalized cavity field factor (CFF)—the ratio between absorption cross section with and without cavity field effects—is plotted against the excitation energy for the EOM TD-PCM simulations with water and the aforementioned faster water-like solvent. The large difference in absorption peak intensities between the equilibrium and the rest of the TD-PCMs in Fig. 8 is also related to the modification of the probing electromagnetic field inside the dielectric. The photo-absorption cross section is, by definition, the ratio between the absorbed and the incoming power of light. In our case, both increase with the dielectric constant: the larger the dielectric constant of a medium, the larger the CFF, but also a larger dielectric constant implies a smaller phase velocity of light, therefore increasing the power of the traveling electromagnetic wave. The impact of the latter effect (⁠$∝refractiveindex=(|ϵ(ω)|+Rϵ(ω))/2$⁠) is stronger than that of the former [∝CFF, e.g., for a spherical cavity⁹³ $∝3ϵ(ω)2ϵ(ω)+12$] for a large enough dielectric constant such as the one corresponding to the equilibrium TD-PCM for water (ϵ_s = 80).

In conclusion, Octopus is now capable of including an implicit dielectric continuum model as an environment for quantum mechanical calculations of excited states in the time-domain. The different versions of the TD-PCM scheme implemented allow us to select the most suitable to capture the relevant physics, accounting for a full range of different relaxation and response times.

In the last couple of years, the first studies investigating magnetization dynamics from first principles in real time have emerged.^95–100 Here, we are interested in transverse magnetic excitations, specifically magnons, which are long wavelength collective excitations with a typical energy of tens to hundreds of millielectron volts. We recently developed an alternative to the linear-response TDDFT formulation^98,99,101 to compute the spin susceptibilities of magnetic systems based on real-time TDDFT.¹⁰² Our approach follows the work of Bertsch et al. for optical excitations. In the original work of Bertsch et al., the system is perturbed by a sudden change in the vector potential, which induces a charge- and current-density response. To investigate magnons, we employ a “transverse magnetic kick,” which induces magnetic fluctuations in the system. To be more precise, our perturbation corresponds to an infinitely short application of a Zeeman term,

where q is the momentum of the spin wave we are exciting, $B$ is the strength of the perturbation, and σ^± = σ_x ± iσ_y are Pauli matrices. In this expression, the z axis is taken to be along the direction of the magnetization of the system before excitation. In case the system has a preferred magnetization direction due to the presence of spin–orbit coupling, usually referred as an “easy axis,” we perform a kick in the transverse direction with respect to this “easy axis” of the system.

The subsequent time evolution of the spin magnetization m(r, t), governed by the time-dependent Schrödinger equation, is then computed and analyzed in Fourier space. If we perturb our system from its ground state and we assume linear response, we have directly that

where χ₊₋(q; ω) is the spin susceptibility we want to extract. A similar expression is obtained for χ₋₊(q; ω).

In order to access finite momenta, one typically has to employ large supercells to perform the dynamics, which can be computationally very expensive when a few meV energy resolution is needed. Fortunately, there is a way to circumvent the construction of supercells by using the so-called generalized Bloch theorem (GBT). The GBT has been introduced by Sandradskii¹⁰³ for the calculation of ground-state spin waves and requires the implementation of specific boundary conditions. We implemented the GBT for the real-time calculation of magnons, taking advantage that Octopus is a real-space finite difference code, for which we can easily specify any type of boundary conditions. However, it is important to note that this applies only in the absence of spin–orbit coupling. The boundary condition, described below, depends on the momentum q of the perturbation and acts differently depending on whether a state was originally “up” or “down” with respect to the unperturbed magnetization. This is determined by the sign of ⟨Φ_nk|S_z|Φ_nk⟩ just before the perturbation for each spinor state |Φ_nk⟩. If we label these states α and β, the boundary condition reads

We checked that performing the simulation using the GBT or the supercell approach leads to the same results up to numerical precision. We tested our approach for cubic Ni, Fe, and Co, which have been widely studied using linear response TDDFT, and we found that our results are in very good agreement with previous work,⁹⁶ thus validating our implementation. Figure 9 shows the results obtained for bulk nickel using the adiabatic local-density approximation. We used here a real-space grid spacing of 0.27 bohr, norm-conserving pseudopotentials, and we employed a 16 × 16 × 16 k-point grid shifted 4 times in order to resolve momenta q that are multiples of 2π/(16a), where a is the lattice parameter of Ni, which we took to be 3.436 Å. In order to obtain the linear response, we took $B=0.02$ and we propagated for 435.4 fs using a time step of 1.81 as. To reduce the numerical burden, symmetries were employed.

Let us now comment on the interest of the proposed approach. So far, we only used this new method for investigating weak magnetic kicks, where we recover the results from linear response theory. However, our approach does not rely on the assumptions of small perturbations and can be used to investigate nonlinear phenomena induced by a strong magnetic field as well as out-of-equilibrium situations where the system is kicked from an excited state. This is the strength of a real-time method.^105–108 Moreover, we can investigate directly the coupling with other degrees of freedom, such as phonons or photons, without any new theory or code development. Using the real-space method, we benefit from the favorable scaling of the time propagation, which is linear in the number of states, whereas sum-over-state approaches usually scale quadratically with the number of states. Finally, our approach offers the great advantage of not requiring the use of an exchange–correlation kernel, as only the exchange–correlation potential is needed to perform a time propagation. This is very interesting in order to test new functionals and theory levels for which deriving the expression of the exchange–correlation kernel can become complicated. These different aspects will be investigated in future works.

In this section, we review the implemented routines for magneto-optical phenomena, which arise from the loss of symmetry between left and right circularly polarized light in the presence of a magnetic field.^109,110 While magneto-optical spectra can be straightforwardly computed for molecules,^111–115 the theory for solids has been developed only recently.¹¹⁶ The reason is that external electromagnetic fields break the translational symmetry of periodic systems, which is formally expressed through the unboundness of the position operator.^117–119 The orbital response to magnetic fields is especially complicated to describe, as such fields lead to non-perturbative changes in wavefunctions and introduce vector coupling to electron dynamics.^120–122 Here, we will focus on calculations of changes in the linear optical response in the presence of the magnetic field within the Sternheimer approach.^19,123–125

To treat uniform magnetic fields in periodic systems, we use the approach based on perturbation theory for the one-particle density matrix.^116,120,122 Such an approach allows us to work under purely periodic boundary conditions and to take into account automatically gauge invariance. A periodic and gauge-invariant counterpart $Õ$ is distinguished for any operator $O=Or1r2$ defined for two points r₁ and r₂ in real space

where A is the vector potential associated with external electric (E = −c⁻¹∂/∂tA, where c is the speed of light) and magnetic (B = ∇ × A) fields and the integral is taken along the straight line between points r₂ and r₁.

The time-dependent Liouville equation for the gauge-invariant counterpart $ρ̃$ of the one-particle density matrix to the first order in E and B takes the form¹¹⁶ 

Here, the commutator and anticommutator of operators $O(1)$ and $O(2)$ are denoted by $[O(1),O(2)]$ and ${O(1),O(2)}$⁠, respectively, the velocity operator $V=−i[r,H̃]$ is computed taking into account all non-local contributions to the Hamiltonian (e.g., non-local pseudopotentials), and the Hamiltonian is represented as $H̃=H0+δH̃$⁠, where the difference $δH̃$ between the gauge-invariant counterpart $H̃$ of the Hamiltonian and unperturbed Hamiltonian H₀ is related to the local-field effects. Unlike the singular position operator r, the commutator $[r,ρ̃]$ of the position operator with the periodic function $ρ̃$ is well defined in Eq. (38) and corresponds to the derivative with respect to the wave vector, $i∂kρ̃k$⁠, in reciprocal space. Differentiating Eq. (38), one finds the derivatives of the density matrix $ρ̃(P)=∂ρ̃/∂P$ with respect to external perturbations P (E, B, etc.).

For TDDFT calculations in Octopus, ρ is the Kohn–Sham density matrix. The nth order derivative $ρ̃(P)$ describing the joint response to the perturbations P = P₁P₂ … P_n is divided into four blocks within and between the occupied (V) and unoccupied subspaces (C): $ρ̃VV(P)=Pvρ̃(P)Pv$⁠, $ρ̃CC(P)=Pcρ̃(P)Pc$⁠, $ρ̃VC(P)=Pvρ̃(P)Pc$⁠, and $ρ̃CV(P)=Pcρ̃(P)Pv$⁠, where P_$v$ = ρ⁽⁰⁾ and P_c = 1 − P_$v$ are the projectors onto the occupied and unoccupied bands. In accordance with the density matrix perturbation theory,¹²⁶ to get the elements $ρ̃CV(P)$⁠, Eq. (38) is projected onto unperturbed Kohn–Sham wavefunctions $|ψvk(0)$ of occupied bands $v$ to give an equation for the response function $|ηvk(P)=ρ̃CV(P)(Ω)|ψvk(0)$⁠,

Here, the operator on the left-hand side is given by L_$v$k(Ω) = Ω + H₀ −ϵ_$v$k, where Ω is the frequency considered and ϵ_$v$k is the energy of the unperturbed state $|ψvk(0)$⁠. The operator R comes from the right-hand side of Eq. (38) and is determined by the derivatives of the density matrix of the previous orders. If the local-field effects are taken into account, the right-hand side R also depends on the derivative of the electron density n^(P)(r₁) = ρ^(P)(r₁, r₂)δ(r₁ − r₂). In this case, Eq. (39) needs to be solved self-consistently. The calculations, in practice, work with the periodic parts of the wavefunctions, $|uvk(0)$⁠. The commutator $[r,ρ̃]$ corresponding to $i∂kρ̃k$ in the reciprocal space is computed within the k · p theory.^19,124,125 Equation (39) is solved using the efficient Sternheimer approach,^19,123–125 where the function $|ηvk(P)(Ω)$ that fits into Eq. (39) is found iteratively at each frequency Ω. To avoid divergences at resonances, a small but finite imaginary frequency iδ is added to the frequency Ω₀ of the external perturbation so that Ω = Ω₀ + iδ.

Once the solution of Eq. (39) is known, the elements $ρ̃CV(P)$ are obtained as

The elements of $ρ̃(P)$ between the occupied and unoccupied subspaces are found using the relation $ρ̃VC(P)(Ω)=(ρ̃CV(P)(−Ω*))*$⁠, and for that, Eq. (39) is also solved for the frequency −Ω^*.

To find the elements within the occupied $ρ̃VV(P)$ and unoccupied $ρ̃CC(P)$ subspaces, the idempotency condition ρ = ρρ is used. In terms of the periodic counterpart $ρ̃$ of the density matrix and to the first order in the magnetic field, it is written as^120,122