N^{\prime}M^{\prime}\end{subarray}}^{sc}\cdot\mathbf{u}_{q\nu}e^{i[(\mathbf{k}+\mathbf{q})\cdot\mathbf{R}_{N}-\mathbf{k}\cdot\mathbf{R}_{N}^{\prime}]},$
(82)

where $C_{nM}$ are the LCAO coefficients and $\mathbf{u}_{q\nu}$ are mass-
scaled phonon displacement vectors.

###  VI.8 Raman spectrum

The Raman effect describes inelastic light scattering, where vibrational modes
are excited within the material. Resonant and non-resonant Raman spectra of
finite systems such as molecules can be calculated in various approximations
[232] using the corresponding interfaces in ASE. The Stokes Raman intensity is
then written as

$I(\omega)=I_{0}\sum_{\nu}\frac{n_{\nu}+1}{\omega_{\nu}}|\sum_{\alpha,\beta}u_{\mathrm{in}}^{\alpha}R_{\alpha\beta}^{\nu}u_{\mathrm{out}}^{\beta}|^{2}\delta(\omega-\omega_{\nu}),$
(83)

where $\nu$ denotes phonon mode at $\mathbf{q}=0$ with frequency of
$\omega_{\nu}$ and $n_{\nu}$ is the corresponding Bose-Einstein distribution.
Furthermore, $u_{\mathrm{in}}^{\alpha}$ and $u_{\mathrm{out}}^{\beta}$ are the
polarization vectors of the incoming and outgoing light, and
$R_{\alpha\beta}^{\nu}$ denotes the Raman tensor for phonon mode $\nu$.

The predominant approach for calculating $R_{\alpha\beta}^{\nu}$ involves the
use of the Kramers-Heisenberg-Dirac (KHD) method. Within the KHD framework,
the Raman tensor is determined by taking the derivative (utilizing a finite-
difference method) of the electric polarizability concerning the vibrational
normal modes. Alternatively, one can employ time-dependent third-order
perturbation theory to compute Raman tensors. These two approaches are
equivalent when local field effect are negligible [56]. However, each approach
comes with its own set of advantages and drawbacks. The KHD method is
computationally more efficient but is limited to computing first-order Raman
processes. The perturbative approach can be extended to higher-order Raman
processes involving multiple phonons, but it is more computationally
demanding, necessitating a greater number of bands and a finer k-mesh grid to
achieve convergence. The perturbative approach has been implemented in GPAW
and is elaborated below while the KHD method has been implemented in the ASR
package [22], utilizing GPAW as the computational backend.

In the perturbative approach, the Raman tensor $R_{\alpha\beta}^{\nu}$ is
given by [233, 56]

$\displaystyle
R_{\alpha\beta}^{\nu}\equiv\sum_{ijmn\mathbf{k}}\left[\frac{p_{ij}^{\alpha}(M_{jm}^{\nu}\delta_{in}-M_{ni}^{\nu}\delta_{jm})p_{mn}^{\beta}}{(\hbar\omega_{\mathrm{in}}-\varepsilon_{ji})(\hbar\omega_{\mathrm{out}}-\varepsilon_{mn})}\right.$
(84)
$\displaystyle+\frac{p_{ij}^{\alpha}(p_{jm}^{\beta}\delta_{in}-p_{ni}^{\beta}\delta_{jm})M_{mn}^{\nu}}{(\hbar\omega_{\mathrm{in}}-\varepsilon_{ji})(\hbar\omega_{\nu}-\varepsilon_{mn})}$
$\displaystyle+\frac{p_{ij}^{\beta}(M_{jm}^{\nu}\delta_{in}-M_{ni}^{\nu}\delta_{jm})p_{mn}^{\alpha}}{(-\hbar\omega_{\mathrm{out}}-\varepsilon_{ji})(-\hbar\omega_{\mathrm{in}}-\varepsilon_{mn})}$
$\displaystyle+\frac{p_{ij}^{\beta}(p_{jm}^{\alpha}\delta_{in}-p_{ni}^{\alpha}\delta_{jm})M_{mn}^{\nu}}{(-\hbar\omega_{\mathrm{out}}-\varepsilon_{ji})(\hbar\omega_{\nu}-\varepsilon_{mn})}$
$\displaystyle+\frac{M_{ij}^{\nu}(p_{jm}^{\alpha}\delta_{in}-p_{ni}^{\alpha}\delta_{jm})p_{mn}^{\beta}}{(-\hbar\omega_{\nu}-\varepsilon_{ji})(\hbar\omega_{\mathrm{out}}-\varepsilon_{mn})}$
$\displaystyle+\left.\frac{M_{ij}^{\nu}(p_{jm}^{\beta}\delta_{in}-p_{ni}^{\beta}\delta_{jm})p_{mn}^{\alpha}}{(-\hbar\omega_{\nu}-\varepsilon_{ji})(-\hbar\omega_{\mathrm{in}}-\varepsilon_{mn})}\right]$
$\displaystyle\times f_{i}(1-f_{j})f_{n}(1-f_{m}),$

where the first term is referred to as the resonant part and the remaining
terms represent different time orderings of the interaction in terms of
Feynman diagrams.
$p_{nm}^{\alpha}=\langle\psi_{n\mathbf{k}}|\hat{p}^{\alpha}|\psi_{m\mathbf{k}}\rangle$
is the momentum matrix element between electronic bands $m$ and $n$, with
transition energy $\varepsilon_{nm}=E_{n}-E_{m}$ in polarization direction
$\alpha$ and $M_{nm}^{\nu}$ is the electron–phonon coupling strength in the
optical limit $|\mathbf{q}|=0$ as defined in Eq. (82).

Fig. 11 shows the polarization-resolved Raman spectrum of bulk MoS2 in the 2H
phase as computed with a laser frequency of $\omega_{\mathrm{in}}=488$ nm.
This example uses only the resonant term, as the other contributions are small
in this case. It’s worth noting that we have conducted a comparison of the
calculated spectra using the ASR package and observed a high level of
agreement for several materials, for example MoS2. The results obtained from
both methods closely align with each other in terms of peak positions and
dominant peaks and minor disagreements between the two methods can be
attributed to differences in implementation details and the distinct
approximations employed by each. Specifically, within the ASR package, we
utilized the phonopy package [234] to compute phonon frequencies and
eigenvectors, whereas in GPAW, we directly computed phonon frequencies and
eigenvectors using ASE’s phonon module. Furthermore, in the ASR
implementation, we rigorously enforced the symmetry of the polarizability
tensor and the ASR results therefore typically exhibit a more accurate
adherence to the required symmetry of the Raman tensor compared to the GPAW
implementation.

Figure 11:  Polarization-resolved Raman spectrum of bulk MoS2 in the 2H phase
at $\omega=488$ nm. Phonons and potential changes were computed using a 700 eV
plane wave cutoff and $2\times 2\times 2$ $\mathbf{k}$-point mesh in a
$3\times 3\times 2$ supercell. Each peak has been labeled according to its
irreducible representation.

###  VI.9 Quadratic optical response functions

The nonlinear optical response of materials can be obtained by going beyond
first order perturbation theory. Presently, the GPAW implementation is
restricted to second-order response within the dipole approximation and
without inclusion of local field effects. We apply the independent particle
approximation, which cannot capture collective behavior such as excitonic
effects [57]. A spatially homogeneous incident electric field can be written
in terms of its Fourier components as

$\bm{\mathcal{E}}(t)=\sum_{\alpha,\omega_{1}}\mathcal{E}_{\alpha}(\omega_{1})\mathbf{e}_{\alpha}e^{-i\omega_{1}t}\,,$
(85)

where $\omega_{1}$ runs over positive and negative frequencies,
$\mathbf{e}_{\alpha}$ denotes the unit vector along the $\alpha$-direction,
and $\mathcal{E}_{\alpha}(\omega_{1})$ is the electric field at frequency
$\omega_{1}$. The induced quadratic polarization density
$\mathcal{P}_{\gamma}^{(2)}(t)$ can then be expressed as

$\displaystyle\mathcal{P}_{\gamma}^{(2)}(t)=$
$\displaystyle\epsilon_{0}\sum_{\omega_{1}\omega_{2}}\sum_{\alpha\beta}\chi_{\gamma\alpha\beta}^{(2)}(\omega_{1},\omega_{2})$
$\displaystyle\qquad\qquad\times\mathcal{E}_{\alpha}(\omega_{1})\mathcal{E}_{\beta}(\omega_{2})e^{-i(\omega_{1}+\omega_{2})t}\,,$
(86)

where $\chi_{\gamma\alpha\beta}^{(2)}$ is the rank-3 quadratic susceptibility
tensor. Due to intrinsic permutation symmetry, i.e.
$\chi_{\gamma\alpha\beta}^{(2)}(\omega_{1},\omega_{2})=\chi_{\gamma\beta\alpha}^{(2)}(\omega_{2},\omega_{1})$,
$\chi_{\gamma\alpha\beta}^{(2)}$ has at most 18 independent elements which may
be further reduced by the Neumann principle and point group symmetries [235].
We note that the corresponding quadratic conductivity tensor is readily
derived from $\chi_{\gamma\alpha\beta}^{(2)}$ as
$\sigma_{\gamma\alpha\beta}^{(2)}(\omega_{1},\omega_{2})=-i\epsilon_{0}\chi_{\gamma\alpha\beta}^{(2)}(\omega_{1},\omega_{2})(\omega_{1}+\omega_{2})$
due to relationship between current density and polarization density [236].

Among the various response functions that can be calculated from
$\chi_{\gamma\alpha\beta}^{(2)}(\omega_{1},\omega_{2})$, we have implemented
second-harmonics generation (SHG) and the shift current tensor. The
implementation currently requires time-reversal symmetry, which limits the
application to non-magnetic systems. For SHG, the susceptibility tensor is
separated into a pure interband term $\chi_{\gamma\alpha\beta}^{(2e)}$ and a
mixed term $\chi_{\gamma\alpha\beta}^{(2i)}$ that read

$\displaystyle\chi_{\gamma\alpha\beta}^{(2e)}\equiv
C_{0}\sideset{}{{}^{\prime}}{\sum}_{\mathbf{k},nml}\dfrac{\mathrm{Re}\\{r_{nm}^{\gamma}(r_{ml}^{\alpha}r_{ln}^{\beta}+r_{ml}^{\beta}r_{ln}^{\alpha})\\}}{2(\varepsilon_{ln}-\varepsilon_{ml})}$
$\displaystyle\qquad\times\bigg{(}\frac{2f_{nm}}{2\hbar\omega-\varepsilon_{mn}}-\frac{f_{nl}}{\hbar\omega-\varepsilon_{ln}}+\frac{f_{ml}}{\hbar\omega-\varepsilon_{ml}}\bigg{)}\,,$
(87) $\displaystyle\chi_{\gamma\alpha\beta}^{(2i)}\equiv
iC_{0}\sideset{}{{}^{\prime}}{\sum}_{\mathbf{k},nm}f_{nm}\Bigg{[}\frac{2\mathrm{Im}\\{r_{nm}^{\gamma}(r_{mn;\beta}^{\alpha}+r_{mn;\alpha}^{\beta})\\}}{\varepsilon_{mn}(2\hbar\omega-\varepsilon_{mn})}$
$\displaystyle\qquad+\frac{\mathrm{Im}\\{r_{mn}^{\alpha}r_{nm;\beta}^{\gamma}+r_{mn}^{\beta}r_{nm;\alpha}^{\gamma}\\}}{\varepsilon_{mn}(\hbar\omega-\varepsilon_{mn})}$
$\displaystyle\qquad+\frac{\mathrm{Im}\\{r_{nm}^{\gamma}(r_{mn}^{\alpha}\Delta_{mn}^{\beta}+r_{mn}^{\beta}\Delta_{mn}^{\alpha})\\}}{\varepsilon_{mn}^{2}}$
$\displaystyle\qquad\qquad\qquad\times\bigg{(}\frac{1}{\hbar\omega-\varepsilon_{mn}}-\frac{4}{2\hbar\omega-\varepsilon_{mn}}\bigg{)}$
$\displaystyle\qquad-\frac{\mathrm{Im}\\{r_{mn}^{\alpha}r_{nm;\gamma}^{\beta}+r_{mn}^{\beta}r_{nm;\gamma}^{\alpha}\\}}{2\varepsilon_{mn}(\hbar\omega-\varepsilon_{mn})}\Bigg{]}\,.$
(88)

Here, $C_{0}\equiv e^{3}/(2\epsilon_{0}V)$ where $V$ is the crystal volume,
$\Delta_{mn}^{\alpha}\equiv(p_{mm}^{\alpha}-p_{nn}^{\alpha})/m$ denotes the
velocity difference between bands $n$ and $m$, and $r_{nm}^{\alpha}$ is the
interband ($n\neq m$) position matrix element, obtained from
$imr_{nm}^{\alpha}=\hbar p_{nm}^{\alpha}/\varepsilon_{nm}$. All energies,
occupations and matrix elements in the preceding expressions depend on the
$\mathbf{k}$-vector. Also, the summation over $\mathbf{k}$ implies an integral
over the first BZ, i.e. $(2\pi)^{3}\sum_{\mathbf{k}}\rightarrow
V\int_{\mathrm{BZ}}\mathrm{d}^{3}{\mathbf{k}}$. The primed summation signs
indicate omission of terms with two or more identical indices. Finally, the
generalized derivative $r_{nm;\alpha}^{\beta}$ (for $n\neq m$) are evaluated
from the sum rule [237]

$\displaystyle r_{nm;\alpha}^{\beta}$
$\displaystyle=\hbar\frac{r_{nm}^{\alpha}\Delta_{mn}^{\beta}+r_{nm}^{\beta}\Delta_{mn}^{\alpha}}{\varepsilon_{nm}}$
$\displaystyle+\frac{i\hbar}{\varepsilon_{nm}}\sum_{l\neq
n,m}(\varepsilon_{lm}r_{nl}^{\alpha}r_{lm}^{\beta}-\varepsilon_{nl}r_{nl}^{\beta}r_{lm}^{\alpha})\,.$
(89)

Here, infinite sums have been substituted with finite sums over a limited, yet
sizable, set of bands. It’s important to emphasize that both sides of the
equation are dependent on the $\mathbf{k}$-vector, and the summation on the
right-hand side pertains exclusively to bands. It should be mentioned that
another implementation for the quadratic susceptibility tensor in the velocity
gauge is also available in the code, but is not documented here. For
sufficiently many bands, the results of the two implementations are identical
[236].

Regarding the shift current, where a DC current is induced in response to an
incident AC field, one needs to compute the quadratic conductivity tensor
$\sigma_{\gamma\alpha\beta}^{(2)}$:

$\displaystyle\sigma_{\gamma\alpha\beta}^{(2)}$ $\displaystyle=\frac{\pi
e^{3}}{\hbar^{2}V}\sideset{}{{}^{\prime}}{\sum}_{\mathbf{k},nm}f_{nm}\mathrm{Im}\\{r_{mn}^{\alpha}r_{nm;\gamma}^{\beta}+r_{mn}^{\beta}r_{nm;\gamma}^{\beta}\\}$
$\displaystyle\qquad\qquad\times\delta(\hbar\omega-\varepsilon_{mn})\,.$ (90)

In practice, the delta function is replaced by a Lorentzian with a finite
broadeining $\eta$. To avoid numerical instabilities, the implementation of
Eqs. (VI.9)-(VI.9) uses a tolerance, such that terms are neglected if the
associated energy differences or differences in Fermi levels are smaller than
a tolerance. The default value of the tolerance is $10^{-6}$ eV, and $10^{-4}$
for the energy and Fermi level differences, respectively.

##  VII Real-time TDDFT

Figure 12: Photoabsorption spectrum for a Ag147 icosahedral nanoparticle and
the transition contribution map at 3.8 eV. The map reveals that transitions
between KS states near the Fermi level contribute constructively to the
plasmon resonance, while transitions from occupied states at the $d$-band edge
contribute destructively (screening).

The real-time TDDFT (RT-TDDFT) scheme, also known as time-propagation TDDFT,
is implemented in FD [34] and LCAO [52] modes. It requires non-periodic
boundary conditions, but is not restricted to the linear regime and can be
applied to model molecules in strong fields. The metod may combined with
hybrid quantum-classical modeling to simulate the dynamical interaction
between molecules and plasmon resonances at metal surfaces [238]. The LCAO-RT-
TDDFT is the more recent implementation supporting versatile analyses and
enabling modeling of large systems thanks to its efficiency [52, 53, 239, 54,
240, 241, 242, 243, 244, 245]. We focus on the capabilities of the LCAO
version in this section, but some of the described functionalities are also
available in FD mode.

The time-dependent KS equation in the PAW formalism is

$\displaystyle
i\hat{S}\frac{\mathrm{d}}{\mathrm{d}t}|\tilde{\psi}_{n}(t)\rangle=\left(\hat{H}[n(\mathbf{r})]+\hat{v}(t)\right)|\tilde{\psi}_{n}(t)\rangle,$
(91)

where the Kohn–Sham Hamiltonian $\hat{H}[n(\mathbf{r})]$ is implicitly
dependent on time through the time-dependent density and $\hat{v}(t)$ is an
explicitly time-dependent external potential. We have additionally assumed
that the overlap matrix $\hat{S}$ is independent of time, i.e. there are no
ion dynamics.

Starting from the ground state, Eq. (91) is propagated forward numerically.
After each step, a new density is computed and $\hat{H}[n(\mathbf{r})]$
updated accordingly. The user can can freely define the external potential and
implemented standard potentials include the delta-kick
$\hat{v}(t)=\hat{x}\delta(t)$ and a Gaussian pulse
$\hat{v}(t)=\hat{x}\exp[-\sigma^{2}(t-t_{0})^{2}/2]\cos[\omega(t-t_{0})]$,
where $\hat{x}$ is the dipole operator in the direction $x$.

During the propagation, different time-dependent variables can be recorded,
and after the propagation they can be post-processed to quantities of physical
and chemical interest. As a basic example, the time-dependent dipole moment
recorded for a delta-kick perturbation can be Fourier-transformed to yield the
photoabsorption spectrum [246]. Observables are recorded by attaching
observers to the calculation and implemented observers include writers for the
dipole moment, magnetic moment, KS density matrix in frequency space and wave
functions.

RT-TDDFT calculations can be started from the ground state or continued from
the last state of a previous time propagation and the time-limiter feature
allows one to limit jobs to a predefined amount of wall time. Together with
the continuation capability, this facilitates time propagation in short
chunks, efficiently using shared high-performance resources.

In the LCAO-RT-TDDFT implementation, Eq. (91) is cast into a matrix equation
and solved with ScaLAPACK [52]. The intermediate step of updating the Hartree
and XC potential is performed on the real-space grid.

###  VII.1 Kohn–Sham decomposition

The time-dependent KS density matrix can be written as

$\displaystyle\begin{split}\rho_{nn^{\prime}}(t)=\sum_{m}f_{m}\times&\int\left(\psi^{(0)}_{n}(\mathbf{r})\right)^{*}\psi_{m}(\mathbf{r},t)\mathrm{d}{\mathbf{r}}\\\
\times&\int\psi^{*}_{m}(\mathbf{r}^{\prime},t)\psi^{(0)}_{n^{\prime}}(\mathbf{r^{\prime}})\mathrm{d}{\mathbf{r}^{\prime}},\end{split}$
(92)

where $\psi^{*}_{m}(\mathbf{r}^{\prime},t)$ are the time-dependent KS orbitals
with ground state occupation factors $f_{m}$. The KS density matrix is a
central quantity enabling computation of observables and may be evaluated
efficiently in the LCAO mode [54].

The Fourier transform of the induced density matrix
$\delta\rho_{nn^{\prime}}(\omega)=\mathcal{F}[\rho_{nn^{\prime}}(t)-\rho_{nn^{\prime}}(0)]$
can be built on the fly during time propagation through the density-matrix
observer. Details on the implementation are described in Ref. [54]. The KS
density matrix in frequency space is related to the Casida eigenvectors and
gives similar information as the solution of the Casida equation [54].
Observables such as the polarizability can be decomposed into a sum over the
electron–hole part of the KS density matrix $\rho_{nn}$, where
$f_{n}>f_{n^{\prime}}$. This enables illustrative analyses, e.g., by
visualizing $\rho_{nn^{\prime}}(\omega)$ on energy axes as a transition
contribution map [247], from which the nature of the localized surface plasmon
resonance can be understood (Fig. 12, see Ref. [54] for detailed discussion).

###  VII.2 Hot-carrier analysis

The KS density matrix is a practical starting point for analyzing hot-carrier
generation in plasmonic nanostructures. In the regime of weak perturbations,
the KS density matrix resulting from arbitrary pulses can be obtained from
delta-kick calculations by a post-processing routine [243, 244]. By
decomposing the matrix into different spatial and energetic contributions,
hot-carrier generation in nanostructures [243, 245] and across nanoparticle-
semiconductor [241] and nanoparticle-molecule [240, 244] interfaces can be
studied. Computational codes and workflows for such hot-carrier analyses are
provided in Refs. [248, 249].

###  VII.3 Circular dichroism for molecules

Electronic circular dichroism (ECD) is a powerful spectroscopic method for
investigating chiral properties at the molecular level. The quantity that
characterizes the ECD is the rotatory strength, which is defined through the
magnetic dipole moment $m(\omega)$ as

$R(\omega)=\frac{1}{\pi\kappa}\mathrm{Re}\bigg{[}\sum_{\alpha}m_{\alpha}^{(\alpha)}(\omega)\bigg{]},$
(93)

where the index $\alpha$ enumerates Cartesian directions, the $\alpha$
superscript in parenthesis indicates the $\delta$-kick direction and $\kappa$
is the strength of the kick. To resolve $R(\omega)$, one needs to perform the
$\delta$-kick in all three Cartesian directions using a perturbing electric
field $\mathbf{E}^{\alpha}(t)=\kappa\hat{\mathbf{\alpha}}\delta(t)$. The
frequency components of the dipole moment $m_{i}^{(i)}(\omega)$ are calculated
by Fourier transforming $m_{\alpha}^{(\alpha)}(t)$, which is recorded during
the propagation. Finally, the time-dependent magnetic dipole moment is
obtained as the expectation value of the operator
$\hat{\mathbf{m}}=-\frac{i}{2c}\hat{\mathbf{r}}\times\hat{\nabla}$,

$\mathbf{m}(t)=\sum_{n}f_{n}\int\psi_{n}^{*}({\mathbf{r}},t)\,\hat{\mathbf{m}}\,\psi_{n}({\mathbf{r}},t)~{}{\rm
d}\mathbf{r},$ (94)

where $f_{n}$ is the occupation number of KS orbital $n$ and $\psi_{n}({\bf
r},t)$ is the time-evolved KS state. The current GPAW implementation supports
both the FD and LCAO mode and the computational efficiency of the LCAO mode
enables calculation of the ECD of nanoscale metal-organic nanoclusters. More
details on the implementation can be found in Ref. [55].

###  VII.4 Radiation reaction potential

Plasmonic or collective molecular excitations are strongly susceptible to any
kind of optical interaction. Induced currents will couple via Maxwell’s
equation of motion to the optical environment and result in radiative losses,
i.e., decay towards the ground state. It is possible to solve the Maxwell
problem formally by obtaining Green’s tensor $\textbf{G}_{\perp}(\omega)$. The
dipolar interaction between electric field and electronic dipole can be
absorbed into a local potential
$v_{\text{rr}}(\textbf{r},t)=e\textbf{r}\cdot\textbf{E}_{r,\perp}(t)=e\textbf{r}\cdot\big{[}\mathcal{F}_{t}^{-1}(i\mu_{0}\omega\textbf{G}_{\perp}(\omega))\ast\partial_{t}\textbf{R}(t)\big{]}$
suitable for Kohn–Sham TDDFT, where $\textbf{R}(t)$ is the total electronic
dipole moment. A detailed discussion can be found in Ref. [250].

For many simple structures, such as free-space, one-dimensional wave-guides,
or dielectric spheres, $\textbf{G}_{\perp}(\omega)$ is analytically known and
radiative losses can then be included in TDDFT without additional
coimputational cost. The tutorials on the GPAW web-page [61] includes an
example for one-dimensional wave-guides for which the user can specify the
cross sectional area and the polarization of the propagating modes. Extending
the functionality of the radiation-reaction potential to 3D free-space and
collective interaction of large ensembles in Fabry-Pérot cavities from first
principles [251] is essential for the understanding of polaritonic chemistry.
This functionality is currently under development.

###  VII.5 Ehrenfest dynamics

Molecular dynamics (MD) simulations usually rely on the Born-Oppenheimer
approximation, where the electronic system is assumed to react so much faster
than the ionic system that it reaches its ground state at each time step.
Thus, forces for the dynamics are calculated from the DFT ground-state
density. While this approximation is sufficiently valid in most situations,
there are cases where the explicit dynamics of the electronic system can
affect the molecular dynamics, or the movement of the atoms can affect
averaged spectral properties. These cases can be handled using so-called
Ehrenfest dynamics, ie. time- dependent density functional theory molecular
dynamics (TDDFT/MD).

Ehrenfest dynamics is implemented in the FD mode [252]. A description of the
theory and a tutorial is available on the GPAW web-page [61]. This
functionality has been used to model the electronic stopping of ions including
core-electron excitations [253], study charge transfer at hybrid interfaces in
the presence of water [254], simulate the coherent diffraction of neutral atom
beams from graphene [255], model the dependence of carbon bond breaking under
Ar+-ion irradiation on $sp$ hybridization [256], and reveal charge-transfer
dynamics at electrified sulfur cathodes [257]. An LCAO implementation,
inspired by recent work in the Siesta code [258, 259], is currently under
development.

##  VIII Excited-state DFT methods

###  VIII.1 Improved virtual orbitals

The linear response TDDFT approach generally provides reasonably accurate
excitation energies for low-lying valence excited states, where the orbitals
associated with the holes and excited electrons overlap significantly.
However, it tends to fail for excitations involving spatial rearrangment of
the electrons, such as charge transfer [260, 113], Rydberg [261] and doubly
excited states [262].

Some of these problems can be alleviated by using range separated functionals
(see section III.3.4). However, these functionals come with a significantly
increased computational cost due to the evaluation of exchange integrals.
Moreover, due to the missing cancellation of Coulomb and exchange terms for
canonical unoccupied orbitals within Hartree-Fock theory, one obtains spurious
unoccupied orbitals. This leads to slow covergence of linear-response TDDFT
calculations with respect to the number of unoccuppied orbitals when hybrid
and range-separated functionals are used [122, 123].

Substantial improvement in convergence with respect to unoccupied orbitals can
be obtained using improved virtual orbitals as devised by Huzinaga and Arnau
[263, 264]. In this approach, a modified Fock operator is used for the
unoccupied orbitals, which mimics the interaction between a hole arbitrarily
chosen among the occupied ground state orbitals and the excited electron. This
leads to faster convergence in linear-response TDDFT calculations with hybrid
and range-separated functionals, and also makes it possible to evaluate
excited state properties. For example, the energetics of long-range charge
transfer can be obtained by means of a ground state calculations because the
difference between the energy of an improved virtual orbital and a hole tends
to approximate the excitation energy. The improved virtual orbitals approach
is available in GPAW and details on the implementation are described in Ref.
[123].

###  VIII.2 Variational excited-state calculations

GPAW also offers the possibility to perform excited-state calculations using
an alternative time-independent density functional approach [265] (sometimes
referred to as the “$\Delta$SCF” method), which does not suffer from the
limitations of linear-response TDDFT mentioned in the previous section. The
method involves variational optimization of the orbitals corresponding to a
specific excited state by optimizing the electronic energy to a stationary
point other than the ground state. The computational effort is similar to that
of ground-state calculations and the variational optimization guarantees that
the Hellmann–Feynman theorem is fulfilled. Therefore, all the ground-state
machinery available in GPAW to evaluate atomic forces can be used for geometry
optimization and for simulating the dynamics of atoms on the excited state.
Furthermore, coupling this time-independent, variational approach for excited-
state calculations with external MM potentials (see section IV) does not
involve additional implementation efforts compared to ground-state
calculations and provides a means for performing excited-state QM/MM molecular
dynamics simulations that include the state-specific response of a solvent,
i.e. the response due to changes in the electron density of the solute.

Variationally optimized excited states correspond to single Slater
determinants of optimal orbitals with a non-aufbau occupation and are
typically saddle points on the electronic energy surface. Hence, variational
calculations of excited states are prone to collapsing to lower-energy
solutions, which preserve the symmetry of the initial guess. A simple maximum
overlap method (MOM) [266, 267] is available in GPAW to address this problem.
At each SCF step, the MOM occupies those orbitals that overlap the most with
the orbitals of a non-aufbau initial guess, usually obtained from a ground-
state calculation. The MOM, however, does not guarantee that variational
collapse is avoided, and convergence issues are common when using SCF
eigensolvers with density-mixing schemes developed for ground-state
calculations.

####  VIII.2.1 Direct orbital optimization

To alleviate the issues leading to variational collapse and achieve more
robust convergence to excited-state solutions, GPAW contains two alternative
strategies that are more reliable than conventional SCF eigensolvers with the
MOM. They are based on direct optimization of the orbitals and use saddle
point search algorithms akin to those for transition-state searches on
potential energy surfaces of atomic rearrangements. These approaches also
facilitate variational excited-state calculations of nonunitary invariant
functionals, such as self-interaction corrected functionals (see section
III.3.5).

The first of these methods is a direct orbital optimization approach
supplemented with the MOM (DO-MOM). This method is an extension of the direct
minimization approach using the exponential transformation illustrated in
section III.2.5 where the search is for a generic stationary point of
$E[\mathbf{\Psi}]$ instead of a minimum:

$\underset{\mathbf{\Psi}}{\rm
stat\,}E[\mathbf{\Psi}]=\underset{\mathbf{\Psi}_{0}}{\rm
min\,}\underset{\mathbf{A}}{\rm
stat\,}\mathcal{F}[\mathbf{A},{\mathbf{\Psi}_{0}}]\,.$ (95)

The DO-MOM is available in GPAW for both LCAO [51, 268], real-space grid, and
plane wave basis sets [50]. For the LCAO basis, the excited-state optimization
only necessitates making the energy stationary with respect to the elements of
the anti-Hermitian matrix $\mathbf{A}$ (the orbital rotation angles), while
calculations using the real-space grid and plane-wave basis include an outer-
loop minimization with respect to the reference orbitals
${\mathbf{\Psi}_{0}}$. The optimization in the linear space of anti-Hermitian
matrices uses efficient quasi-Newton algorithms that can handle negative
Hessian eigenvalues and therefore converge on saddle points. GPAW implements a
novel limited-memory SR1 (L-SR1) algorithm, which has proven to be robust for
calculations of excitations in molecules [50, 51].

The DO-MOM relies on estimating the degrees of freedom along which the energy
needs to be maximized starting from an initial guess. For valence and Rydberg
excitations, the initial guess consisting of ground-state canonical orbitals
with non-aufbau occupation numbers and preconditioning with a diagonal
approximation of the electronic Hessian using the orbital energies [79] can be
sufficient. However, if the excitation involves significant charge transfer,
large rearrangements of the energy ordering of the orbitals can occur, and DO-
MOM can struggle to converge. A second direct optimization method with
generalized mode following (DO-GMF) [49] alleviates these problems and is also
implemented. In DO-GMF, the components of the energy gradient $\mathbf{g}$
along the modes $\mathbf{v}_{i}$ corresponding to the $n$ lowest eigenvalues
of the electronic Hessian are inverted yielding

$\mathbf{g}^{\mathrm{\,mod}}=\mathbf{g}-2\sum_{i=1}^{n}\mathbf{v}_{i}\mathbf{v}_{i}^{\mathrm{T}}\mathbf{g}\,,$
(96)

and a minimization using the modified gradient $\mathbf{g}^{\mathrm{\,mod}}$
is performed by following the $n$ modes simultaneously. This procedure
guarantees convergence to an $n^{\rm th}$-order saddle point, eliminating the
risk of variational collapse altogether. While it is computationally more
expensive than DO-MOM due to the need for a partial Hessian diagonalization,
DO-GMF is more robust. Hence, it is particularly useful in the exploration of
potential energy surfaces because it is able to follow an excited state
through bond-breaking configurations, where broken-symmetry solutions appear,
by targeting the solution that preserves the saddle point order. This
important advantage is exemplified by the challenging double-bond twisting in
ethylene [49, 133], where DO-GMF calculations of the lowest doubly excited
state provide an avoided crossing with the ground state whereas other methods
fail to do so.

####  VIII.2.2 Example applications of direct optimization

Figure 13: Applications of time-independent, variational calculations of
excited states to a molecule in vacuum (left), a molecule in solution (middle)
and a solid-state system (right). Left: Deviation of the calculated excitation
energy of a charge-transfer excited state in the N-phenylpyrrole molecule from
the theoretical best estimate of 5.58 eV [269]. The results of linear-response
TDDFT calculations with hybrid functionals are from Ref. [269]. Middle: Time-
evolution of the interligand angles of the [Cu(dmphen)2]+ complex upon
photoexcitation to the lowest metal-to-ligand charge-transfer (MLCT) state in
acetonitrile. The experimental results from femtosecond X-ray scattering
measurements [270] are compared to the average over excited-state molecular
dynamics trajectories obtained using the QM/MM electrostatic embedding scheme
in GPAW [165] (see section IV) and convoluted with the experimental
instrument-response function [270]. Right: Vertical excitation energy for
excitations in the negative nitrogen-vacancy center in diamond obtained with
the r2SCAN functional as compared to the results of previous calculations
using an advanced quantum embedding approach [271]. The orbitals involved in
the electronic transitions are visualized in the inset (C atoms are grey and
the N atom is orange).

The efficiency and robustness of the direct-optimization approaches combined
with the possibility of choosing different basis-set types make variational
calculations of excitated states in GPAW applicable to a great variety of
systems ranging from molecules in gas phase or solution to solids.

State-specific orbital relaxation enables the description of challenging
excitations characterized by large density rearrangements. Figure 13 shows the
error on the vertical excitation energy of a charge-transfer excitation in the
twisted N-phenylpyrrole molecule [49] obtained with direct optimization in
GPAW using the LDA, PBE and BLYP functionals and an sz+aug-cc-pVDZ [53] basis
set as compared to the results of linear-response TDDFT calculations with the
same basis set and functionals, as well as the hybrid functionals PBE0 and
B3LYP (results from Ref. [269]). For the variational calculations, the energy
of the singlet excited state is computed using the spin-purification formula
$E_{\mathrm{s}}=2E_{\mathrm{m}}-E_{\mathrm{t}}$, where $E_{\mathrm{m}}$ is the
energy of a spin-mixed state obtained by promoting an electron in one spin
channel and $E_{\mathrm{t}}$ is the energy of the triplet state with the same
character. The variational calculations underestimate the theoretical best-
estimate value of the excitation energy (5.58 eV) in Ref. [269] by 0.15-0.3
eV, an error that is significantly smaller than that of linear-response TDDFT
calculations with the same functionals ($-2.0$ eV) or with the more
computationally intensive PBE0 hybrid functional ($-0.85$ eV) [269].

The method has also been used to simulate the photoinduced structural changes
of photocatalytic metal complexes and concomitant solvation dynamics [270,
165, 182, 272, 181]. Figure 13 shows an application to the prototypical copper
complex photosensitizer [Cu(dmphen)2]+
(dmphen=2,9-dimethyl-1,10-phenanthroline) in acetonitrile, where QM/MM
molecular dynamics simulations elucidated an intricate interplay between
deformation of the ligands and rearrangement of the surrounding solvent
molecules following a photoexcitation to a metal-to-ligand charge-transfer
state [270, 165].

The last example of an application shown in Figure 13 is a calculation of the
excited states of a solid state system [273], the negatively charged nitrogen-
vacancy center in diamond, which comprises a prototypical defect for quantum
applications. The system is described with a large supercell of up to 511
atoms, and the calculations use a plane-wave basis set. In contrast to
previous reports, a range of different density functionals is found to give
the correct energy ordering of the excited states, with the r2SCAN functional
providing the best agreement with high-level many-body quantum embedding
calculations with an error of less than $0.06$ eV [273, 271]. This example
shows that the direct optimization methods in GPAW are promising tools for
simulating excited states in extended systems, where alternative approaches
are either computationally expensive or lack accuracy.

##  IX Other features

###  IX.1 Electric polarization

The formal polarization of bulk materials may be calculated from the modern
theory of polarization [274, 275] as

$\mathbf{P}_{\mathrm{F}}=\mathbf{P}_{\mathrm{F}}^{\mathrm{el}}+\mathbf{P}_{\mathrm{F}}^{\mathrm{n}},$
(97)

where

$\mathbf{P}_{\mathrm{F}}^{\mathrm{el}}=\frac{e}{(2\pi)^{3}}\mathrm{Im}\int_{\mathrm{BZ}}\mathrm{d}\mathbf{k}\sum_{n}f_{\mathbf{k}n}\langle
u_{\mathbf{k}n}|\nabla_{\mathbf{k}}|u_{\mathbf{k}n}\rangle$ (98)

is the electronic contribution and

$\mathbf{P}_{\mathrm{F}}^{\mathrm{n}}=\frac{1}{V_{\mathrm{cell}}}\sum_{a}Z^{a}\mathbf{r}^{a}$
(99)

is the contribution from the nuclei. Here the sums run over atoms in the unit
cell and $Z^{a}$ is the charge of nucleus $a$ (including core electrons),
situated at position $\mathbf{r}^{a}$. The electronic contribution can be
viewed as a Brillouin-zone integral of $\mathbf{k}$-space Berry phases and may
be evaluated from a finite-difference version of Eq. (98) [276]. This involves
the overlaps between Bloch orbitals at neighbouring $\mathbf{k}$-points, which
are straightforward to evaluate in the PAW formalism [46].

Eq. (97) is only defined modulo $e\mathbf{R_{i}}/V_{\mathrm{cell}}$, which
follows from the arbitrary choice of unit cell for the atomic positions as
well as the choice of phases for $u_{\mathbf{k}n}$, which can shift the Berry
phase by $2\pi$. The change in polarization under any adiabatic deformation
is, however, well defined and may be calculated as

$\Delta\mathbf{P}=\int_{0}^{\lambda}\mathrm{d}\lambda\frac{\mathrm{d}\mathbf{P}_{\mathrm{F}}}{\mathrm{d}\lambda},$
(100)

where $\lambda$ is some dimensionless variable parameterizing the adiabatic
path. In particular, for ferroelectrics the spontaneous polarization
$\mathbf{P}_{\mathrm{S}}$ can be evaluated by choosing a path that deforms the
structure from the polar ground state at $\lambda=1$ to a non-polar structure
at $\lambda=0$.

In Fig. 14, we show an example of this for tetragonal KNbO3, which is a well-
known ferroelectric [277]. The polar structure was relaxed under the
constraint of tetragonal symmetry using PBE ($\lambda=1$) and then linearly
interpolated to the inverted structure ($\lambda=-1$) passing through a
centrosymmetric point at $\lambda=0$. There are infinitely many polarization
branches differing by the polarization quantum $ec/V_{\mathrm{cell}}$ ($c$
being the lattice constant in the $z$-direction) and the spontaneous
polarization is obtained by choosing a single branch and evaluating the
difference in formal polarization at $\lambda=1$ and $\lambda=0$.
Interestingly, the centrosymmetric point has a non-vanishing polarization
given by half the polarization quantum. This is allowed due to the multi-
valued nature of the formal polarization and such a “topological polarization”
in non-polar materials has been shown to yield gapless surface states [278,
279]. Here, however, we merely use the topological polarization to emphasize
the importance of evaluating the spontaneous polarization as the change in
$\mathbf{P}_{\mathrm{F}}$ along an adiabatic path.

Figure 14: Formal polarization along an adiabatic path connecting two states
of polarization and the energy along the path in tetragonal KNbO3. The
spontaneous polarization is obtained as the difference in polarization between
a polar ground state ($\lambda=1$) and a non-polar reference structure
($\lambda=0$). The energy along the path is also shown.

The expressions (97-100) can be applied to extract various properties of non-
polar materials as well, for example the Born effective charge tensors

$Z^{a}_{\alpha\beta}=V_{\mathrm{cell}}\frac{\partial P_{\beta}}{\partial
r_{\alpha}^{a}},$ (101)

which yield the change in polarization resulting from small shifts in atomic
positions. In GPAW, these are obtained by a simple call to a module that
introduces a (user-defined) shift of all atoms in the unit cell and calculates
the resulting change in polarization from Eq. (97-99). The Born charges may be
combined with the atomic force matrix to calculate equilibrium positions under
an applied static electric field and the lattice contribution to the dynamic
polarizability can be calculated from the eigenvalues and eigenvectors of the
force matrix [280]. In addition, the piezoelectric response can be obtained by
calculating the change in polarization in response to external strain [45].
The lattice contribution to the polarizability is typically orders of
magnitude smaller than the electronic part, but for ferroelectrics the soft
phonon modes associated with spontaneous polarization can give rise to
significant lattice polarizabilities.

We exemplify this by the well known case of 2D ferroelectric GeS [281, 282,
283] where we obtain a spontaneous polarization of 490 pC/m [284, 44] in
excellent agreement with previous calculations [283]. The lattice and
electronic 2D in-plane polarizabilities are
$\alpha^{\mathrm{2D}}_{\mathrm{lat}}=4.32$ Å and
$\alpha^{\mathrm{2D}}_{\mathrm{el}}=3.75$ Å respectively [45]. For comparison,
the non-polar case of 2D MoS2 yields
$\alpha^{\mathrm{2D}}_{\mathrm{lat}}=0.09$ Å and
$\alpha^{\mathrm{2D}}_{\mathrm{el}}=6.19$ Å [45].

###  IX.2 Berry phases and band topology

Topological phases such as the quantum spin Hall state and Chern insulator
depend crucially on presence of spin–orbit coupling and the band topology may
be obtained from the evolution of $\mathbf{k}$-space Berry phases across the
Brillouin zone. In particular, for insulators the eigenvalues of the Berry
phase matrix of occupied states

$\gamma_{mn}^{i}(\mathbf{k}_{\perp})=i\int_{0}^{2\pi}\mathrm{d}k_{i}\langle
u_{m}(\mathbf{k})|\partial_{k_{i}}|u_{n}(\mathbf{k})\rangle$ (102)

must change by an integer multiple of $2\pi$ when a component of
$\mathbf{k}_{\perp}$ (components of $\mathbf{k}$ orthogonal to $k_{i}$) is
cycled through the Brillouin zone [285]. Here, $u_{n}(\mathbf{k})$ are spinor
Bloch states, which are typically not smooth functions of $\mathbf{k}$. The
evaluation of Eq. (102) thus requires the construction of a smooth gauge and
in GPAW this is handled by the parallel-transport algorithm of Ref. [286].

The method has been applied to high throughput search for new topological two-
dimensional materials [46] and in Fig. 15 we show the calculated Berry phases
of 1T’-MoS2 [287]. Due to time-reversal symmetry, the Berry phases at $\Gamma$
and $Y$ are two-fold degenerate, and therefore no perturbation that conserves
time-reversal symmetry can open a gap in the Berry-phase spectrum. This
property is diagnostic for the quantum spin Hall insulating state and closely
related to the presence of gapless edge states [285].

Figure 15: Berry phases of the quantum spin Hall insulator 1T’-MoS2 obtained
from PBE with non-selfconsistent spin–orbit coupling. The colors indicate the
expectation value of $S_{z}$ for each state as defined in Ref. [46] .

The eigenvalues of Eq. (102) may also be used to calculate the electronic
contribution to the formal polarization, since the sum of all individual Berry
phases yields the same value as one may obtain from Eq. (98). The present
approach is, however, more involved since it requires the construction of a
smooth gauge, which is not needed in Eq. (98).

###  IX.3 Wannier functions

Wannier functions (WFs) provide a localized representation of the electronic
states of a solid. The WFs are defined by a unitary transformation of the
Bloch eigenstates that minimises the spatial extent of the resultant orbitals.
Specifically, the $n$th Wannier function in unit cell $i$ is written as

$w_{in}(\mathbf{r})=\frac{1}{N_{\mathbf{k}}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}_{i}}\tilde{\psi}_{n\mathbf{k}}(\mathbf{r}),$
(103)

where $\tilde{\psi}_{n\mathbf{k}}$ is a generalized Bloch function (a
superposition of Bloch eigenstates at $\mathbf{k}$).

The minimization of the spatial extent of a set of WFs
$\\{w_{n}(\mathbf{r})\\}_{n=1}^{N_{w}}$ is equivalent to the maximization of
the spread functional [288]

$\Omega=\sum_{n=1}^{N_{w}}\sum_{\alpha=1}^{N_{q}}W_{\alpha}|Z_{\alpha,nn}|^{2},$
(104)

where

$Z_{\alpha,nn}=\langle
w_{n}|e^{-i\mathbf{q}_{\alpha}\cdot\mathbf{r}}|w_{n}\rangle.$ (105)

The $\\{\mathbf{q}_{\alpha}\\}$ is a set of at most 6 reciprocal vectors
connecting a $\mathbf{k}$-point to its neighbors and $W_{\alpha}$ are
corresponding weights accounting for the shape of the unit cell [289].

The generalized Bloch functions of Eq. (103) are determined by minimizing
$\Omega$ using e.g. a conjugate gradient scheme as implemented in the ASE
Wannier module. The inputs to this Wannierization algorithm are the matrices

$\displaystyle Z_{\alpha,ij}^{(0),\mathbf{k}}$ $\displaystyle=$
$\displaystyle\langle\tilde{\psi}_{\mathbf{k}i}|e^{-i\mathbf{q}_{\alpha}\cdot\mathbf{r}}|\tilde{\psi}_{\mathbf{k}+\mathbf{q}_{\alpha},j}\rangle$
(106) $\displaystyle+$
$\displaystyle\sum_{aii^{\prime}}e^{-i\mathbf{q}_{\alpha}\cdot\mathbf{R}^{a}}\Delta
S_{ii^{\prime}}^{a}$
$\displaystyle\times\langle\tilde{\psi}_{\mathbf{k}i}|\tilde{p}_{i}^{a}\rangle\langle\tilde{p}_{i^{\prime}}^{a}|\tilde{\psi}_{\mathbf{k}+\mathbf{q}_{\alpha},j}\rangle,$

where $\Delta S_{ii^{\prime}}^{a}$ are the PAW corrections from Eq. 8. From
these matrices, the ASE Wannier module can be used to construct partially
occupied Wannier functions [290, 291], which are a generalization of
maximally-localized Wannier functions [286] to entangled bands and non-
periodic systems. Recently, a further improvement in terms of robustness of
the Wannierisation procedure was achieved using a modified spread functional
containing a penalty term proportional to the variance of the spread
distribution of the WFs, which leads to a more uniform spread distribution
[292].

###  IX.4 Point defect calculations with hybrid functionals

Point defects play a crucial role in many applications of semiconductors[293,
294]. First-principles calculations can be used to determine the atomic
structure, formation energy, and charge-transition levels of points defects.
It is well established that the best description of point defects in
semiconductors/insulators is obtained using range separated hybrids, such as
the HSE06 xc-functional[5, 295]. To illustrate the use of GPAW for point
defect calculations, we determine the formation energy diagrams of the CN and
CB defects in the hexagonal boron nitride (hBN) crystal with the HSE06
functional. These defects have been proposed to be responsible for the deep-
level luminescence signal with a zero-phonon line (ZPL) around 4.1 eV.[296,
295]. The results are compared to similar results obtained with the VASP
software package.

For a point-defect $\mathrm{D}$ in charge state $q$, the formation energy Ef
is calculated from the formula,

$E^{\mathrm{f}}\left[\mathrm{D}^{q}\right]=E_{\mathrm{tot}}\left[\mathrm{D}^{q}\right]-E_{\mathrm{tot}}^{\mathrm{bulk}}-\sum_{i}n_{i}\mu_{i}+qE_{\mathrm{F}}+E_{\mathrm{corr}}.$
(107)

Here $E_{\mathrm{tot}}\left[\mathrm{D}^{q}\right]$ and
$E_{\mathrm{tot}}^{\mathrm{bulk}}$ are the total energies of the crystal with
the point defect in charge state $q$ and of the neutral pristine crystal,
respectively. $\mu_{i}$ is the chemical potential of the element $i$ while
$n_{i}$ is the number of atoms added to ($n_{i}>0$) or removed from
($n_{i}<0$) the crystal to create the defect. $E_{\mathrm{F}}$ is the chemical
potential of the electrons, i.e. the Fermi level, which is written as
$E_{\mathrm{F}}=E_{\mathrm{VBM}}+\Delta E_{\mathrm{F}}$, where VBM is the
valence band maximum. Finally, Ecorr is a correction term which accounts for:
(i) the spurious electrostatic interaction between the periodic images of the
defect and their interaction with the compensating homogeneous background
charge and (ii) the potential shift between the pristine and defect system.

For more details on the methodology of point defect calculations we refer the
reader to the excellent review papers on the topic [297, 298, 299, 300].

Figure 16: Defect formation energies for CN and CB, under N-rich and N-poor
conditions, respectively. The dashed lines are reproduced from Ref. 295.

All calculations have been performed using the HSE06 functional with the
default mixing parameter $\alpha$=0.25, plane wave cut-off of 800 eV, and
forces converged to 0.01 eV/Å. The lattice of the hBN crystal was fixed at the
experimental parameters ($a=2.50$ Å and $c=6.64$ Å) [301]. The band gap of the
pristine crystal was determined to be 5.58 eV (using $8\times 8\times 4$
$k$-points) in good agreement with the experimental band gap of 6.08 eV [302].
The structure of the point defects were relaxed in a $4\times 4\times 2$ (128
atom) supercell using $\Gamma$-point $k$-point sampling. For each defect three
different charge states ($q=1,0,-1$) were considered. The corrections
($E_{\mathrm{corr}}$) due to image charges and potential alignment are
evaluated following Freysoldt, Neugebauer and Van de Walle [303] as
implemented in GPAW.

Figure 16 shows the defect formation energies as a function of Fermi level for
CN and CB, at N-rich and N-poor conditions, respectively. We can see that
under N-rich conditions, CB is energetically lower, whereas CN is favorable
under N-poor conditions. CB shows a 1/0 charge transition at 3.73 eV above the
VBM, whereas CN has a 0/-1 charge transition 3.26 eV deep inside the band gap.
We find good agreement with a similar study [295] also employing plane waves
and the HSE06 functional (VASP calculations). Minor discrepancies can be
attributed to use of a different supercell size and slightly higher fraction
of the non-local mixing parameter ($\alpha$=0.31).

###  IX.5 Point-group symmetry representations

GPAW allows for the automatized assignment of point-group symmetry
representations for the pre-computed Kohn–Sham wavefunctions. This can be used
for determining the wave-function symmetry representations for both molecules
[304] and extended structures [305] to analyze, for example, the symmetry-
involved degeneracy of the bands and selection rules for dipole transitions.

The analysis follows directly from group theory [306], stating that the
solutions to the Schrödinger equation inherit the symmetry group of the
respective Hamiltonian, or essentially the external potential invoked by the
atomic configuration. The representation matrices $\Gamma$ are computed as

$\Gamma_{n}(T)=\int\phi_{n}^{\dagger}(\mathbf{r})P(T)\phi_{n}(\mathbf{r})\mathrm{d}\mathbf{r},$
(108)

where $\phi$ is a normalized wavefunction, $n$ is the eigenstate label, and
$P(T)$ is an operation that corresponds to the transformation $T$ of the
symmetry group of the Hamiltonian. The operations $P(T)$ include rotations
that are non-trivial for the rectangular grid that the computed wavefunctions
are projected onto. The wavefunction rotations are performed on the grid by
cubic interpolation. The output of the analysis contains the irreducible
representation weights $c_{\alpha,n}$ for each eigenstate $n$ as solved from

$\Gamma_{n}(T)=\sum_{\alpha}c_{\alpha,n}\chi_{\alpha}(T),$ (109)

where $\chi_{\alpha}$ are the character vectors of the group (i.e. rows of the
character table).

When doing the analysis, the user needs to input the coordinates of the center
of symmetry (typically the coordinates of a single atom), and the point group
for which the analysis is run. It is possible to analyze only a part of the
wave function by selecting a cutoff radius from the center of symmetry beyond
which the parts of the wave function are neglected. This enables the
investigation of the purity of the local symmetry even if the symmetry of the
Hamiltonian is broken far from the center of symmetry [304, 305]. To date,
point groups of C2, C2v, C3v, D2d, D3h, D5, D5h, I, Ih, Oh, Td, and Th are
implemented.

###  IX.6 Band-structure unfolding

When studying defect formation, charge-ordered phases or structural phase
transitions, it is often needed to perform DFT calculations on a super-cell. A
super-cell (SC) calculation comes with the cost of having to account for many
more electrons in the unit cell when compared to the primitive cell (PC). This
implies that besides the increased computational effort, the band structure of
a SC contains more bands in a smaller Brillouin zone as compared to the PC. In
order to compare electronic band structures between SC and PC, unfolding the
band structure of the SC into the one of the primitive cell (PC) becomes
convenient.

GPAW features the possibility of performing band-structure unfolding in the
real-space grid, plane wave, and LCAO modes. The implementation allows to
unfold the SC band structure without the explicit calculation of the overlap
between SC and PC wavefunctions, following the procedure described in Ref.
[307]. The unfolded band structure is given in terms of the spectral function

$A(\mathbf{k},\epsilon)=\sum_{m}P_{\mathbf{K}m}(\mathbf{k})\delta(\epsilon_{\mathbf{K}m}-\epsilon),$
(110)

with $\mathbf{k}$, $\mathbf{K}$ momenta in the PC and SC Brillouin zone
respectively, $\epsilon_{\mathbf{K}m}$ the SC eigenvalues obtained for
momentum $\mathbf{K}$ and band index $m$, and $P_{\mathbf{K}m}(\mathbf{k})$ is
calculated as

$P_{\mathbf{K}m}(\mathbf{k})=\sum_{\\{\mathbf{G}\\}}|C_{\mathbf{K}m}(\mathbf{G}+\mathbf{k}-\mathbf{K})|^{2},$
(111)

where $C_{\mathbf{K}m}$ are the Fourier coefficients of the eigenstate
$|\mathbf{K}m\rangle$ and $\\{\mathbf{G}\\}$ the subset of reciprocal space
vectors of the SC that match the reciprocal space vectors of the PC. A more
detailed explanation and technical details on how to perform a band-structure
unfolding can be found on the GPAW web-page [61].

###  IX.7 The QEH model

The quantum-electrostatic heterostructure (QEH) [308] model is an add-on GPAW
feature for calculating the dielectric response and excitations in vertical
stacks of 2D materials, also known as van der Waals (vdW) heterostructures.
The QEH model can be used independently from the GPAW code, but it relies on
the GPAW implementation for the calculation of the fundamental building blocks
used by the model, as elaborated below.

The dielectric screening in 2D materials is particularly sensitive to changes
in the environment and depends on the stacking order and thickness of the 2D
heterostructure, providing a means to tune the electronic excitations,
including quasi-particle band gaps and excitons. While the dielectric response
of freestanding layers can be explicitly represented _ab initio_ in GPAW at
the linear-response TDDFT, GW, and BSE level of theory, lattice mismatch
between different 2D layers often results in large supercells that make these
many-body approaches infeasible. Since the interaction between stacked layers
is generally governed by van der Waals interactions, the main modification to
the non-interacting layers’ dielectric response arises from the long-range
electrostatic coupling between layers.

Therefore, in the QEH model, the dielectric function of the vdW
heterostructure is obtained through an electrostatic coupling of the quantum
dielectric building blocks of individual 2D layers [309]. The dielectric
building blocks consist of monopole and dipole components of the density-
response function of the freestanding layers,
$\tilde{\chi_{i}}(\mathbf{q}_{\parallel},\omega)$, calculated from _ab initio_
on the RPA level. Subsequently, the full density-response function
$\chi_{i,j}(\mathbf{q}_{\parallel},\omega)$ (density perturbation on layer $i$
due to a monopole or dipole component of the perturbing field acting on layer
$j$) is calculated by solving the Dyson equation

$\displaystyle\chi_{i,j}(\mathbf{q}_{\parallel},\omega)$
$\displaystyle=\tilde{\chi_{i}}(\mathbf{q}_{\parallel},\omega)\delta_{i,j}+$
$\displaystyle\tilde{\chi_{i}}(\mathbf{q}_{\parallel},\omega)\sum_{k\neq
i}V_{i,k}(\mathbf{q}_{\parallel})\chi_{k,j}(\mathbf{q}_{\parallel},\omega),$
(112)

where the Coulomb matrix element is obtained as the real-space overlap over
the density, $\rho$, and potential, $\phi$, basis functions on the different
layers.

$V_{i,k}(\mathbf{q}_{\parallel})=\int\rho_{i}(z,\mathbf{q}_{\parallel})\phi_{k}(z,\mathbf{q}_{\parallel})\mathrm{d}z.$
(113)

From the density-response function, the dielectric function is obtained in the
basis of monopole/dipole perturbations on each layer in the heterostructure.
While building blocks pre-computed with GPAW for a large variety of 2D
materials are provided with the QEH package, GPAW offers the possibility of
calculating custom building blocks for any 2D material as explained in Ref.
[310].

As an illustrative example of the QEH model, we show how the static dielectric
function of a heterostructure can be engineered by multi-layer stacking. Fig.
17 shows the static dielectric function of a vdW heterostructure made up by
stacking $N$ MoS2 layers and $N$ WSe2 layers. We see that the dielectric
function increases significantly as a function of the number of layers,
eventually approaching a bulk limit. The knowledge of the layer-dependence of
the dielectric response for such a heterostructeres could be further exploited
to investigate inter- and intra-layer excitonic properties and band-edge
renormalization effects.

Figure 17: The static macroscopic dielectric function of a vdW heterostructure
interface as a function of the number of layers. The heterostructure is made
up of $N$ MoS2 layers on one half and $N$ WSe2 layers on the other half (see
inset). Increasing the number of layers eventually leads to a bulk-like limit
for the chosen stacking configuration.

###  IX.8 Solvent models

The presence of a solvent has a large effect on the energetics and the
electronic structure of molecules or extended surfaces. In particular, the
arguably most important solvent, water, is able to stabilize ions or
zwitterions that would not form in the gas phase. The main effect relates to
the large permittivity of water ($\varepsilon_{r}=78$) that effectively
screens Coulomb interactions.

A convenient and computationally lean method to describe this effect is the
inclusion of a position-dependent solvent permittivity
$\varepsilon(\mathbf{r})$ in the electrostatics via the Poisson solver [311].
The solvent is represented solely as a polarizable continuum that averages out
all movements and re-arrangements of the solvent molecules and their
electrons. The computational cost is, therefore, practically the same as a
calculation in vacuum.

This implementation allows the calculation of solvation free energies of
neutral and ionic species in solution [311]. Further, it can be applied to
periodic surfaces interfaced with an electrolyte to reproduce reasonable
potential drops within the simulation of electrochemical reaction processes,
as we will elaborate in the following.

###  IX.9 Charged electrochemical interfaces

Simulating atomistic processes at a solid–liquid interface held at a
controlled electrode potential is most appropriately performed in the
electronically grand-canonical ensemble [312, 313, 314, 315]. Here, electrons
can be exchanged dynamically with an external electron reservoir at a well-
defined electrochemical potential. In a periodic system, a non-zero net charge
would lead to divergence of the energy; therefore, any fractional electrons
that are added to or removed from the system must be compensated by an equal
amount of counter charge. Several approaches able to account for this change
in boundary conditions have recently been brought forward [316, 317, 318, 319,
320]. In GPAW, this is conveniently accomplished with the introduction of a
jellium slab of equal and opposite charge to the needed electronic charge; the
jellium is embedded in an implicit solvent localized in the vacuum region
above the simulated surface (cf. IX.8). As a particular highlight of this
Solvated Jellium Method (SJM) [321], as it is known in GPAW, we are able to
localize the excess charge on only the top side of the simulated atomistic
surface, which occurs naturally by introducing the jellium region solely in
the top-side vacuum and electrostatically decoupling the two sides of the cell
via a dipole correction. Both a purely implicit and hybrid explicit-implicit
solvent can be applied.

In the SJM, the simulated electrode potential is a monotonic function of the
number of electrons in the simulation; calculations can be run in either a
constant-charge or a constant-potential ensemble. The electrode potential
($\phi_{\mathrm{e}}$) within SJM is defined as the Fermi level ($\mu$)
referenced to an electrostatic potential deep in the solvent (the solution
inner potential $\Phi_{\mathrm{w}}$), where the whole charge on the electrode
has been screened and no electric field is present,

$\phi_{\mathrm{e}}=\Phi_{\mathrm{w}}-\mu.$ (114)

We can relate $\phi_{\mathrm{e}}$ to the commonly used reference potentials,
for example to the standard hydrogen electrode, by subtracting its absolute
potential such as the experimental value of 4.44 V as reported by Trasatti
[322]. In practice, the reference potentials depend on the used solvent model
[323], and the reference can be calibrated using computed and measured
potentials of zero charge [324].

The energy used in the analysis of electrode reactions is the grand-potential
energy $\Omega$,

$\Omega\equiv E_{\mathrm{tot}}+\phi_{\mathrm{e}}N_{\mathrm{e}},$ (115)

where $N_{\mathrm{e}}$ are the excess electrons in the simulation. This allows
for energetic comparisons between calculations with different numbers of
electrons, and is the default energy returned to atomistic methods by SJM.
While $E_{\mathrm{tot}}$ is consistent with the forces in traditional
electronic-structure calculations, the grand-potential energy $\Omega$ is
consistent with the forces in constant-potential simulations [325, 315]. This
means that relaxations that follow forces will correctly find local minima in
$\Omega$, and any kind of structure optimization or dynamical routine can be
performed on the grand-canonical potential energy surface, such as the search
for saddle points [326, 327, 328] or molecular-dynamics simulations [329].

In constant-potential mode, the potential is controlled by a damped iterative
technique that varies $N_{\mathrm{e}}$ to find the target $\phi_{\mathrm{e}}$;
in practice, a trajectory (such as a relaxation or nudged elastic band) is
run, where in a first series of SCF cycles the potential equilibrates. Upon
achieving a target potential within the given threshold, the code will conduct
the chosen geometry-optimization routine under constant potential control.

###  IX.10 Constrained DFT

Constrained DFT (cDFT) [330, 331, 332] is a computationally efficient method
for constructing diabatic or charge/spin-localized states. GPAW includes a
real-space implementation of cDFT [333], which can be used in both the FD and
LCAO modes. Compared to most cDFT implementations, in GPAW the periodicity can
be chosen flexibly between isolated molecules and periodic systems in one-,
two-, or three-dimensions with $\mathbf{k}$-point sampling.

The key difference of cDFT compared to normal DFT is the introduction of an
auxiliary potential to force a certain region (in real space around a
molecule, molecular fragment, or atom) to carry a predefined charge or spin.
This leads to a modified energy functional

$\begin{split}F[n(\mathbf{r}),\\{V_{i}\\}]&=E^{\mathrm{KS}}[n(\mathbf{r})]+\\\
&\sum_{i}V_{i}\sum_{\sigma}\left(\int\mathrm{d}\mathbf{r}w^{\sigma}_{i}(\mathbf{r})n^{\sigma}(\mathbf{r})-N_{i}\right),\end{split}$
(116)

where $E^{\mathrm{KS}}$ is the Kohn–Sham energy functional, $\sigma$ denotes
the spin, $n^{\sigma}(\mathbf{r})$ is the spin-dependent electron density, and
$N_{i}$ is the predefined charge or spin constraint. $V_{i}$ acts as Lagrange
multiplier, which determines the strength of the auxiliary potential and needs
to be determined self-consistently as discussed below.
$w^{\sigma}_{i}(\mathbf{r})$ is the weight function that defines how the
charge or spin are to be partitioned, i.e. the regions where charge/spin is to
be localised. This will be discussed in some detail below.

Introducing the constraining term in Eq. (116) leads to a new localized, spin-
dependent external potential

$v_{\text{eff}}^{\sigma}(\mathbf{r})=\dfrac{\delta F[n(\mathbf{r})]}{\delta
n(\mathbf{r})}=\dfrac{\delta E^{\mathrm{KS}}[n(\mathbf{r})]}{\delta
n(\mathbf{r})}+\sum_{i}V_{i}\sum_{\sigma}w^{\sigma}_{i}(\mathbf{r}).$ (117)

The constraint is further enforced by demanding that the ${V_{i}}$ need to
satisfy the chosen constraints

$C\geq\bigg{\lvert}\sum_{\sigma}\int\mathrm{d}\mathbf{r}w^{\sigma}_{i}(\mathbf{r})n^{\sigma}(\mathbf{r})-N_{i}\bigg{\rvert}.$
(118)

In practice, the strength of the constraining potential ${V_{i}}$ is found
through a self-consistent two-stage optimization of both $\\{V_{i}\\}$ and
$n(\mathbf{r})$. As the derivatives of $F[n(\mathbf{r}),\\{V_{i}\\}]$ with
respect to ${V_{i}}$ are readily available [333], gradient-based optimization
algorithms in SciPy [68] are used for optimizing $\\{V_{i}\\}$. The weight
functions are defined by a Hirshfeld-type partitioning scheme with Gaussian
atomic densities, and $w^{\sigma}_{i}(\mathbf{r})$ and the resulting external
potential are presented on the grid. With these definitions, the forces
resulting from the cDFT external potential can be analytically computed and
used in e.g. geometry optimization or molecular-dynamics simulations.

cDFT has been widely used for computing transfer rates within Marcus theory,
which depends on the reorganization and reaction (free) energies and the
diabatic coupling matrix element; the GPAW-cDFT implementation includes all
the needed tools for obtaining these parameters for bulk, surface, and
molecular systems [333, 334]. Recently, the cDFT approach has been combined
with molecular-dynamics methods to compute the reorganization energy at
electrochemical interfaces [335] as well as with grand-canonical DFT methods
(see section IX.9) to construct fixed electron-potential diabatic states
[336].

###  IX.11 Orbital-free DFT

Orbital-free DFT (OFDFT) approximates the DFT energy functional by modelling
the kinetic energy as a direct functional of the density

$\begin{split}E_{\text{OF}}[n]=\underbrace{\int\text{d}\textbf{r}\,n^{1/2}(\textbf{r})\left(-\frac{1}{2}\nabla^{2}\right)\,n^{1/2}(\textbf{r})}_{T_{\text{W}}[n]}+J[n]+V[n]\\\
+E_{\text{xc}}[n]+T_{\text{s}}[n]-T_{\text{W}}[n].\end{split}$ (119)

Levy and colleagues showed that a Kohn–Sham-like equation derived
variationally from the equation above holds for the square-root of the density
[337]

$\left(-\frac{1}{2}\nabla^{2}+V_{\text{eff}}(\textbf{r})\right)\,n^{1/2}(\textbf{r})=\mu\,n^{1/2}(\textbf{r}).$
(120)

OFDFT approximately enforces the Pauli principle, partially accounting for
quantum effects in an averaged way.

The OFDFT scheme implemented in GPAW offers the advantage of accessing all-
electron values while maintaining computational linear-scaling time with
respect to system size. To achieve this, we employ the PAW method in
conjunction with real-space methods obtaining a mean absolute error of 10 meV
per atom when compared to reference all-electron values [338].

While OFDFT functionals perform better using local pseudopotentials in bulk
materials, the OFDFT PAW implementation can be interesting for assessing
density functionals. For example, in studying large-$Z$ limits or
semiclassical limits of density functionals, the all-electron values allow us
to find highly-performing OFDFT functionals [339].

###  IX.12 Zero-field splitting

The zero-field splitting (ZFS) refers to the energetic splitting of the
magnetic sub-levels of a localized triplet state in the absence of a magnetic
field [340]. The origin of the ZFS is the magnetic dipole-dipole interactions
between the two electrons of the triplet. This interaction is described by a
spin Hamiltonian of the form ($\alpha,\beta=x,y,z$) [341, 342]

$\hat{H}_{\mathrm{ZFS}}=\sum_{\alpha\beta}\hat{S}_{a}D_{\alpha\beta}\hat{S}_{\beta},$
(121)

where $\mathbf{\hat{S}}$ is the total spin operator and $\mathbf{D}$ is the
ZFS tensor given by

$\displaystyle D_{\alpha\beta}$
$\displaystyle=\frac{1}{2}\frac{\mu_{0}}{4\pi}g_{e}^{2}\mu_{\mathrm{B}}^{2}\int\rho_{2}(\mathbf{r}_{1},\mathbf{r}_{2})\frac{\delta_{\alpha\beta}{r}^{2}-3r_{\alpha}r_{\beta}}{r^{5}}\mathrm{d}\mathbf{r}_{1}\mathrm{d}\mathbf{r}_{2}.$
(122)

where $r_{\alpha}$ and $r_{\beta}$ denote the Cartesian components of
$\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$, $\rho_{2}$ is the two-particle
density matrix of the Kohn–Sham ground-state Slater determinant,
$\mu_{\mathrm{B}}$ is the Bohr magneton and $g_{e}$ is the Landé splitting
factor. GPAW computes the $\mathbf{D}$-tensor by evaluating the double
integral in reciprocal space using the pseudo density including compensation
charges [305].

###  IX.13 Hyperfine coupling

The hyperfine coupling describes the interaction between a magnetic dipole of
a nuclear spin, $\mathbf{\hat{I}}^{N}$, and the magnetic dipole of the
electron-spin distribution, $\mathbf{\hat{S}}(\mathbf{r})$. The interaction is
described by the spin Hamiltonian ($\alpha,\beta=x,y,z$)

$\hat{H}_{\mathrm{HF}}^{N}=\sum_{\alpha\beta}\hat{S}_{\alpha}A_{\alpha\beta}^{N}\hat{I}^{N}_{\beta}$
(123)

where the hyperfine tensor of nucleus $N$ at $\mathbf{R}=\mathbf{0}$ is given
by [343]

$\displaystyle A^{N}_{\alpha\beta}$
$\displaystyle=\frac{2\mu_{0}}{3}g_{e}\mu_{B}g_{N}\mu_{N}\int\delta_{T}(\mathbf{r})\rho_{s}(\mathbf{r})\mathrm{d}\mathbf{r}$
(124)
$\displaystyle+\frac{\mu_{0}}{4\pi}g_{e}\mu_{B}g_{N}\mu_{N}\int\frac{3r_{\alpha}r_{\beta}-\delta_{\alpha\beta}r^{2}}{r^{5}}\rho_{s}(\mathbf{r})\mathrm{d}\mathbf{r}.$

The first term is the isotropic Fermi contact term, which is proportional to
the spin density, $\rho_{s}(\mathbf{r})$, at the centre of the nucleus.
$\delta_{T}(\mathbf{r})$ is a smeared-out $\delta$-function. $g_{e}$ and
$g_{N}$ are the gyromagnetic ratios for the electron and nucleus, and
$\mu_{N}$ is the nuclear magneton. The second term represents the anisotropic
part of the hyperfine coupling tensor and results from dipole-dipole
interactions between nuclear and electronic magnetic moments. GPAW evaluates
$A^{N}$ using the pseudo spin density with compensation charges [305].

##  X Outlook

As described in this review, GPAW is a highly versatile code that is both
maintenance-, user-, and developer-friendly at the same time. The continued
expansion of the code requires substantial efforts and can be lifted only
because of the dedicated team of developers contributing at all levels. There
are currently a number of ongoing as well as planned developments for GPAW,
which will further improve performance and applicability of the code. We are
currently finishing a major refactoring of the code, which will make it even
more developer friendly and facilitate easier implementation of new
functionality.

Another priority is to improve the parallelization of hybrid functional
calculations in plane wave mode by enabling parallelization over bands and
$\mathbf{k}$-points. In the same vein, there is ongoing work to support LCAO-
based hybrid functional calculations using a resolution-of-identity approach.
A natural next step would then be LCAO-based GW calculations. Such a method
could potentially be very efficient compared to plane wave calculations, but
it is currently unclear if the accuracy can be maintained with the limited
LCAO basis. In relation to quasiparticle calculations, there are plans to
implement (quasiparticle) self-consistent GW and vertex corrected GW using
nonlocal xc-kernels from TDDFT. Constrained RPA calculations that provide a
partially screened Coulomb interaction useful for ab initio calculation of
interaction parameters in low-energy model Hamiltonians is currently being
implemented.

Underlying any GPAW calculation are the PAW potentials. The current potentials
date back to 2009. A new set of potentials, including both soft and norm-
conserving potentials (for response function calculations), is being worked
on.

As described herein, GPAW already has an efficient implementation of real time
TDDFT in the LCAO basis while Ehrenfest dynamics is supported only the
comparatively slower in grid mode. Work to enable Ehrenfest dynamics in LCAO
mode is ongoing.

The current version of GPAW supports GPU acceleration only for standard ground
state calculations. The CuPy library greatly simplifies the task of porting
GPAW to GPU and we foresee that large parts of the code, including more
advanced features such as linear response and GW calculations, will become GPU
compatible.

##  XI Acknowledgements

K. S. T. acknowledges funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation program Grant No.
773122 (LIMA) and Grant agreement No. 951786 (NOMAD CoE). K. S. T. is a Villum
Investigator supported by VILLUM FONDEN (grant no. 37789). Funding for A. O.
D. and G. L. was provided by the Icelandic Research Fund (grant nos. 196279
and 217734, respectively). F. N has received funding from the European Union’s
Horizon 2020 research and innovation program under the Marie Skłodowska-Curie
grant agreement No. 899987\. M. M. M. was supported by the Academy of Finland
(grant #338228). T. O. acknowledges support from the Villum foundation Grant
No. 00029378. S. K. and K. W. J. acknowledge support from the VILLUM Center
for Science of Sustainable Fuels and Chemicals, which is funded by the VILLUM
Fonden research grant (9455). T. B. was funded by the Danish National Research
Foundation (DNRF 146). J. S. acknowledges funding from the Independent
Research Fund Denmark (DFF- FTP) through grant no. 9041-00161B. C. S.
acknowledges support from the Swedish Research Council (VR) through Grant No.
2016-06059 and funding from the Horizon Europe research and innovation program
of the European Union under the Marie Skłodowska-Curie grant agreement no.
101065117. Partially funded by the European Union. Views and opinions
expressed are, however, those of the author(s) only and do not necessarily
reflect those of the European Union or REA. Neither the European Union nor the
granting authority can be held responsible for them. T. S. received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (Grant agreement No. 756277-ATMEN). O.
L.-A. have been supported by Minciencias and University of Antioquia
(Colombia). K. T. W. was supported by the U.S. Department of Energy, Office of
Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and
Biosciences Division, Catalysis Science Program to the SUNCAT Center for
Interface Science and Catalysis. G. K. acknowledges funding from V-Sustain:
The VILLUM Centre for the Science of Sustainable Fuels and Chemicals (grant
no. 9455). Additional funding: Knut and Alice Wallenberg Foundation
(2019.0140; J. F. and P. E.), the Swedish Research Council (2020-04935; J. F.
and P. E.), the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No 838996 (T. R.)
and 101065117 (C. S.). Computations were enabled by resources provided by the
National Academic Infrastructure for Supercomputing in Sweden (NAISS) at NSC,
PDC, and C3SE partially funded by the Swedish Research Council through grant
agreement no. 2022-06725.

## References

  * [1] Pierre Hohenberg and Walter Kohn.  Inhomogeneous electron gas.  Phys. Rev., 136(3B):B864, 1964.
  * [2] Walter Kohn and Lu Jeu Sham.  Self-consistent equations including exchange and correlation effects.  Phys. Rev., 140(4A):A1133, 1965.
  * [3] John P Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple.  Phys. Rev. Lett., 77(18):3865, 1996.
  * [4] Axel D Becke.  Density-functional thermochemistry. iii. the role of exact exchange.  The Journal of Chemical Physics, 98:5648–5652, 1993.
  * [5] Jochen Heyd, Gustavo E Scuseria, and Matthias Ernzerhof.  Hybrid functionals based on a screened coulomb potential.  The Journal of Chemical Physics, 118(18):8207–8215, 2003.
  * [6] Glenn M Torrie and John P Valleau.  Nonphysical sampling distributions in monte carlo free-energy estimation: Umbrella sampling.  J. Comput. Phys., 23(2):187–199, 1977.
  * [7] Petros Souvatzis, Olle Eriksson, MI Katsnelson, and SP Rudin.  Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory.  Phys. Rev. Lett., 100(9):095901, 2008.
  * [8] Mark S. Hybertsen and Steven G. Louie.  Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies.  Phys. Rev. B, 34:5390–5413, Oct 1986.
  * [9] Dorothea Golze, Marc Dvorak, and Patrick Rinke.  The GW compendium: A practical guide to theoretical photoemission spectroscopy.  Front. Chem., 7:377, 2019.
  * [10] Mark E Casida, Christine Jamorski, Kim C Casida, and Dennis R Salahub.  Molecular excitation energies to high-lying bound states from time-dependent density-functional response theory: Characterization and correction of the time-dependent local density approximation ionization threshold.  The Journal of Chemical Physics, 108(11):4439–4449, 1998.
  * [11] Giovanni Onida, Lucia Reining, and Angel Rubio.  Electronic excitations: density-functional versus many-body green’s-function approaches.  Rev. Mod. Phys., 74(2):601, 2002.
  * [12] Jens Kehlet Nørskov, Thomas Bligaard, Jan Rossmeisl, and Claus Hviid Christensen.  Towards the computational design of solid catalysts.  Nat. Chem., 1(1):37–46, 2009.
  * [13] Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, et al.  Commentary: The materials project: A materials genome approach to accelerating materials innovation.  APL Mater., 1(1):011002, 2013.
  * [14] Stefano Curtarolo, Gus LW Hart, Marco Buongiorno Nardelli, Natalio Mingo, Stefano Sanvito, and Ohad Levy.  The high-throughput highway to computational materials design.  Nat. Mater., 12(3):191–201, 2013.
  * [15] Sten Haastrup, Mikkel Strange, Mohnish Pandey, Thorsten Deilmann, Per S Schmidt, Nicki F Hinsche, Morten N Gjerding, Daniele Torelli, Peter M Larsen, Anders C Riis-Jensen, et al.  The computational 2d materials database: high-throughput modeling and discovery of atomically thin crystals.  2D Mater., 5(4):042002, 2018.
  * [16] Ivano E Castelli, Thomas Olsen, Soumendu Datta, David D Landis, Søren Dahl, Kristian S Thygesen, and Karsten W Jacobsen.  Computational screening of perovskite metal oxides for optimal solar light capture.  Energ. Environ. Sci., 5(2):5814–5819, 2012.
  * [17] Andreas Marek, Volker Blum, Rainer Johanni, Ville Havu, Bruno Lang, Thomas Auckenthaler, Alexander Heinecke, Hans-Joachim Bungartz, and Hermann Lederer.  The elpa library: scalable parallel eigenvalue solutions for electronic structure theory and computational science.  Journal of Physics: Condensed Matter, 26(21):213201, 2014.
  * [18] Peitao Liu, Merzuk Kaltak, Jiří Klimeš, and Georg Kresse.  Cubic scaling GW: Towards fast quasiparticle calculations.  Phys. Rev. B, 94(16):165109, 2016.
  * [19] Anubhav Jain, Shyue Ping Ong, Wei Chen, Bharat Medasani, Xiaohui Qu, Michael Kocher, Miriam Brafman, Guido Petretto, Gian-Marco Rignanese, Geoffroy Hautier, et al.  Fireworks: a dynamic workflow system designed for high-throughput applications.  Concurr. Comp-pract. E., 27(17):5037–5059, 2015.
  * [20] Giovanni Pizzi, Andrea Cepellotti, Riccardo Sabatini, Nicola Marzari, and Boris Kozinsky.  Aiida: automated interactive infrastructure and database for computational science.  Nato. Sc. S. Ss. Iii. C. S., 111:218–230, 2016.
  * [21] Jens Jørgen Mortensen, Morten Gjerding, and Kristian Sommer Thygesen.  MyQueue: Task and workflow scheduling system.  J. Open Source Softw., 5(45):1844, 2020.
  * [22] Morten Gjerding, Thorbjørn Skovhus, Asbjørn Rasmussen, Fabian Bertoldo, Ask Hjorth Larsen, Jens Jørgen Mortensen, and Kristian Sommer Thygesen.  Atomic simulation recipes: A python framework and library for automated workflows.  Nato. Sc. S. Ss. Iii. C. S., 199:110731, 2021.
  * [23] Jonathan Schmidt, Mário RG Marques, Silvana Botti, and Miguel AL Marques.  Recent advances and applications of machine learning in solid-state materials science.  npj Comput. Mater., 5(1):83, 2019.
  * [24] Albert P Bartók, Mike C Payne, Risi Kondor, and Gábor Csányi.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  Phys. Rev. Lett., 104(13):136403, 2010.
  * [25] Matthias Rupp, Alexandre Tkatchenko, Klaus-Robert Müller, and O Anatole Von Lilienfeld.  Fast and accurate modeling of molecular atomization energies with machine learning.  Phys. Rev. Lett., 108(5):058301, 2012.
  * [26] Andrew M Teale, Trygve Helgaker, Andreas Savin, Carlo Adamo, Bálint Aradi, Alexei V Arbuznikov, Paul W Ayers, Evert Jan Baerends, Vincenzo Barone, Patrizia Calaminici, et al.  DFT exchange: sharing perspectives on the workhorse of quantum chemistry and materials science.  Phys. Chem. Chem. Phys., 24(47):28700–28781, 2022.
  * [27] Xinguo Ren, Patrick Rinke, Christian Joas, and Matthias Scheffler.  Random-phase approximation and its applications in computational chemistry and materials science.  J. Mater. Sci., 47:7447–7471, 2012.
  * [28] Jianwei Sun, Adrienn Ruzsinszky, and John P Perdew.  Strongly constrained and appropriately normed semilocal density functional.  Phys. Rev. Lett., 115(3):036402, 2015.
  * [29] Thomas Olsen and Kristian S Thygesen.  Random phase approximation applied to solids, molecules, and graphene-metal interfaces: From van der Waals to covalent bonding.  Phys. Rev. B, 87(7):075111, 2013.
  * [30] Alexandre Tkatchenko and Matthias Scheffler.  Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data.  Phys. Rev. Lett., 102(7):073005, 2009.
  * [31] Max Dion, Henrik Rydberg, Elsebeth Schröder, David C Langreth, and Bengt I Lundqvist.  Van der Waals density functional for general geometries.  Phys. Rev. Lett., 92(24):246401, 2004.
  * [32] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen.  Real-space grid implementation of the projector augmented wave method.  Phys. Rev. B, 71:035109, Jan 2005.
  * [33] J Enkovaara, C Rostgaard, J J Mortensen, J Chen, M Dułak, L Ferrighi, J Gavnholt, C Glinsvad, V Haikola, H A Hansen, H H Kristoffersen, M Kuisma, A H Larsen, L Lehtovaara, M Ljungberg, O Lopez-Acevedo, P G Moses, J Ojanen, T Olsen, V Petzold, N A Romero, J Stausholm-Møller, M Strange, G A Tritsaris, M Vanin, M Walter, B Hammer, H Häkkinen, G K H Madsen, R M Nieminen, J K Nørskov, M Puska, T T Rantala, J Schiøtz, K S Thygesen, and K W Jacobsen.  Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method.  Journal of Physics: Condensed Matter, 22(25):253202, jun 2010.
  * [34] Michael Walter, Hannu Häkkinen, Lauri Lehtovaara, Martti Puska, Jussi Enkovaara, Carsten Rostgaard, and Jens Jørgen Mortensen.  Time-dependent density-functional theory in the projector augmented-wave method.  The Journal of Chemical Physics, 128(24):244101, June 2008.
  * [35] A. H. Larsen, M. Vanin, J. J. Mortensen, K. S. Thygesen, and K. W. Jacobsen.  Localized atomic basis set in the projector augmented wave method.  Phys. Rev. B, 80:195112, Nov 2009.
  * [36] Jun Yan, Jens J Mortensen, Karsten W Jacobsen, and Kristian S Thygesen.  Linear density response function in the projector augmented wave method: Applications to solids, surfaces, and interfaces.  Phys. Rev. B, 83(24):245122, 2011.
  * [37] Thomas Olsen and Kristian S Thygesen.  Extending the random-phase approximation for electronic correlation energies: The renormalized adiabatic local density approximation.  Phys. Rev. B, 86(8):081103, 2012.
  * [38] Falco Hüser, Thomas Olsen, and Kristian S Thygesen.  Quasiparticle gw calculations for solids, molecules, and two-dimensional materials.  Phys. Rev. B, 87(23):235132, 2013.
  * [39] Jun Yan, Karsten W Jacobsen, and Kristian S Thygesen.  Optical properties of bulk semiconductors and graphene/boron nitride: The bethe-salpeter equation with derivative discontinuity-corrected density functional energies.  Phys. Rev. B, 86(4):045208, 2012.
  * [40] Joachim Sødequist and Thomas Olsen.  Type ii multiferroic order in two-dimensional transition metal halides from first principles spin-spiral calculations.  2D Mater., 10(3):035016, may 2023.
  * [41] Daniele Torelli and Thomas Olsen.  Calculating critical temperatures for ferromagnetic order in two-dimensional materials.  2D Mater., 6:015028, 2018.
  * [42] Frederik L Durhuus, Thorbjørn Skovhus, and Thomas Olsen.  Plane wave implementation of the magnetic force theorem for magnetic exchange constants: application to bulk fe, co and ni.  Journal of Physics: Condensed Matter, 35(10):105802, jan 2023.
  * [43] Thorbjørn Skovhus and Thomas Olsen.  Dynamic transverse magnetic susceptibility in the projector augmented-wave method: Application to Fe, Ni, and Co.  Phys. Rev. B, 103:245110, 2021.
  * [44] Mads Kruse, Urko Petralanda, Morten N. Gjerding, Karsten W. Jacobsen, Kristian S. Thygesen, and Thomas Olsen.  Two-dimensional ferroelectrics from high throughput computational screening.  npj Comput. Mater., 9:45, 2023.
  * [45] Morten Niklas Gjerding, Alireza Taghizadeh, Asbjørn Rasmussen, Sajid Ali, Fabian Bertoldo, Thorsten Deilmann, Nikolaj Rørbæk Knøsgaard, Mads Kruse, Ask Hjorth Larsen, Simone Manti, Thomas Garm Pedersen, Urko Petralanda, Thorbjørn Skovhus, Mark Kamper Svendsen, Jens Jørgen Mortensen, Thomas Olsen, and Kristian Sommer Thygesen.  Recent progress of the computational 2D materials database (C2DB).  2D Mater., 8:044002, 2021.
  * [46] Thomas Olsen, Erik Andersen, Takuya Okugawa, Daniele Torelli, Thorsten Deilmann, and Kristian S. Thygesen.  Discovering two-dimensional topological insulators from high-throughput computations.  Phys. Rev. Mater., 3:024005, 2019.
  * [47] Ask Hjorth Larsen, Jens Jørgen Mortensen, Jakob Blomqvist, Ivano E Castelli, Rune Christensen, Marcin Dułak, Jesper Friis, Michael N Groves, Bjørk Hammer, Cory Hargus, Eric D Hermes, Paul C Jennings, Peter Bjerre Jensen, James Kermode, John R Kitchin, Esben Leonhard Kolsbjerg, Joseph Kubal, Kristen Kaasbjerg, Steen Lysgaard, Jón Bergmann Maronsson, Tristan Maxson, Thomas Olsen, Lars Pastewka, Andrew Peterson, Carsten Rostgaard, Jakob Schiøtz, Ole Schütt, Mikkel Strange, Kristian S Thygesen, Tejs Vegge, Lasse Vilhelmsen, Michael Walter, Zhenhua Zeng, and Karsten W Jacobsen.  The atomic simulation environment—a python library for working with atoms.  Journal of Physics: Condensed Matter, 29(27):273002, jun 2017.
  * [48] Arash A Mostofi, Jonathan R Yates, Young-Su Lee, Ivo Souza, David Vanderbilt, and Nicola Marzari.  wannier90: A tool for obtaining maximally-localised wannier functions.  Comput. Phys. Commun., 178(9):685–699, 2008.
  * [49] Yorick L. A. Schmerwitz, Gianluca Levi, and Hannes Jónsson.  Calculations of Excited Electronic States by Converging on Saddle Points Using Generalized Mode Following.  J. Chem. Theory Comput., 19(12):3634–3651, 2023.
  * [50] Aleksei V. Ivanov, Gianluca Levi, Elvar Ö. Jónsson, and Hannes Jónsson.  Method for calculating excited electronic states using density functionals and direct orbital optimization with real space grid or plane-wave basis set.  J. Chem. Theory Comput., 17:5034–5049, 8 2021.
  * [51] Gianluca Levi, Aleksei V. Ivanov, and Hannes Jónsson.  Variational density functional calculations of excited states via direct optimization.  J. Chem. Theory Comput., 16:6968–6982, 11 2020.
  * [52] M. Kuisma, A. Sakko, T. P. Rossi, A. H. Larsen, J. Enkovaara, L. Lehtovaara, and T. T. Rantala.  Localized surface plasmon resonance in silver nanoparticles: Atomistic first-principles time-dependent density-functional theory calculations.  Phys. Rev. B, 91(11):115431, 2015.
  * [53] Tuomas P. Rossi, Susi Lehtola, Arto Sakko, Martti J. Puska, and Risto M. Nieminen.  Nanoplasmonics simulations at the basis set limit through completeness-optimized, local numerical basis sets.  The Journal of Chemical Physics, 142(9):094114, 2015.
  * [54] Tuomas P. Rossi, Mikael Kuisma, Martti J. Puska, Risto M. Nieminen, and Paul Erhart.  Kohn–Sham Decomposition in Real-Time Time-Dependent Density-Functional Theory: An Efficient Tool for Analyzing Plasmonic Excitations.  J. Chem. Theory Comput., 13(10):4779–4790, 2017.
  * [55] Esko Makkonen, Tuomas P Rossi, Ask Hjorth Larsen, Olga Lopez-Acevedo, Patrick Rinke, Mikael Kuisma, and Xi Chen.  Real-time time-dependent density functional theory implementation of electronic circular dichroism applied to nanoscale metal–organic clusters.  The Journal of Chemical Physics, 154(11):114102, 2021.
  * [56] Alireza Taghizadeh, Ulrik Leffers, Thomas G Pedersen, and Kristian S Thygesen.  A library of ab initio raman spectra for automated identification of 2d materials.  Nat. Commun., 11(1):3011, 2020.
  * [57] Alireza Taghizadeh, Kristian S Thygesen, and Thomas G Pedersen.  Two-dimensional materials with giant optical nonlinearities near the theoretical upper limit.  Acs Nano, 15(4):7155–7167, 2021.
  * [58] Mikkel Ohm Sauer, Alireza Taghizadeh, Urko Petralanda, Martin Ovesen, Kristian Sommer Thygesen, Thomas Olsen, Horia Cornean, and Thomas Garm Pedersen.  Shift current photovoltaic efficiency of 2d materials.  npj Comput. Mater., 9(1):35, 2023.
  * [59] GPAW’s git-repository: https://gitlab.com/gpaw/gpaw/.
  * [60] Holger Krekel, Bruno Oliveira, Ronny Pfannschmidt, Floris Bruynooghe, Brianna Laugher, and Florian Bruhin.  pytest, 2004.
  * [61] GPAW tutorials and exercises: https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/tutorialsexercises.html.
  * [62] Matteo Frigo and Steven G. Johnson.  The design and implementation of FFTW3.  P. Ieee, 93(2):216–231, 2005.  Special issue on “Program Generation, Optimization, and Platform Adaptation”.
  * [63] L. S. Blackford, J. Choi, A. Cleary, E. D’Azevedo, J. Demmel, I. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker, and R. C. Whaley.  ScaLAPACK Users’ Guide.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 1997\.
  * [64] Miguel A.L. Marques, Micael J.T. Oliveira, and Tobias Burnus.  Libxc: A library of exchange and correlation functionals for density functional theory.  Comput. Phys. Commun., 183(10):2272–2281, 2012.
  * [65] Susi Lehtola, Conrad Steigemann, Micael J. T. Oliveira, and Miguel A. L. Marques.  Recent developments in libxc — A comprehensive library of functionals for density functional theory.  SoftwareX, 7:1–5, January 2018.
  * [66] Ask Hjorth Larsen, Mikael Kuisma, Joakim Löfgren, Yann Pouillon, Paul Erhart, and Per Hyldgaard.  libvdwxc: a library for exchange–correlation functionals in the vdw-df family.  Model. Simul. Mater. Sc., 25(6):065004, jun 2017.
  * [67] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [68] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nat. Methods, 17:261–272, 2020.
  * [69] P. E. Blöchl.  Projector augmented-wave method.  Phys. Rev. B, 50:17953–17979, Dec 1994.
  * [70] Toma Susi, Jacob Madsen, Ursula Ludacka, Jens Jørgen Mortensen, Timothy J. Pennycook, Zhongbo Lee, Jani Kotakoski, Ute Kaiser, and Jannik C. Meyer.  Efficient first principles simulation of electron scattering factors for transmission electron microscopy.  Ultramicroscopy, 197:16–22, February 2019.
  * [71] Jacob Madsen, Timothy J. Pennycook, and Toma Susi.  ab initio description of bonding for transmission electron microscopy.  Ultramicroscopy, 231:113253, December 2021.
  * [72] Jacob Madsen and Toma Susi.  abTEM: transmission electron microscopy from first principles.  Open Research Europe, 1(24):13015, 2021.
  * [73] G. Kresse and J. Furthmüller.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  Phys. Rev. B, 54:11169–11186, Oct 1996.
  * [74] Nicola Marzari, David Vanderbilt, Alessandro De Vita, and M. C. Payne.  Thermal contraction and disordering of the al(110) surface.  Phys. Rev. Lett., 82:3296–3299, Apr 1999.
  * [75] Peter E. Blöchl, O. Jepsen, and O. K. Andersen.  Improved tetrahedron method for brillouin-zone integrations.  Phys. Rev. B, 49:16223–16233, Jun 1994.
  * [76] Susi Lehtola, Frank Blockhuys, and Christian Van Alsenoy.  An overview of self-consistent field calculations within finite basis sets.  Molecules, 25(5):1–23, 2020.
  * [77] Troy Van Voorhis and Martin Head-Gordon.  A geometric approach to direct minimization.  Mol. Phys., 100:1713–1721, 2002.
  * [78] M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos.  Iterative minimization techniques for ab initio total-energy calculations: Molecular dynamics and conjugate gradients.  Rev. Mod. Phys., 64(4):1045–1097, 1992.
  * [79] Martin Head-Gordon and John A. Pople.  Optimization of wave function and geometry in the finite basis Hartree-Fock method.  J. Phys. Chem-us., 92(11):3063–3069, 1988.
  * [80] Aleksei V. Ivanov, Elvar Ö. Jónsson, Tejs Vegge, and Hannes Jónsson.  Direct energy minimization based on exponential transformation in density functional calculations of finite and extended systems.  Comput. Phys. Commun., 267:108047, 2021.
  * [81] PAW-XML specification.  https://esl.cecam.org/data/paw-xml/.
  * [82] Xavier Gonze, Bernard Amadon, Gabriel Antonius, Frédéric Arnardi, Lucas Baguet, Jean-Michel Beuken, Jordan Bieder, François Bottin, Johann Bouchet, Eric Bousquet, Nils Brouwer, Fabien Bruneval, Guillaume Brunin, Théo Cavignac, Jean-Baptiste Charraud, Wei Chen, Michel Côté, Stefaan Cottenier, Jules Denier, Grégory Geneste, Philippe Ghosez, Matteo Giantomassi, Yannick Gillet, Olivier Gingras, Donald R. Hamann, Geoffroy Hautier, Xu He, Nicole Helbig, Natalie Holzwarth, Yongchao Jia, François Jollet, William Lafargue-Dit-Hauret, Kurt Lejaeghere, Miguel A.L. Marques, Alexandre Martin, Cyril Martins, Henrique P.C. Miranda, Francesco Naccarato, Kristin Persson, Guido Petretto, Valentin Planes, Yann Pouillon, Sergei Prokhorenko, Fabio Ricci, Gian-Marco Rignanese, Aldo H. Romero, Michael Marcus Schmitt, Marc Torrent, Michiel J. van Setten, Benoit Van Troeye, Matthieu J. Verstraete, Gilles Zérah, and Josef W. Zwanziger.  The ABINIT project: Impact, environment and recent developments.  Comput. Phys. Commun., 248:107042, 2020.
  * [83] François Jollet, Marc Torrent, and Natalie Holzwarth.  Generation of projector augmented-wave atomic data: A 71 element validated table in the XML format.  Comput. Phys. Commun., 185(4):1246–1254, 2014.
  * [84] C. Hartwigsen, S. Goedecker, and J. Hutter.  Relativistic separable dual-space gaussian pseudopotentials from h to rn.  Phys. Rev. B, 58:3641–3662, Aug 1998.
  * [85] Martin Schlipf and François Gygi.  Optimization algorithm for the generation of oncv pseudopotentials.  Comput. Phys. Commun., 196:36–44, 2015.
  * [86] CuPy.  https://cupy.dev/.
  * [87] Header only porting.  https://github.com/mlouhivu/hop/.
  * [88] Samuli Hakala, Ville Havu, Jussi Enkovaara, and Risto M. Nieminen.  Parallel electronic structure calculations using multiple graphics processing units (gpus).  Lect. Notes Comput. Sc., 7782:63–76, 2013.
  * [89] Samuli Hakala, Jussi Enkovaara, Ville Havu, Jun Yan, Lin Li, Chris O’Grady, and Risto M. Nieminen.  Grid-based projector-augmented wave method.  Method in Electronic Structure Calculations on Graphics Processing Units, 2015.
  * [90] John P. Perdew and Karla Schmidt.  Jacob’s ladder of density functional approximations for the exchange-correlation energy.  In AIP Conference Proceedings, volume 577, pages 1–20. AIP Publishing, July 2001.
  * [91] Jianmin Tao, John P. Perdew, Viktor N. Staroverov, and Gustavo E. Scuseria.  Climbing the Density Functional Ladder: Nonempirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids.  Phys. Rev. Lett., 91(14):146401, September 2003.
  * [92] John P. Perdew and Yue Wang.  Accurate and simple analytic representation of the electron-gas correlation energy.  Phys. Rev. B, 45(23):13244–13249, June 1992.
  * [93] Guillermo Román-Pérez and José M. Soler.  Efficient implementation of a van der waals density functional: Application to double-wall carbon nanotubes.  Phys. Rev. Lett., 103:096102, Aug 2009.
  * [94] Oleg Gritsenko, Robert van Leeuwen, Erik van Lenthe, and Evert Jan Baerends.  Self-consistent approximation to the Kohn-Sham exchange potential.  Phys. Rev. A, 51(3):1944, 1995.
  * [95] M. Kuisma, J. Ojanen, J. Enkovaara, and T. T. Rantala.  Kohn-Sham potential with discontinuity for band gap materials.  Phys. Rev. B, 82(11):115106, 2010.
  * [96] Jun Yan, Karsten W. Jacobsen, and Kristian S. Thygesen.  First-principles study of surface plasmons on Ag(111) and H/Ag(111).  Phys. Rev. B, 84(23):235430, 2011.
  * [97] Jun Yan, Karsten W. Jacobsen, and Kristian S. Thygesen.  Conventional and acoustic surface plasmons on noble metal surfaces: A time-dependent density functional theory study.  Phys. Rev. B, 86(24):241404, 2012.
  * [98] J. Magnus Rahm, Christopher Tiburski, Tuomas P. Rossi, Ferry Anggoro Ardy Nugroho, Sara Nilsson, Christoph Langhammer, and Paul Erhart.  A library of late transition metal alloy dielectric functions for nanophotonic applications.  Adv. Funct. Mater., 30(35):2002122, 2020.
  * [99] Keshab Bashyal, Christopher K Pyles, Sajjad Afroosheh, Aneer Lamichhane, and Alexey T Zayak.  Empirical optimization of DFT+U and HSE for the band structure of ZnO.  Journal of Physics: Condensed Matter, 30(6):065501, jan 2018.
  * [100] M. K. Y. Chan and G. Ceder.  Efficient band gap prediction for solids.  Phys. Rev. Lett., 105:196403, Nov 2010.
  * [101] John P. Perdew.  Density functional theory and the band gap problem.  Int. J. Quantum Chem., 28(S19):497–523, 1985.
  * [102] Kevin J. May and Alexie M. Kolpak.  Improved description of perovskite oxide crystal structure and electronic properties using self-consistent Hubbard $u$ corrections from ACBN0.  Phys. Rev. B, 101:165117, Apr 2020.
  * [103] 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, Jan 1998.
  * [104] 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–R5470, Aug 1995.
  * [105] F Aryasetiawan, T Miyake, and R Sakuma.  The constrained RPA method for calculating the Hubbard U from first-principles.  The LDA+ DMFT approach to strongly correlated materials, 2011\.
  * [106] Matteo Cococcioni and Stefano de Gironcoli.  Linear response approach to the calculation of the effective interaction parameters in the LDA+U method.  Phys. Rev. B, 71(3), jan 2005.
  * [107] Guy C. Moore, Matthew K. Horton, Alexander M. Ganose, Martin Siron, Edward Linscott, David D. O’Regan, and Kristin A. Persson.  High-throughput determination of Hubbard U and Hund J values for transition metal oxides via linear response formalism, 2022.
  * [108] Lei Wang, Thomas Maxisch, and Gerbrand Ceder.  Oxidation energies of transition metal oxides within the $\mathrm{GGA}+\mathrm{U}$ framework.  Phys. Rev. B, 73:195107, May 2006.
  * [109] M F M Taib, D T Mustaffa, N H Hussin, M H Samat, A M M Ali, O H Hassan, and M Z A Yahya.  First principles study on Zn doped MgO using Hubbard U correction.  Materials Research Express, 6(9):094012, jul 2019.
  * [110] Maituo Yu, Shuyang Yang, Chunzhi Wu, and Noa Marom.  Machine learning the Hubbard U parameter in DFT+U using Bayesian optimization.  npj Comput. Mater., 6(1):180, 2020.
  * [111] Yoshihiro Tawada, Takao Tsuneda, Susumu Yanagisawa, Takeshi Yanai, and Kimihiko Hirao.  A long-range-corrected time-dependent density functional theory.  The Journal of Chemical Physics, 120(18):8425–8433, May 2004.
  * [112] Aron J. Cohen, Paula Mori-Sánchez, and Weitao Yang.  Insights into Current Limitations of Density Functional Theory.  Science, 321(5890):792–794, August 2008.
  * [113] Andreas Dreuw, Jennifer L. Weisman, and Martin Head-Gordon.  Long-range charge-transfer excited states in time-dependent density functional theory require non-local exchange.  J. Chem. Phys., 119(6):2943–2946, 2003.
  * [114] E. J. Baerends, O. V. Gritsenko, and R. van Meer.  The Kohn–Sham gap, the fundamental gap and the optical gap: The physical meaning of occupied and virtual Kohn–Sham orbital energies.  Phys. Chem. Chem. Phys., 15(39):16408–16425, September 2013.
  * [115] Stephan Kümmel.  Charge-Transfer Excitations: A Challenge for Time-Dependent Density Functional Theory That Has Been Met.  Adv. Energy Mater., 7(16):1700440, August 2017.
  * [116] Roi Baer, Ester Livshits, and Ulrike Salzner.  Tuned Range-Separated Hybrids in Density Functional Theory.  Annu. Rev. Phys. Chem., 61(1):85–109, 2010.
  * [117] Carlo Adamo and Vincenzo Barone.  Toward chemical accuracy in the computation of NMR shieldings: The PBE0 model.  Chem. Phys. Lett., 298(1–3):113–119, December 1998.
  * [118] Takeshi Yanai, David P Tew, and Nicholas C Handy.  A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP).  Chem. Phys. Lett., 393(1–3):51–57, July 2004.
  * [119] Yoshinobu Akinaga and Seiichiro Ten-no.  Range-separation by the Yukawa potential in long-range corrected density functional theory with Gaussian-type basis functions.  Chem. Phys. Lett., 462(4–6):348–351, September 2008\.
  * [120] Michael Seth and Tom Ziegler.  Range-Separated Exchange Functionals with Slater-Type Functions.  J. Chem. Theory Comput., 8(3):901–907, March 2012.
  * [121] Ester Livshits and Roi Baer.  A well-tempered density functional theory of electrons in molecules.  Phys. Chem. Chem. Phys., 9(23):2932–2941, June 2007.
  * [122] Rolf Würdemann.  Berechnung optischer Spektren und Grundzustandseigenschaften neutraler und geladener Moleküle mittels Dichtefunktionaltheorie.  PhD thesis, Universität Freiburg, Freiburg, Germany, July 2016.
  * [123] Rolf Würdemann and Michael Walter.  Charge Transfer Excitations with Range Separated Functionals Using Improved Virtual Orbitals.  J. Chem. Theory Comput., 14(7):3667–3676, July 2018.
  * [124] Perdew J. P. and A. Zunger.  Self-interaction correction to density-functional approximations for many-electron systems.  Phys. Rev. B, 23:5048, 1981.
  * [125] S. Klüpfel, P. J. Klüpfel, and H. Jónsson.  Importance of complex orbitals in calculating the self-interaction corrected ground state of atoms.  Phys. Rev. A, 84:050501(R), 2011.
  * [126] S. Lehtola and H. Jónsson.  Variational, self-consistent implementation of the perdew-zunger self-interaction correction with complex optimal orbitals.  J. Chem. Theory Comput., 10:5324, 2014.
  * [127] S. Lehtola, M. Head-Gordon, and H. Jónsson.  Complex orbitals, multiple local minima, and symmetry breaking in perdew-zunger self-interaction corrected density functional theory calculations.  J. Chem. Theory Comput., 12:3195, 2016.
  * [128] H. Jónsson.  Simulation of surface processes.  Proceedings of the National Academy of Sciences, 108:944, 2011.
  * [129] S. Klüpfel, P. Klüpfel, and H. Jónsson.  The effect of the perdew-zunger self-interaction correction to density functionals on the energetics of small molecules.  J. Chem. Phys., 137:124102, 2012.
  * [130] H. Gudmundsdóttir, Y. Zhang, P. M. Weber, and H. Jónsson.  Self-interaction corrected density functional calculations of molecular rydberg states.  J. Chem. Phys., 139:194102, 2013.
  * [131] A. V. Ivanov, T. Ghosh, E. Ö. Jónsson, and H. Jónsson.  Mn dimer can be described accurately with density functional calculations when self-interaction correction is applied.  J. Phys. Chem. Lett., 12:4240, 2021.
  * [132] H. Gudmundsdóttir, E. Ö. Jónsson, and H. Jónsson.  Calculations of al dopant in alpha-quartz using a variational implementation of the perdew-zunger self-interaction correction.  New J. Phys., 17:083006, 2015.
  * [133] Yorick L. A. Schmerwitz, Aleksei V. Ivanov, Elvar Ö. Jónsson, Hannes Jónsson, and Gianluca Levi.  Variational density functional calculations of excited states: Conical intersection and avoided crossing in ethylene bond twisting.  J. Phys. Chem. Lett., 13:3990–3999, 2022.
  * [134] J. J. Mortensen, K. Kaasbjerg, S. L. Frederiksen, J. K. Nørskov, J. P. Sethna, and K. W. Jacobsen.  Bayesian Error Estimation in Density-Functional Theory.  Phys. Rev. Lett., 95(21):216401, November 2005.
  * [135] Jess Wellendorff, Keld T. Lundgaard, Andreas Møgelhøj, Vivien Petzold, David D. Landis, Jens K. Nørskov, Thomas Bligaard, and Karsten W. Jacobsen.  Density functionals for surface science: Exchange-correlation model development with Bayesian error estimation.  Phys. Rev. B, 85:235149, 2012.
  * [136] Jess Wellendorff, Keld T. Lundgaard, Karsten W. Jacobsen, and Thomas Bligaard.  mBEEF: An accurate semi-local Bayesian error estimation density functional.  The Journal of Chemical Physics, 140(14):144107, April 2014.
  * [137] Keld T. Lundgaard, Jess Wellendorff, J. Voss, Karsten W. Jacobsen, and Thomas Bligaard.  mBEEF-vdW: Robust fitting of error estimation density functionals.  Phys. Rev. B, 93(23):235162, 2016.
  * [138] Andrew J. Medford, Jess Wellendorff, Aleksandra Vojvodic, Felix Studt, Frank Abild-Pedersen, Karsten W. Jacobsen, Thomas Bligaard, and Jens K. Nørskov.  Assessing the reliability of calculated catalytic ammonia synthesis rates.  Science, 345(6193):197–200, 2014.
  * [139] G Ciccotti, M Ferrario, and J P Ryckaert.  Molecular dynamics of rigid systems in cartesian coordinates a general formulation.  Mol. Phys., 47:1253–1264, 1982.
  * [140] Andrew A Peterson.  Global optimization of adsorbate–surface structures while preserving molecular identity.  Top. Catal., 57:40–53, 2014.
  * [141] E B Tadmor, G S Smith, N Bernstein, and E Kaxiras.  Mixed finite element and atomistic formulation for complex crystals.  Phys. Rev. B, 59:235, 1999.
  * [142] C G Broyden.  The convergence of a class of double-rank minimization algorithms.  J. I. Math. Appl., 6:76–90, 1970.
  * [143] R Fletcher.  A new approach to variable metric algorithms.  Comput. J., 13:317–322, 1970.
  * [144] D Goldfarb.  A family of variable metric updates derived by variational means.  Math. Comput., 24:23–26, 1970.
  * [145] D F Shanno.  Conditioning of quasi-newton methods for function minimization.  Math. Comput., 24:647–656, 1970.
  * [146] D Liu and J Nocedal.  On the limited memory method for large scale optimization.  Mathematical Programming B, 45:503–528, 1989.
  * [147] Erik Bitzek, Pekka Koskinen, Franz Gahler, Michael Moseler, and Peter Gumbsch.  Structural relaxation made simple.  Phys. Rev. Lett., 97:170201, 2006.
  * [148] D Packwood, J R Kermode, L Mones, N Bernstein, J Woolley, N Gould, C Ortner, and G Csanyi.  A universal preconditioner for simulating condensed phase materials.  J. Chem. Phys., 144:164109, 2016.
  * [149] Estefania Garijo del Rio, Jens J Mortensen, and Karsten W Jacobsen.  Local bayesian optimizer for atomic structures.  Phys. Rev. B, 100:104103, 2019.
  * [150] G Henkelman and H Jonsson.  A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives.  J. Chem. Phys., 111:7010, 1999.
  * [151] H Jonsson, G Mills, and K W Jacobsen.  in ‘classical and quantum dynamics in condensed phase systems’, edited by b. j. berne, g. cicotti, and d. f. coker.  World Scientific, 1998.
  * [152] G Henkelman and H Jonsson.  Improved tangent estimate in the neb method for finding minimum energy paths and saddle points.  J. Chem. Phys., 113:9978, 2000.
  * [153] G Henkelman, B P Uberuaga, and H Jonsson.  A climbing-image neb method for finding saddle points and minimum energy paths.
  * [154] S Smidstrup, A Pedersen, K Stokbro, and H Jonsson.  Improved initial guess for minimum energy path calculations.  J. Chem. Phys., 140:214106, 2014.
  * [155] E L Kolsbjerg, M N Groves, and B Hammer.  An automated nudged elastic band method.  J. Chem. Phys., 145:094107, 2016.
  * [156] P Lindgren, G Kastlunger, and A A Peterson.  Scaled and dynamic optimizations of nudged elastic bands.  J. Chem. Theory Comput., 15:5787–5793, 2019.
  * [157] S Makri, C Ortner, and J R Kermode.  A preconditioning scheme for minimum energy path finding methods.  J. Chem. Phys., 150:094109, 2019.
  * [158] A A Peterson.  Acceleration of saddle-point searches with machine learning.  J. Chem. Phys., 145:074106, 2016.
  * [159] O Koistinen, F B Dagbjartsdottir, V Asgeirsson, A Vehtari, and H Jonsson.  Nudged elastic band calculations accelerated with gaussian process regression.  J. Chem. Phys., 147:152720, 2017.
  * [160] J A G Torres, P C Jennings, M H Hansen, J R Boes, and T Bligaard.  Low–scaling algorithm for nudged elastic band calculations using a surrogate machine learning model.  Phys. Rev. Lett., 122:156001, 2019.
  * [161] Manuel Šarić, Jan Rossmeisl, and Poul Georg Moses.  Modeling the adsorption of sulfur containing molecules and their hydrodesulfurization intermediates on the Co-promoted MoS2 catalyst by DFT.  J. Catal., 358:131–140, February 2018.
  * [162] Thomas M. Østergaard, Livia Giordano, Ivano E. Castelli, Filippo Maglia, Byron K. Antonopoulos, Yang Shao-Horn, and Jan Rossmeisl.  Oxidation of ethylene carbonate on li metal oxide surfaces.  The Journal of Physical Chemistry C, 122(19):10442–10449, April 2018.
  * [163] Logi Arnarson, Per S. Schmidt, Mohnish Pandey, Alexander Bagger, Kristian S. Thygesen, Ifan E. L. Stephens, and Jan Rossmeisl.  Fundamental limitation of electrocatalytic methane conversion to methanol.  Phys. Chem. Chem. Phys., 20(16):11152–11159, 2018.
  * [164] Hao Wan, Anders W. Jensen, María Escudero-Escribano, and Jan Rossmeisl.  Insights in the oxygen reduction reaction: From metallic electrocatalysts to diporphyrins.  ACS Catal., 10(11):5979–5989, April 2020.
  * [165] Gianluca Levi, Elisa Biasin, Asmus Ougaard Dohn, and Hannes Jónsson.  On the interplay of solvent and conformational effects in simulated excited-state dynamics of a copper phenanthroline photosensitizer.  Phys. Chem. Chem. Phys., 22:748–757, 2020.
  * [166] D A Wales and J P K Doye.  Global optimization of clusters, crystals, and biomolecules.  Science, 285:1368, 1999.
  * [167] S Goedecker.  Minima hopping: An efficient search method for the global minimum of the potential energy surface of complex molecular systems.  J. Chem. Phys., 120:9911, 2004.
  * [168] L B Vilhelmsen and B Hammer.  A genetic algorithm for first principles global structure optimization of supported nano structures.  The Journal of Chemical Physics, 141:044711, 2014.
  * [169] Mads-Peter V. Christiansen, Nikolaj Rønne, and Bjørk Hammer.  Atomistic global optimization x: A python package for optimization of atomistic structures.  The Journal of Chemical Physics, 157(5):054701, August 2022.
  * [170] Paul C. Jennings, Steen Lysgaard, Jens Strabo Hummelshøj, Tejs Vegge, and Thomas Bligaard.  Genetic algorithms for computational materials discovery accelerated by machine learning.  npj Comput. Mater., 5(1), April 2019.
  * [171] Malthe K. Bisbo and Bjørk Hammer.  Efficient Global Structure Optimization with a Machine-Learned Surrogate Model.  Phys. Rev. Lett., 124(8):086102, February 2020.
  * [172] Sami Kaappa, Casper Larsen, and Karsten Wedel Jacobsen.  Atomic structure optimization with machine-learning enabled interpolation between chemical elements.  Phys. Rev. Lett., 127(16), October 2021.
  * [173] Sami Kaappa, Estefanía Garijo del Río, and Karsten Wedel Jacobsen.  Global optimization of atomic structures with gradient-enhanced gaussian process regression.  Phys. Rev. B, 103(17), May 2021.
  * [174] Casper Larsen, Sami Kaappa, Andreas Lynge Vishart, Thomas Bligaard, and Karsten Wedel Jacobsen.  Machine-learning enabled optimization of atomic structures using atoms with fractional existence, 2022.
  * [175] Ralf Wanzenböck, Marco Arrigoni, Sebastian Bichelmaier, Florian Buchner, Jesús Carrete, and Georg K. H. Madsen.  Neural-network-backed evolutionary search for SrTiO$\\_3$(110) surface reconstructions.  Digital Discovery, 1(5):703–710, 2022.
  * [176] Andreas Møgelhøj, André K. Kelkkanen, K. Thor Wikfeldt, Jakob Schiøtz, Jens Jørgen Mortensen, Lars G. M. Pettersson, Bengt I. Lundqvist, Karsten W. Jacobsen, Anders Nilsson, and Jens K. Nørskov.  Ab initio van der Waals interactions in simulations of water alter structure from mainly tetrahedral to high-density-like.  The Journal of Physical Chemistry B, 115(48):14149–14160, August 2011.
  * [177] Martin Hangaard Hansen and Jan Rossmeisl.  pH in grand canonical statistics of an electrochemical interface.  The Journal of Physical Chemistry C, 120(51):29135–29143, December 2016.
  * [178] A. O. Dohn, E. Ö. Jónsson, G. Levi, J. J. Mortensen, O. Lopez-Acevedo, K. S. Thygesen, K. W. Jacobsen, J. Ulstrup, N. E. Henriksen, K. B. Møller, and H. Jónsson.  Grid-based projector augmented wave (GPAW) implementation of quantum mechanics/molecular mechanics (QM/MM) electrostatic embedding and application to a solvated diplatinum complex.  J. Chem. Theory Comput., 13(12):6010–6022, 2017.
  * [179] Asmus Ougaard Dohn, Elvar Örn Jónsson, Kasper Skov Kjær, Tim Brandt van Driel, Martin Meedom Nielsen, Karsten Wedel Jacobsen, Niels Engholm Henriksen, and Klaus Braagaard Møller.  Direct dynamics studies of a binuclear metal complex in solution: The interplay between vibrational relaxation, coherence, and solvent effects.  J. Phys. Chem. Lett., 5(14):2414–2418, June 2014.
  * [180] Tim B. van Driel, Kasper S. Kjær, Robert W. Hartsock, Asmus O. Dohn, Tobias Harlang, Matthieu Chollet, Morten Christensen, Wojciech Gawelda, Niels E. Henriksen, Jong Goo Kim, Kristoffer Haldrup, Kyung Hwan Kim, Hyotcherl Ihee, Jeongho Kim, Henrik Lemke, Zheng Sun, Villy Sundström, Wenkai Zhang, Diling Zhu, Klaus B. Møller, Martin M. Nielsen, and Kelly J. Gaffney.  Atomistic characterization of the active-site solvation dynamics of a model photocatalyst.  Nat. Commun., 7(1), 2016.
  * [181] Gianluca Levi, Mátyás Pápai, Niels E. Henriksen, Asmus O. Dohn, and Klaus B. Møller.  Solution structure and ultrafast vibrational relaxation of the ptpop complex revealed by $\delta$scf-qm/mm direct dynamics simulations.  J. Phys. Chem. C, 122:7100–7119, 4 2018.
  * [182] Kristoffer Haldrup, Gianluca Levi, Elisa Biasin, Peter Vester, Mads Goldschmidt Laursen, Frederik Beyer, Kasper Skov Kjær, Tim Brandt Van Driel, Tobias Harlang, Asmus O. Dohn, Robert J. Hartsock, Silke Nelson, James M. Glownia, Henrik T. Lemke, Morten Christensen, Kelly J. Gaffney, Niels E. Henriksen, Klaus B. Møller, and Martin M. Nielsen.  Ultrafast X-Ray Scattering Measurements of Coherent Structural Dynamics on the Ground-State Potential Energy Surface of a Diplatinum Molecule.  Phys. Rev. Lett., 122(6):63001, 2019.
  * [183] Asmus O. Dohn, Kasper S. Kjær, Tobias B. Harlang, Sophie E. Canton, Martin M. Nielsen, and Klaus B. Møller.  Electron transfer and solvent-mediated electronic localization in molecular photocatalysis.  Inorg. Chem., 55(20):10637–10644, October 2016.
  * [184] Elvar Örn Jónsson, Asmus Ougaard Dohn, and Hannes Jónsson.  Polarizable embedding with a transferable H$\\_2$O potential function I: Formulation and tests on dimer.  J. Chem. Theory Comput., 15(12):6562–6577, 2019.
  * [185] Asmus Ougaard Dohn, Elvar Örn Jónsson, and Hannes Jónsson.  Polarizable embedding with a transferable H$\\_2$O potential function II: Application to H$\\_2$O clusters and liquid water.  J. Chem. Theory Comput., 15(12):6578–6587, 2019.
  * [186] Thomas Olsen.  Designing in-plane heterostructures of quantum spin Hall insulators from first principles: 1T’-MoS$\\_2$ with adsorbates.  Phys. Rev. B, 94:235106, 2016.
  * [187] A.I. Liechtenstein, M.I. Katsnelson, V.P Antropov, and V.A. Gubanov.  Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys.  J. Magn. Magn. Mater., 67:65–74, 1987.
  * [188] N. D. Mermin and H. Wagner.  Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models.  Phys. Rev. Lett., 17:1133–1136, 1966.
  * [189] J. L. Lado and J. Fernández-Rossier.  On the origin of magnetic anisotropy in two dimensional CrI 3.  2D Mater., 4:035002, 2017.
  * [190] Daniele Torelli, Kristian Sommer Thygesen, and Thomas Olsen.  High throughput computational screening for 2D ferromagnetic materials: the critical role of anisotropy and local correlations.  2D Mater., 6:045018, 2019.
  * [191] Daniele Torelli, Hadeel Moustafa, Karsten W. Jacobsen, and Thomas Olsen.  High-throughput computational screening for two-dimensional magnetic materials based on experimental databases of three-dimensional compounds.  npj Comput. Mater., 6:158, 2020.
  * [192] U von Barth and L. Hedin.  A local exchange-correlation potential for the spin polarized case. i.  J. Phys. C Solid State Phys., 5:1629–1642, 1972.
  * [193] Giovanni Scalmani and Michael J. Frisch.  A New Approach to Noncollinear Spin Density Functional Theory beyond the Local Density Approximation.  J. Chem. Theory Comput., 8:2193–2196, 2012.
  * [194] Tae Yun Kim and Cheol-Hwan Park.  Magnetic Anisotropy and Magnetic Ordering of Transition-Metal Phosphorus Trisulfides.  Nano Lett., 21:10114–10121, 2021.
  * [195] D. Ceresoli, U. Gerstmann, A. P. Seitsonen, and F. Mauri.  First-principles theory of orbital magnetization.  Phys. Rev. B, 81(6):060409, 2010.
  * [196] A. J. P. Meyer and G. Asch.  Experimental g’ and g values of Fe, Co, Ni, and their alloys.  J. Appl. Phys., 32(3):S330–S333, 1961.
  * [197] Simon Foner.  High-field antiferromagnetic resonance in Cr$\\_{2}$O$\\_{3}$.  Phys. Rev., 130:183–197, 1963.
  * [198] T. A. Kaplan.  Classical Spin-Configuration Stability in the Presence of Competing Exchange Forces.  Phys. Rev., 116(4):888–889, nov 1959.
  * [199] L. M. Sandratskii.  Energy band structure calculations for crystals with spiral magnetic structure.  Phys. Status Solidi B., 136(1):167–180, 1986.
  * [200] LM Sandratskii.  Insight into the dzyaloshinskii-moriya interaction through first-principles study of chiral magnetic structures.  Phys. Rev. B, 96(2):024450, 2017.
  * [201] Ryogo Kubo.  Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems.  J. Phys. Soc. Jpn., 12(6):570–586, 6 1957.
  * [202] E. K U Gross and Walter Kohn.  Local density-functional theory of frequency-dependent linear response.  Phys. Rev. Lett., 55:2850–2852, 1985.
  * [203] Erich Runge and E. K. U. Gross.  Density-Functional Theory for Time-Dependent Systems.  Phys. Rev. Lett., 52:997–1000, 1984.
  * [204] S. Sharma, J. K. Dewhurst, A. Sanna, and E. K. U. Gross.  Bootstrap approximation for the exchange-correlation kernel of time-dependent density-functional theory.  Phys. Rev. Lett., 107:186401, Oct 2011.
  * [205] AH MacDonald, SH Vosko, and PT Coleridge.  Extensions of the tetrahedron method for evaluating spectral properties of solids.  Journal of Physics C: Solid State Physics, 12(15):2991, 1979.
  * [206] Carlo A. Rozzi, Daniele Varsano, Andrea Marini, Eberhard K. U. Gross, and Angel Rubio.  Exact coulomb cutoff technique for supercell calculations.  Phys. Rev. B, 73:205119, May 2006.
  * [207] Falco Hüser, Thomas Olsen, and Kristian S Thygesen.  How dielectric screening in two-dimensional crystals affects the convergence of excited-state calculations: Monolayer MoS$\\_2$.  Phys. Rev. B, 88(24):245309, 2013.
  * [208] S. Latini, T. Olsen, and K. S. Thygesen.  Excitons in van der waals heterostructures: The important role of dielectric screening.  Phys. Rev. B, 92:245123, Dec 2015.
  * [209] Christopher E. Patrick and Kristian S. Thygesen.  Hubbard-U corrected Hamiltonians for non-self-consistent random-phase approximation total-energy calculations: A study of ZnS, TiO$\\_{2}$, and NiO.  Phys. Rev. B, 93:035133, Jan 2016.
  * [210] Per S. Schmidt and Kristian S. Thygesen.  Benchmark database of transition metal surface and adsorption energies from many-body perturbation theory.  The Journal of Physical Chemistry C, 122(8):4381–4390, 2018.
  * [211] Thomas Olsen, Jun Yan, Jens J. Mortensen, and Kristian S. Thygesen.  Dispersive and covalent interactions between graphene and metal surfaces from the random phase approximation.  Phys. Rev. Lett., 107:156401, Oct 2011.
  * [212] Filipp Furche and Troy Van Voorhis.  Fluctuation-dissipation theorem density-functional theory.  The Journal of Chemical Physics, 122(16), 04 2005.  164106\.
  * [213] Thomas Olsen and Kristian S. Thygesen.  Beyond the random phase approximation: Improved description of short-range correlation by a renormalized adiabatic local density approximation.  Phys. Rev. B, 88:115131, Sep 2013.
  * [214] Thomas Olsen and Kristian S. Thygesen.  Accurate ground-state energies of solids and molecules from time-dependent density-functional theory.  Phys. Rev. Lett., 112:203001, May 2014.
  * [215] Christopher E. Patrick and Kristian S. Thygesen.  Adiabatic-connection fluctuation-dissipation DFT for the structural properties of solids—The renormalized ALDA and electron gas kernels.  The Journal of Chemical Physics, 143(10), 04 2015.  102802\.
  * [216] Thomas Olsen, Christopher E. Patrick, Jefferson E. Bates, Adrienn Ruzsinszky, and Kristian S. Thygesen.  Beyond the RPA and GW methods with adiabatic xc-kernels for accurate ground state and quasiparticle energies.  npj Comput. Mater., 5(1):106, Nov 2019.
  * [217] Paweł Buczek, Arthur Ernst, and Leonid 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.
  * [218] Léon Van Hove.  Time-Dependent Correlations between Spins and Neutron Scattering in Ferromagnetic Crystals.  Phys. Rev., 95(6):1374–1384, 9 1954.
  * [219] Thorbjørn Skovhus, Thomas Olsen, and Henrik M Rønnow.  Influence of static correlation on the magnon dynamics of an itinerant ferromagnet with competing exchange interactions: First-principles study of mnbi.  Phys. Rev. Mater., 6:054402, 5 2022.
  * [220] Thorbjørn Skovhus and Thomas Olsen.  Magnons in antiferromagnetic bcc cr and ${\mathrm{cr}}\\_{2}{\mathrm{o}}\\_{3}$ from time-dependent density functional theory.  Phys. Rev. B, 106:085131, Aug 2022.
  * [221] HN Rojas, Rex William Godby, and RJ Needs.  Space-time method for ab initio calculations of self-energies and dielectric response functions of solids.  Phys. Rev. Lett., 74(10):1827, 1995.
  * [222] Christoph Friedrich, Stefan Blügel, and Arno Schindlmayr.  Efficient implementation of the g w approximation within the all-electron flapw method.  Phys. Rev. B, 81(12):125102, 2010.
  * [223] Ivan Duchemin and Xavier Blase.  Robust analytic-continuation approach to many-body gw calculations.  J. Chem. Theory Comput., 16(3):1742–1756, 2020.
  * [224] Dario A Leon, Claudia Cardoso, Tommaso Chiarotti, Daniele Varsano, Elisa Molinari, and Andrea Ferretti.  Frequency dependence in g w made simple using a multipole approximation.  Phys. Rev. B, 104(11):115157, 2021.
  * [225] Filip A Rasmussen, Per S Schmidt, Kirsten T Winther, and Kristian S Thygesen.  Efficient many-body calculations for two-dimensional materials using exact limits for the screened potential: Band gaps of MoS$\\_2$, H-BN, and phosphorene.  Phys. Rev. B, 94(15):155406, 2016.
  * [226] Andrea Marini, Conor Hogan, Myrta Grüning, and Daniele Varsano.  Yambo: an ab initio tool for excited state calculations.  Comput. Phys. Commun., 180(8):1392–1403, 2009.
  * [227] Thomas Olsen, Simone Latini, Filip Rasmussen, and Kristian S. Thygesen.  Simple Screened Hydrogen Model of Excitons in Two-Dimensional Materials.  Phys. Rev. Lett., 116:056401, 2016.
  * [228] Thomas Olsen.  Unified Treatment of Magnons and Excitons in Monolayer CrI$\\_3$ from Many-Body Perturbation Theory.  Phys. Rev. Lett., 127:166402, 2021.
  * [229] Zhen Li, Patrizio Graziosi, and Neophytos Neophytou.  Deformation potential extraction and computationally efficient mobility calculations in silicon from first principles.  Phys. Rev. B, 104(19):195201, 2021.
  * [230] Amy Y Liu and Andrew A Quong.  Linear-response calculation of electron-phonon coupling parameters.  Phys. Rev. B, 53(12):R7575, 1996.
  * [231] Feliciano Giustino.  Electron-phonon interactions from first principles.  Rev. Mod. Phys., 89:015003, Feb 2017.
  * [232] Michael Walter and Michael Moseler.  Ab Initio Wavelength-Dependent Raman Spectra: Placzek Approximation and Beyond.  J. Chem. Theory Comput., 16(1):576–586, January 2020.
  * [233] Ado Jorio, Mildred S Dresselhaus, Riichiro Saito, and Gene Dresselhaus.  Raman spectroscopy in graphene related systems.  John Wiley & Sons, 2011.
  * [234] Atsushi Togo, Laurent Chaput, Terumasa Tadano, and Isao Tanaka.  Implementation strategies in phonopy and phono3py.  J. Phys. Condens. Matter, 35(35):353001, 2023.
  * [235] Robert W. Boyd.  Nonlinear Optics.  Elsevier Science Publishing Co Inc., Amsterdam, 3rd edition, 2008.
  * [236] Alireza Taghizadeh, F. Hipolito, and T. G. Pedersen.  Linear and nonlinear optical response of crystals using length and velocity gauges: Effect of basis truncation.  Phys. Rev. B, 96:195413, 2017.
  * [237] Claudio Aversa and J. E. Sipe.  Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis.  Phys. Rev. B, 52:14636–14645, 1995.
  * [238] Arto Sakko, Tuomas P Rossi, and Risto M Nieminen.  Dynamical coupling of plasmons and molecular excitations by hybrid quantum/classical calculations: time-domain approach.  Journal of Physics: Condensed Matter, 26(31):315013, jul 2014.
  * [239] Tuomas P. Rossi, Kirsten T. Winther, Karsten W. Jacobsen, Risto M. Nieminen, Martti J. Puska, and Kristian S. Thygesen.  Effect of edge plasmons on the optical properties of MoS$\\_{2}$ monolayer flakes.  Phys. Rev. B, 96:155407, Oct 2017.
  * [240] Priyank V. Kumar, Tuomas P. Rossi, Mikael Kuisma, Paul Erhart, and David J. Norris.  Direct hot-carrier transfer in plasmonic catalysis.  Faraday Discuss., 214:189–197, 2019.
  * [241] Priyank V. Kumar, Tuomas P. Rossi, Daniel Marti-Dafcik, Daniel Reichmuth, Mikael Kuisma, Paul Erhart, Martti J. Puska, and David J. Norris.  Plasmon-Induced Direct Hot-Carrier Transfer at Metal–Acceptor Interfaces.  Acs Nano, 13(3):3188–3195, 2019.
  * [242] Tuomas P. Rossi, Timur Shegai, Paul Erhart, and Tomasz J. Antosiewicz.  Strong plasmon-molecule coupling at the nanoscale revealed by first-principles modeling.  Nat. Commun., 10(1):3336, 2019.
  * [243] Tuomas P. Rossi, Paul Erhart, and Mikael Kuisma.  Hot-Carrier Generation in Plasmonic Nanoparticles: The Importance of Atomic Structure.  Acs Nano, 14(8):9963–9971, 2020.
  * [244] Jakub Fojt, Tuomas P. Rossi, Mikael Kuisma, and Paul Erhart.  Hot-Carrier Transfer across a Nanoparticle–Molecule Junction: The Importance of Orbital Hybridization and Level Alignment.  Nano Lett., 22(21):8786–8792, 2022.
  * [245] Daniel Sorvisto, Patrick Rinke, and Tuomas P. Rossi.  Single-atom dopants in plasmonic nanocatalysts.  The Journal of Physical Chemistry C, 127(18):8585–8590, 2023.
  * [246] K. Yabana and G. F. Bertsch.  Time-dependent local-density approximation in real time.  Phys. Rev. B, 54(7):4484–4487, 1996.
  * [247] Sami Malola, Lauri Lehtovaara, Jussi Enkovaara, and Hannu Häkkinen.  Birth of the localized surface plasmon resonance in monolayer-protected gold nanoclusters.  Acs Nano, 7(11):10263–10270, 2013.
  * [248] Tuomas Rossi, Paul Erhart, and Mikael Kuisma.  Code for "hot-carrier generation in plasmonic nanoparticles: The importance of atomic structure".  Zenodo.