# GPAW: open Python package for electronic-structure calculations

Jens Jørgen Mortensen  CAMd, Department of Physics, Technical University of
Denmark, 2800 Kgs. Lyngby, Denmark  Ask Hjorth Larsen  CAMd, Department of
Physics, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Mikael
Kuisma  CAMd, Department of Physics, Technical University of Denmark, 2800
Kgs. Lyngby, Denmark  Aleksei V. Ivanov  Riverlane Ltd, St Andrews House, 59
St Andrews Street, Cambridge, CB2 3BZ, United Kingdom  Alireza Taghizadeh
CAMd, Department of Physics, Technical University of Denmark, 2800 Kgs.
Lyngby, Denmark  Andrew Peterson  School of Engineering, Brown University,
Providence, Rhode Island, USA  Anubhab Haldar  Department of Electrical and
Computer Engineering, Boston University, United States of America.  Asmus
Ougaard Dohn  Department of Physics, Technical University of Denmark, 2800
Lyngby, Denmark; Science Institute and Faculty of Physical Sciences, VR-III,
University of Iceland, Reykjavík 107, Iceland  Christian Schäfer  Department
of Physics, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden
Elvar Örn Jónsson  Science Institute and Faculty of Physical Sciences,
University of Iceland, VR-III, 107 Reykjavík, Iceland  Eric D. Hermes
Quantum-Si, 29 Business Park Drive, Branford, CT 06405, USA  Fredrik Andreas
Nilsson  CAMd, Department of Physics, Technical University of Denmark, 2800
Kgs. Lyngby, Denmark  Georg Kastlunger  CatTheory, Department of Physics,
Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Gianluca Levi
Science Institute and Faculty of Physical Sciences, University of Iceland, VR-
III, 107 Reykjavík, Iceland  Hannes Jónsson  Science Institute and Faculty of
Physical Sciences, University of Iceland, VR-III, 107 Reykjavík, Iceland
Hannu Häkkinen  Departments of Physics and Chemistry, Nanoscience Center,
University of Jyväskylä, FI-40014 Jyväskylä, Finland  Jakub Fojt  Department
of Physics, Chalmers University of Technology, SE-412 96 Gothenburg, Sweden
Jiban Kangsabanik  CAMd, Department of Physics, Technical University of
Denmark, 2800 Kgs. Lyngby, Denmark  Joachim Sødequist  CAMd, Department of
Physics, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Jouko
Lehtomäki  Department of Applied Physics, Aalto University, P.O. Box 11100,
00076 Aalto, Finland  Julian Heske  CAMd, Department of Physics, Technical
University of Denmark, 2800 Kgs. Lyngby, Denmark  Jussi Enkovaara  CSC – IT
Center for Science Ltd., P.O. Box 405, FI-02101 Espoo, Finland  Kirsten
Trøstrup Winther  SUNCAT Center for Interface Science and Catalysis, SLAC
National Accelerator Laboratory, Menlo Park, California 94025, United States
Marcin Dulak  Marko M. Melander  Department of Chemistry, Nanoscience Center,
University of Jyväskylä, FI-40014 Jyväskylä, Finland  Martin Ovesen  CAMd,
Department of Physics, Technical University of Denmark, 2800 Kgs. Lyngby,
Denmark  Martti Louhivuori  CSC – IT Center for Science Ltd., P.O. Box 405,
FI-02101 Espoo, Finland  Michael Walter  University of Freiburg, FIT Freiburg
Centre for Interactive Materials and Bioinspired Technologies, Georges-Köhler-
Allee 105, 79110 Freiburg, Germany  Morten Gjerding  CAMd, Department of
Physics, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Olga
Lopez-Acevedo  Biophysics of Tropical Diseases, Max Planck Tandem Group,
University of Antioquia UdeA, 050010 Medellin, Colombia  Paul Erhart
Department of Physics, Chalmers University of Technology, SE-412 96
Gothenburg, Sweden  Robert Warmbier  School of Physics and Mandelstam
Institute for Theoretical Physics, University of the Witwatersrand, 1 Jan
Smuts Avenue, 2001, Johannesburg, South Africa  Rolf Würdemann  Freiburger
Materialforschungszentrum, Universität Freiburg, Stefan-Meier-Straße 21,
D-79104 Freiburg, Germany  Sami Kaappa  Computational Physics Laboratory,
Tampere University, P.O. Box 692, FI-33014 Tampere, Finland  Simone Latini
Nanomade, Department of Physics, Technical University of Denmark, 2800 Kgs.
Lyngby, Denmark  Tara Maria Boland  CAMd, Department of Physics, Technical
University of Denmark, 2800 Kgs. Lyngby, Denmark  Thomas Bligaard  Department
of Energy Conversion and Storage, Technical University of Denmark, DK-2800
Lyngby  Thorbjørn Skovhus  CAMd, Department of Physics, Technical University
of Denmark, 2800 Kgs. Lyngby, Denmark  Toma Susi  University of Vienna,
Faculty of Physics, Boltzmanngasse 5, 1090 Vienna, Austria  Tristan Maxson
Department of Chemical and Biological Engineering, The University of Alabama,
Tuscaloosa, AL 35487, USA  Tuomas Rossi  CSC – IT Center for Science Ltd.,
P.O. Box 405, FI-02101 Espoo, Finland  Xi Chen  School of Physical Science and
Technology, Lanzhou University, Lanzhou, Gansu 730000, China  Yorick Leonard
A. Schmerwitz  Science Institute and Faculty of Physical Sciences, University
of Iceland, VR-III, 107 Reykjavík, Iceland  Jakob Schiøtz  CAMd, Department of
Physics, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Thomas
Olsen  CAMd, Department of Physics, Technical University of Denmark, 2800 Kgs.
Lyngby, Denmark  Karsten Wedel Jacobsen  CAMd, Department of Physics,
Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Kristian Sommer
Thygesen  CAMd, Department of Physics, Technical University of Denmark, 2800
Kgs. Lyngby, Denmark

###### Abstract

We review the GPAW open-source Python package for electronic structure
calculations. GPAW is based on the projector-augmented wave method and can
solve the self-consistent density functional theory (DFT) equations using
three different wave-function representations, namely real-space grids, plane
waves, and numerical atomic orbitals. The three representations are
complementary and mutually independent and can be connected by transformations
via the real-space grid. This multi-basis feature renders GPAW highly
versatile and unique among similar codes. By virtue of its modular structure,
the GPAW code constitutes an ideal platform for implementation of new features
and methodologies. Moreover, it is well integrated with the Atomic Simulation
Environment (ASE) providing a flexible and dynamic user interface. In addition
to ground-state DFT calculations, GPAW supports many-body GW band structures,
optical excitations from the Bethe–Salpeter Equation (BSE), variational
calculations of excited states in molecules and solids via direct
optimization, and real-time propagation of the Kohn–Sham equations within
time-dependent DFT. A range of more advanced methods to describe magnetic
excitations and non-collinear magnetism in solids are also now available. In
addition, GPAW can calculate non-linear optical tensors of solids, charged
crystal point defects, and much more. Recently, support of GPU acceleration
has been achieved with minor modifications of the GPAW code thanks to the CuPy
library. We end the review with an outlook describing some future plans for
GPAW.

###### Contents

  1. I Introduction
  2. II Why GPAW?
    1. II.1 User’s perspective
    2. II.2 Developer’s perspective
  3. III Ground-state DFT
    1. III.1 Projector augmented-wave method
    2. III.2 Numerical implementation
      1. III.2.1 Wave-function representations
      2. III.2.2 All-electron quantities
      3. III.2.3 Solving the Kohn–Sham equation
      4. III.2.4 Updating wavefunctions in dynamics
      5. III.2.5 Direct minimization
      6. III.2.6 Convergence criteria
      7. III.2.7 PAW data sets and pseudopotentials
      8. III.2.8 Parallelization
      9. III.2.9 GPU implementation
    3. III.3 XC-functionals
      1. III.3.1 Libxc and Libvdwxc
      2. III.3.2 GLLB-sc
      3. III.3.3 Hubbard U
      4. III.3.4 Hybrids
      5. III.3.5 SIC
      6. III.3.6 BEEF
  4. IV Ion dynamics
    1. IV.1 Structure relaxation
    2. IV.2 Reaction paths and barriers
    3. IV.3 Global structure optimization
    4. IV.4 Molecular dynamics and QM/MM
  5. V Magnetism and spin
    1. V.1 Spin–orbit coupling
    2. V.2 Self-consistent non-collinear magnetism
    3. V.3 Orbital magnetization
    4. V.4 Constant B-field
    5. V.5 Spin spirals
  6. VI Response functions and excitations
    1. VI.1 Linear-response TDDFT
      1. VI.1.1 Implementation for periodic systems
      2. VI.1.2 Spectral representation
    2. VI.2 Dielectric function
      1. VI.2.1 Screening in low-dimensional systems
    3. VI.3 Adiabatic-connection fluctuation-dissipation theorem
    4. VI.4 Magnetic response
      1. VI.4.1 Transverse magnetic susceptibility
      2. VI.4.2 The spectrum of transverse magnetic excitations
      3. VI.4.3 Liechtenstein MFT
    5. VI.5 GW approximation
    6. VI.6 Bethe–Salpeter Equation (BSE)
    7. VI.7 Electron–phonon coupling
    8. VI.8 Raman spectrum
    9. VI.9 Quadratic optical response functions
  7. VII Real-time TDDFT
    1. VII.1 Kohn–Sham decomposition
    2. VII.2 Hot-carrier analysis
    3. VII.3 Circular dichroism for molecules
    4. VII.4 Radiation reaction potential
    5. VII.5 Ehrenfest dynamics
  8. VIII Excited-state DFT methods
    1. VIII.1 Improved virtual orbitals
    2. VIII.2 Variational excited-state calculations
      1. VIII.2.1 Direct orbital optimization
      2. VIII.2.2 Example applications of direct optimization
  9. IX Other features
    1. IX.1 Electric polarization
    2. IX.2 Berry phases and band topology
    3. IX.3 Wannier functions
    4. IX.4 Point defect calculations with hybrid functionals
    5. IX.5 Point-group symmetry representations
    6. IX.6 Band-structure unfolding
    7. IX.7 The QEH model
    8. IX.8 Solvent models
    9. IX.9 Charged electrochemical interfaces
    10. IX.10 Constrained DFT
    11. IX.11 Orbital-free DFT
    12. IX.12 Zero-field splitting
    13. IX.13 Hyperfine coupling
  10. X Outlook
  11. XI Acknowledgements

##  I Introduction

The electronic-structure (ES) problem, i.e. the solution to the time-
independent Schrödinger equation for a collection of electrons and atomic
nuclei, forms the starting point for the quantum-mechanical treatment of
matter. Indeed, all chemical and physical properties of any substance (solid,
molecule, surface, etc.) can in principle be obtained from the energies and
wave functions that constitute the solution. A pioneering step towards solving
the many-body ES problem was the formulation and formal proof of density
functional theory (DFT) by Hohenberg and Kohn in 1964 [1], and a practical
scheme for its solution by Kohn and Sham in 1965 [2]. Today, most codes
solving the ES problem from first principles are based on DFT. Such codes are
extremely powerful and allow one to determine the atomic structure of solids
and molecules containing hundreds of atoms with a relative error below 1% [3,
4, 5]. Once the atomic structure of the compound has been solved, its
properties (electronic, magnetic, optical, topological, etc.) can in principle
be determined. The evaluation of properties often involves theories beyond the
formal DFT framework to account for effects such as temperature and lattice
vibrations [6, 7], many-body interactions in excited states [8, 9], or time
dependence [10, 11]. As such, first-principles atomistic calculations often
involve two successive phases: the solution of the ground-state ES problem
(including ion dynamics) and the subsequent evaluation of physical properties.
This review is structured accordingly as Secs. III-IV deal with the first
phase while Secs. V-IX are devoted to the second.

In recent years, the scientific significance of ES codes has shifted from a
useful tool to describe and understand matter at the atomic scale to an
independent driver of the discovery and development of new materials [12, 13,
14, 15, 16]. This change in scope has been fueled by the exponential increase
in computer power accompanied by improved numerical algorithms [17, 18] as
well as the use of workflow management software for high-throughput
computations [19, 20, 21, 22] and the adoption of machine-learning techniques
to leverage the rapidly growing data generated by ES codes [23, 24, 25]. In
parallel with these capacity-extending developments, continuous progress in
the fundamental description of exchange–correlation effects has advanced the
predictive power of ES calculations to a level where they rival experiments in
terms of accuracy for many important properties [26, 27, 28, 29, 30, 31].

The GPAW code was originally intended as a Python-based multigrid solver of
the basic DFT equations within the projector-augmented wave (PAW) formalism
[32]. The name GPAW accordingly was an abbreviation for “grid-based projector-
augmented waves”. Today, other choices than regular grids for representations
of the wave functions exist in GPAW, but the name has stuck. During the years
2005–2010, GPAW evolved to a full-blown DFT package [33] supporting most of
the functionality expected from a modern ES code, in addition to a few more
specialised features including real-time propagation of wave functions [34]
and an alternative basis of numerical atomic orbitals (referred to as the LCAO
basis) [35] to supplement the real-space grid. In 2011, a plane-wave (PW)
basis set was also implemented. Today, the possibility to use three different
types of basis sets and even combining them within a single run remains a
unique feature of GPAW, rendering the code very versatile.

The implementation of the PW basis set laid the groundwork for GPAW’s linear-
response module, which today supports the calculation of linear response
functions [36], total energies from the adiabatic connection fluctuation-
dissipation theorem [29, 37], the GW self-energy method for quasiparticle band
structures [38], the Bethe–Salpeter Equation (BSE) for optical excitations
[39], and more. The code also supports a wide range of features related to the
calculation of magnetism and spin–orbit effects. Examples include spin-spiral
calculations using the generalized Bloch theorem [40], external magnetic
fields, orbital magnetization, magnetic anisotropy [41], adiabatic magnon
dispersions from the magnetic force theorem [42], and dynamic magnetic
response from TDDFT [43]. For solids, the k-space Berry phases can be computed
directly from the Bloch orbitals and may be used to obtain the spontaneous
polarization [44], Born effective charges, piezoelectric response tensors [45]
and various indices characterising the band topology [46].

In addition, GPAW can compute the localisation matrices forming the basis for
the construction of Wannier functions with e.g. the Atomic Simulation
Environment (ASE) [47] or Wannier90 [48]. Electrostatic corrections to the
formation energies of charged point defects in insulators are implemented as
are calculations of the hyperfine coupling and zero-field splitting for
localised electron spins. GPAW also offers the possibility to perform time-
independent, variational calculations of localised electronic excitations, in
e.g. molecules or at crystal point defects, using direct orbital optimisation
strategies implemented for all three types of basis sets [49, 50, 51]. This
provides an efficient and robust alternative to traditional “$\Delta$SCF”
approaches. GPAW can also be used to describe ultrafast electron dynamics
within time-dependent density functional theory (TDDFT) with wave functions
represented either on a real space grid [34] or in the LCAO basis [52]. The
latter can provide a significant speed-up due to the relatively small size of
the basis [53, 54, 55]. The LCAO representation also forms the basis for
calculation of electron–phonon couplings as well as non-linear optical spectra
such as Raman scattering [56] (which can alternatively be obtained in the PW
mode as a finite difference of the dielectric tensor), second-harmonics
generation [57], and shift currents [58] using higher-order perturbation
theory.

##  II Why GPAW?

###  II.1 User’s perspective

There are dozens of electronic-structure codes available for the interested
user. The codes differ in their license (in particular, whether open or
proprietary), the underlying programming language (e.g. Fortran, C, Python),
their treatment of core electrons (all-electron versus pseudopotentials), the
employed representations of the wave functions (plane waves, atom-centered
orbitals, real-space grids), and the beyond-DFT features they support. Why
should one choose GPAW?

In this section, we describe some of the features that make GPAW interesting
from the point of view of a common user who wants to perform electronic-
structure calculations. The next section focuses on its possibilities for more
advanced users, who perhaps want to modify the code or implement completely
new functionalities.

A first point to note is that GPAW is written almost exclusively in Python and
is directly integrated with the Atomic Simulation Environment. This
integration with ASE makes the setup, control, and analysis of calculations
easy and flexible. The programming language is of course a key issue for
developers, but also the common user benefits from Python and the ASE/GPAW
integration. A typical stand-alone program only offers a fixed (though of
course possibly large) set of tasks that it can perform, while Python
scripting allows for a more flexible use of the code. This could for example
mean combining several different GPAW calculations in new ways. Another
advantage is that “inner parts” of the code like the density or the Kohn–Sham
eigenvalues are directly accessible in a structured format within Python for
further analysis. It is even possible to “open up the main loop” of GPAW and
have access to, inspect, and also modify key quantities during program
execution (see Fig. 1).

As already mentioned in the introduction, GPAW distinguishes itself from other
available ES codes by supporting three different ways of representing the wave
functions. The most commonly used basis set is plane waves (PW), which is
appropriate for small or medium-size systems, where high precision is
required. Convergence is easily and systematically controlled by tuning the
cut-off energy. A large number of advanced features and “beyond-DFT” methods
are available in the PW mode. These include the calculation of hybrid
functionals, RPA total energies, linear-response TDDFT, and many-body
perturbation theory techniques like GW and the Bethe–Salpeter equations. The
new GPU implementation also uses the PW mode.

The wave functions can alternatively be represented on real-space grids, which
was the original approach in GPAW. The implementation of this so-called
finite-difference (FD) mode relies on multi-grid solutions of the Poisson and
Kohn–Sham equations. The FD mode allows for more flexible boundary conditions
than the PW mode, which is restricted to periodic supercells. The boundary
conditions may for example be taken to reflect the charge distribution in the
unit cell. Calculations in the FD mode can be systematically converged through
lowering of the grid spacing, but the approach to full convergence is slower
than in the PW mode. The FD mode is particularly well suited for large systems
because the wave-function representation allows for large-scale
parallelization through real-space decomposition. Furthermore, it is possible
to perform time-propagation TDDFT including Ehrenfest dynamics in this mode.

The third representation of the wave functions is a basis of numerical atom-
centered orbitals in the linear combination of atomic orbitals (LCAO) mode.
The size of the basis set can be varied through inclusion of more angular
momentum channels, additional orbitals within a channel, or polarization
functions. GPAW comes with a standard set of orbitals, but a basis-set
generator is included with the code so that users may construct different
basis sets depending on their needs and requirements. The LCAO mode is
generally less accurate than the PW and FD modes, but it allows for the
treatment of considerably larger systems – more than ten thousand atoms. It is
also possible to study electron dynamics through a fast implementation of
time-propagation DFT, and Ehrenfest dynamics is under development.

As explained, the different modes have different virtues and limitations, and
it can therefore be an advantage to apply several modes in a project. For
larger systems, it is for example possible to divide a structure optimization
into two steps. First, an optimization is performed with the fast LCAO basis
leading to an approximately correct structure. This is then followed by an
optimization in either the PW or FD mode, which now requires much fewer steps
because of the good initial configuration. Due to the ASE/Python interface
this combined calculation can easily be performed within a single script.

Since GPAW was originally created with the FD mode only, and the LCAO mode was
added next, some features have been implemented for only those modes. Examples
are real-time TDDFT (see section VII) and electron–phonon coupling (see
section VI.7). Conversely, some new features only work for the PW mode, which
was added after the real-space modes. Examples are RPA total energies (see
section VI.3) and calculation of the stress tensor. To summarize, given the
limitations just mentioned, users should most of the time use PW or LCAO mode
and the choice will depend on the accuracy needed and the resources available.

###  II.2 Developer’s perspective

The GPAW source code is written in the Python and C languages and is hosted on
GitLab [59] licensed under the GNU General Public License v3.0. This ensures
the transparency of all features and allows developers to fully customise
their experience and contribute new features to the community.

An advantage of having a Python code is that the Python script you write to
carry out your calculations will have access to everything inside a GPAW
calculation. An example showing the power and flexibility this affords is the
possibility to have user-code inserted inside the self-consistent field (SCF)
loop as demonstrated in Fig. 1.

⬇

import os

import psutil

process = psutil.Process(os.getpid())

atoms = …

calc = …

for ctx in calc.icalculate(atoms):

# Inside SCF loop

mem = process.memory_info().rss

# Write memory usage to log-file:

ctx.log(f’MEM:␣{ctx.niter}␣{mem}’)

if ctx.niter == 15:

# Stop after 15 iterations

break

Figure 1:  The variable calc is the ground-state DFT calculator object and its
icalculate method yields a context object at every self-consistent field (SCF)
step. As seen, one can use this in a for-loop to implement special logic for
termination of the SCF iterations or for diagnostics. In this example, the
memory usage is written to the log-file for the first 15 SCF iterations.

At the time of this writing (July 2023), GPAW has two versions of the ground-
state DFT code in the main branch of the code. There is the older version that
has grown organically since the birth of GPAW: it has many features, but also
a lot of technical debt that makes it harder to maintain and less ideal to
build new features on top of. The newer ground-state code addresses these
issues by having a better overall design.

The new design greatly improves the ease of implementation of new features.
The goal is to have the new code be feature-complete so that it can pass the
complete test suite and then delete the old code once that is achieved. We
recommend that new work on GPAW is done on top of the new code even though
certain features are not yet production-ready. Three new features, not present
in the old code base, have already been implemented based on the new code: GPU
implementation of PW mode calculations (see section III.2.9), reuse of the
wave functions after unit-cell changes during cell optimization, and spin-
spiral calculations (see section V.5).

GPAW uses pytest [60] for its test suite that currently consists of
approximately 1600 unit and integration tests (see Table 1). A subset of those
tests runs as part of GitLab’s continuous integration (CI) thereby checking
the correctness of every code change. Unfortunately, the full test suite is
too time-consuming to run as part of CI, so we run that nightly both in serial
as well as in parallel using MPI.

Many of the code examples in GPAW’s documentation, exercises and tutorials
[61] require resources (time and number of CPUs) beyond what would make sense
to run as part of the pytest test suite. For that, we use MyQueue [21] to
submit those scripts as jobs to a local supercomputer every weekend. At the
moment this amounts to approximately 5200 core-hours of calculations.

As can be seen from Table 1, the majority of the code is written in Python,
which is an interpreted language that is easy to read, write and debug.

Table 1: Number of files and number of lines of code in the git repository of GPAW. The Python source-code files are split into three parts: the actual code, the test suite, and code examples in the documentation.  Type of file | Files | Lines
---|---|---
Python (the code) | 513 | 146,604
C | 80 | 19,719
Python (test suite) | 681 | 47,147
Python (documentation) | 744 | 32,014

Interpreter-executed code will not run as efficiently as code that is compiled
to native machine code. It is therefore important to make sure that the places
in the code where most of the time is spent (hot spots) are in native machine
code and not in the interpreter. GPAW achieves this by implementing the hot
spots in C-code with Python wrappers that can be called from the Python code.
Examples of such computationally intensive tasks are applying a finite-
difference stencil to a uniform grid, interpolating from one uniform grid to
another, or calculation of overlaps between projector functions and wave
functions. In addition, we have Python interfaces to the numerical libraries
FFTW [62], ScaLAPACK [63], ELPA [17], BLAS, Libxc [64, 65], libvdwxc [66], and
MPI. Finally, GPAW makes heavy use of the NumPy [67] and Scipy [68] Python
packages. NumPy provides us with the numpy.ndarray data type which is an
$N$-dimensional array that we use for storing wave functions, electron
densities, potentials, matrices like the overlap matrix or the LCAO wave
function coefficients, and much more. The use of NumPy arrays allows us to use
the many sub-modules of SciPy to manipulate data. This also gives us an
efficient memory layout allowing us to simply pass a pointer to the memory
whenever we need to call the C-code from the Python code. With this strategy,
we can get away with having most of the code written in a relatively slow
interpreted language and still have most of the time spent in highly efficient
C-code or optimized numerical libraries.

The advantage of the original FD-mode, where there are no Fourier transforms
of the wave functions, is that the algorithms should parallelize well for
large systems. In practice, it has turned out that FD-mode has a number of
disadvantages: 1) Due to integrals over the unit cell being done as sums over
grid-points, there will be a small periodic energy variation as you translate
atoms and the period of the variation will be equal to the grid-spacing used
(the so-called egg-box error); 2) The system sizes that are typically most
interesting for applications of DFT are too small for the parallel scalability
to be the decisive advantage; 3) The memory used to store the wave functions
on uniform grids in real-space is significant. In contrast, the PW mode has
practically no egg-box error, is very efficient for the most typical system
sizes and often uses a factor of 10 less memory compared to an FD mode
calculation of similar accuracy. The main advantages for LCAO mode are low
memory usage and high efficiency for large systems; for small unit cells with
many $\mathbf{k}$-points the PW mode is most efficient. The disadvantages of
LCAO-mode are 1) egg-box errors and 2) the difficulty of reaching the complete
basis-set limit. For PW and FD modes, the complete basis set limit is easy to
reach by simply increasing the number of plane-waves or grid-points
respectively.

##  III Ground-state DFT

###  III.1 Projector augmented-wave method

The diverging Coulomb potential causes rapid oscillations in electronic wave
functions near the nuclei, and special care is required to be able to work
with smooth wave functions. The projector augmented-wave (PAW) method by
Blöchl [69] is a widely adopted generalization of pseudopotential methods,
utilizing their strength of smooth pseudo wave functions while retaining a
mapping from all-electron wave functions ($|\psi_{n}\rangle$) to pseudo wave
functions ($|\tilde{\psi}_{n}\rangle$).

The crux is to define a linear transformation $\hat{T}$ from pseudo to all-
electron space,

$\hat{T}|\tilde{\psi}\rangle=|\psi\rangle,$ (1)

where

$\hat{T}=1+\sum_{a}\sum_{i}\left(|\phi^{a}_{i}\rangle-|\tilde{\phi}^{a}_{i}\rangle\right)\langle\tilde{p}^{a}_{i}|.$
(2)

Here $\tilde{p}^{a}_{i}(\mathbf{r})$, $\phi^{a}_{i}(\mathbf{r})$ and
$\tilde{\phi}^{a}_{i}(\mathbf{r})$ are called projectors, partial waves and
pseudo-partial waves, respectively. The pseudo-partial waves and projectors
are chosen to be biorthogonal, $\int{\rm
d}{\mathbf{r}}\,\tilde{\phi}^{a}_{i}(\mathbf{r})\tilde{p}^{a}_{j}(\mathbf{r})=\delta_{ij}$,
allowing for an approximate closure relation

$\sum_{i}\tilde{\phi}^{a}_{i}(\mathbf{r})\tilde{p}^{a}_{i}(\mathbf{r}^{\prime})\approx\delta(\mathbf{r}-\mathbf{r}^{\prime}),$
(3)

which is utilized heavily to obtain an efficient, but all-electron, picture.
In addition to biorthogonality, the pseudo and all-electron partial waves are
chosen to be equal outside the PAW augmentation sphere cutoff radius
$\tilde{\phi}_{i}(\mathbf{r})=\phi_{i}(\mathbf{r})$, $r>r_{c}.$

The basic recipe for converting an operator $\hat{O}$ is
$\langle\psi_{n}|\hat{O}\psi_{n}\rangle=\langle\tilde{\psi}_{n}|\underbrace{\hat{T}^{\dagger}\hat{O}\hat{T}}_{\hat{\tilde{O}}}\tilde{\psi}_{n}\rangle$.
For example, the all-electron Kohn–Sham equations

$\hat{H}[n({\mathbf{r}})]|\psi_{n}\rangle=\epsilon_{n}|\psi_{n}\rangle,$ (4)

where $\hat{H}$ is the single-particle all-electron Kohn–Sham Hamiltonian
operator, can be transformed to their PAW counterparts:

$\hat{T}^{\dagger}\hat{H}[n({\mathbf{r}})]\hat{T}|\tilde{\psi}_{n}\rangle=\epsilon_{n}\hat{T}^{\dagger}\hat{T}|\tilde{\psi}_{n}\rangle.$
(5)

We have used Eq. 1, and also multiplied with $\hat{T}^{\dagger}$ from left to
make its dual space the pseudo one, i.e.
$\langle\tilde{\psi}|\in\mathcal{H}^{*}$ can act from left. This results in a
PAW Hamiltonian and in PAW overlap operators as follows:

$\displaystyle\hat{\tilde{H}}$
$\displaystyle=\hat{T}^{\dagger}\hat{H}[n({\mathbf{r}})]\hat{T}$
$\displaystyle=-\frac{1}{2}\nabla^{2}+\tilde{v}_{\rm
KS}(\mathbf{r})+\sum_{a}\sum_{ii^{\prime}}|\tilde{p}^{a}_{i}\rangle\Delta
H^{a}_{ii^{\prime}}\langle\tilde{p}^{a}_{i^{\prime}}|,$ (6)
$\displaystyle\hat{\tilde{S}}$
$\displaystyle=\hat{T}^{\dagger}\hat{T}=1+\sum_{a}\sum_{ii^{\prime}}|\tilde{p}^{a}_{i}\rangle\Delta
S^{a}_{ii^{\prime}}\langle\tilde{p}^{a}_{i^{\prime}}|.$ (7)

Terms such as $\Delta H^{a}_{ii^{\prime}}$ and $\Delta S^{a}_{ii^{\prime}}$
represent so-called PAW-corrections. In each part of the description, which
handles a particular kind of operators, such as kinetic energy, spin-
operators, or the electrostatic potential, respective PAW-corrections must be
calculated. The most crucial ones are precalculated, such as overlap, kinetic
energy, Coulomb, and stored in the ‘setup’-file, which also stores the partial
waves and projectors. As an example, the overlap PAW corrections are
precalculated to setup as follows:

$\Delta S^{a}_{ii}=\int{\rm
d}{\mathbf{r}}\,\phi_{i}^{a}(\mathbf{r})\phi_{i^{\prime}}^{a}(\mathbf{r})-\int{\rm
d}{\mathbf{r}}\tilde{\phi}_{i}^{a}(\mathbf{r})\tilde{\phi}_{i^{\prime}}^{a}(\mathbf{r}).$
(8)

We further define the atomic density matrices as

$D^{a}_{\sigma ii^{\prime}}=\sum_{n}f_{n}\langle\tilde{\psi}_{\sigma
n}|\tilde{p}^{a}_{i}\rangle\langle\tilde{p}^{a}_{i^{\prime}}|\tilde{\psi}_{\sigma
n}\rangle.$ (9)

The atomic density matrix contains all information required to construct PAW
corrections to any local all-electron expectation value:

$\langle
O\rangle=\sum_{n}\langle\tilde{\psi}_{n}|\tilde{O}|\tilde{\psi}_{n}\rangle+\sum_{a}\sum_{ii^{\prime}}D^{a}_{ii^{\prime}}\Delta
O^{a}_{i^{\prime}i}.$ (10)

The all-electron atomic density can be constructed as

$n^{a}_{\sigma}(\mathbf{r})=\sum_{ii^{\prime}}D^{a}_{\sigma
ii^{\prime}}\phi^{a*}_{i}(\mathbf{r})\phi^{a}_{i^{\prime}}(\mathbf{r})+n_{\text{core}}^{a},$
(11)

and the corresponding equation for pseudo-densities holds with
$n\rightarrow\tilde{n}$ and $\phi\rightarrow\tilde{\phi}$.

Since the exchange–correlation (xc) potential is non-linear, the PAW
corrections must be evaluated explicitly. The xc PAW corrections are performed
by constructing the atomic all-electron and pseudo-electron densities as given
by Eq. 11:

$\displaystyle\Delta H^{a}_{{\rm xc},ii^{\prime}}$ $\displaystyle=$
$\displaystyle\int{\rm d}{\mathbf{r}}\phi_{i}^{a}(\mathbf{r})v_{\rm
xc}[n^{a}](\mathbf{r})\phi_{i^{\prime}}^{a}(\mathbf{r})$ (12)
$\displaystyle-\int{\rm d}{\mathbf{r}}\tilde{\phi}_{i}^{a}(\mathbf{r})v_{\rm
xc}[\tilde{n}^{a}](\mathbf{r})\tilde{\phi}_{i^{\prime}}^{a}(\mathbf{r}).$

This integral is numerically evaluated in a Cartesian product of a Lebedev
angular grid and a non-uniform radial grid (denser mesh closer to the nucleus)
for each atom.

###  III.2 Numerical implementation

####  III.2.1 Wave-function representations

GPAW supports three representations for smooth wave functions. The plane wave
(PW)

$\tilde{\psi}_{{\mathbf{k}}n}(r)=\sum_{\mathbf{G}}C_{n{\mathbf{G}}}e^{i({\mathbf{k}+\mathbf{G}})\cdot\mathbf{r}}$
(13)

and linear combination of atomic orbitals (LCAO)

$\tilde{\psi}_{{\mathbf{k}}n}(r)=\sum_{\mu}C_{n{\mathbf{k}}\mu}\Phi_{\mathbf{k}\mu}(\mathbf{r})$
(14)

representations rely on basis functions. The finite difference (FD) mode
relies on a representation of kinetic energy operator on a uniform Cartesian
grid.

####  III.2.2 All-electron quantities

The beauty of the PAW method is that you never _need_ to transform the pseudo
wave functions to all-electron wave functions, but you _can_ do it if you want
to. GPAW has tools for interpolating the pseudo wave functions to a fine real-
space grid and adding the PAW corrections. A fine real-space grid is needed to
properly represent the cusp and all the oscillations necessary for the wave
function to be orthogonal to all the frozen core states.

GPAW also has tools for calculating the all-electron electrostatic potential.
This is useful for transmission electron microscopy (TEM) simulations [70].
Most TEM simulations have relied on the so-called independent atom model
(IAM), where the specimen potential is described as a superposition of
isolated atomic potentials. While this is often sufficient, there is
increasing interest in understanding the influence of valence bonding [71].
This can be investigated by a TEM simulation code such as _ab_ TEM [72], which
can directly use _ab initio_ scattering potentials from GPAW.

####  III.2.3 Solving the Kohn–Sham equation

The default method for solving the Kohn–Sham equation for PW and FD modes is
to do iterative diagonalization combined with density mixing; for LCAO mode we
do a full diagonalization of the Hamiltonian. Alternatively, one can do direct
minimization as described in the next section.

For PW and FD modes, we need an initial guess of the wave functions. For this,
we calculate the effective potential from a superposition of atomic densities
and diagnonalize an LCAO Hamiltonian in a small basis set consisting of all
the pseudo partial waves corresponding to bound atomic states.

Each step in the self-consistent field (SCF) loop consists of the following
operations: 1) Diagonalization of the Hamiltonian in the subspace of the
current wave functions (skipped for LCAO). 2) One or more steps through the
iterative eigensolver (except for LCAO, where a full diagonalization is
performed). 3) Update of eigenvalues and occupation numbers. 4) Density mixing
and symmetrization. See previous work [32, 33] and Ref. [73] for details.

GPAW has two kinds of Poisson equation solvers: direct solvers based on
Fourier transforms or Fourier-sine transforms, and iterative multi-grid
solvers (Jacobi or Gauss-Seidel). The default is to use a direct solver
whereas the iterative solvers may be chosen for larger systems where they can
be more efficient.

For 0-, 1- and 2-dimensional systems, the default boundary conditions are to
have the potential go to zero on the cell boundaries. This becomes a problem
for systems involving large dipole moments. The potential due to the dipole is
long-ranged and, thus, the converged potential requires large vacuum sizes.
For molecules (0D systems), the boundary conditions can be improved by adding
multipole moment corrections to the density so that the corresponding
multipoles of the density vanish. The potential of these corrections is added
to the obtained potential. The same trick is used to handle charged systems.
For slabs (2D systems), a dipole layer can be added to account for differences
in the work functions on the two sides of the slab.

Methods for calculating occupation numbers are the Fermi-Dirac, Marzari-
Vanderbilt [74] and Methfessel-Paxton distributions as well as the tetrahedron
method and the improved tetrahedron method [75].

####  III.2.4 Updating wavefunctions in dynamics

Simulations commonly move the atoms without changing other parameters. If an
atom moves only slightly, we would expect most of the charge in its immediate
vicinity to move along with it. We use this to compute an improved guess for
the wave functions in the next self-consistency loop with FD or PW mode where
the eigensolver is iterative.

Near the atoms, the dual basis of pseudo-partial waves and projectors is
almost complete, i.e.

$\sum_{i}|\phi_{i}^{a}\rangle\langle\tilde{p}_{i}^{a}|=1\quad\textrm{(near
atom {\hbox{a}})}.$ (15)

If an atom moves by $\Delta\mathbf{R}^{a}$, the wave functions
$\tilde{\psi}_{n}(\mathbf{r})$ are updated by rigidly moving the projection
$\sum_{i}\tilde{\phi}_{i}^{a}(\mathbf{r})\langle\tilde{p}_{i}^{a}|\tilde{\psi}_{n}\rangle$
along with it, i.e.,

$\displaystyle\tilde{\psi}_{n}^{\mathrm{new}}(\mathbf{r})=\tilde{\psi}_{n}(\mathbf{r})$
$\displaystyle+$
$\displaystyle\sum_{ai}\tilde{\phi}_{i}^{a}(\mathbf{r}-\Delta\mathbf{R}^{a})\langle\tilde{p}_{i}^{a}|\tilde{\psi}_{n}\rangle$
(16) $\displaystyle-$
$\displaystyle\sum_{ai}\tilde{\phi}_{i}^{a}(\mathbf{r})\langle\tilde{p}_{i}^{a}|\tilde{\psi}_{n}\rangle.$

As the partial waves on different atoms are not orthonormal, this expression
generally “double-counts” contributions, resulting in wave functions that are
to some extent unphysical. Nevertheless, we have found that this simple method
achieves a significant speedup ($\sim$15% in realistic structure
optimisations) compared to not updating the wave functions.

The method could be further improved by using the LCAO basis set and the
overlap matrix to prevent double-counting.

####  III.2.5 Direct minimization

Direct orbital minimization [76, 77, 78, 79] is a robust alternative to the
conventional eigensolver and density mixing routines. The orbitals can be
expressed as a unitary transformation of a set of reference, or auxiliary,
orbitals $\mathbf{\Psi}_{0}$:

$\mathbf{\Psi}=\mathbf{U}\mathbf{\Psi}_{0}$ (17)

In the direct minimization (DM) method implemented in GPAW [80, 50], the
unitary matrix $\mathbf{U}$ is parametrized as an exponential transformation,
i.e. $\mathbf{U}=e^{\mathbf{A}}$, where ${\mathbf{A}}$ is an anti-Hermitian
matrix (${\mathbf{A}}=-{\mathbf{A}^{\dagger}}$). The energy can be considered
as a functional of both $\mathbf{A}$ and $\mathbf{\Psi}_{0}$:

$E[\mathbf{\Psi}]=\mathcal{F}[\mathbf{A},\mathbf{\Psi}_{0}]$ (18)

Thus, in general, the optimal orbitals corresponding to the minimum of the
energy functional can be found in a dual loop procedure. First, the energy is
minimized with respect to the elements of $\mathbf{A}$, and second, the
functional $\mathcal{L}[\mathbf{\Psi}_{0}]=\underset{\mathbf{A}}{\rm
min\,}\mathcal{F}[\mathbf{A},\mathbf{\Psi}_{0}]$ is minimized with respect to
$\mathbf{\Psi}_{0}$:

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

Since anti-Hermitian matrices form a linear space, the inner loop minimization
can use well-established local minimization strategies such as efficient
quasi-Newton methods with inexact line search, e.g. the limited-memory
Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm. The outer loop
minimization follows the gradient
$\partial\mathcal{L}/\partial\mathbf{\Psi}_{0}$ projected on the tangent space
at $\mathbf{\Psi}_{0}$.

The GPAW formulation of DM is applicable with all representations of the
orbitals available in GPAW, as well as Kohn–Sham (unitary invariant) and
orbital density dependent (nonunitary invariant) energy functionals, and can
be used for both finite and extended systems. In LCAO calculations, the
reference orbitals are expressed as a linear combination of the atomic basis
functions $\mathbf{\Phi}$, $\mathbf{\Psi}_{0}=\mathbf{\Phi}\mathbf{C}_{0}$,
where the matrix of coefficients $\mathbf{C}_{0}$ is fixed. Therefore, only a
minimization with respect to the elements of the matrix $\mathbf{A}$ is
required [80]. For plane wave and real-space grid representations, a
minimization in the space tangent to the reference orbitals is sufficient if
the functional is unitary invariant. Otherwise, if the functional is
nonunitary invariant, such as when self-interaction correction is used (see
section III.3.5), an inner loop minimization in the occupied–occupied block of
the matrix $\mathbf{A}$ is performed to make the energy stationary with
respect to unitary transfomation of the occupied orbitals [50].

The DM method avoids diagonalization of the Hamiltonian at each step and as a
result it usually involves a smaller computational effort. The DM method has
also been shown to be more robust than conventional eigensolvers and density
mixing in calculations of molecules and extended systems [80]. However, the
current implementation does not support finite-temperature distribution of
occupation numbers and thus can only be used for systems with a finite band
gap.

####  III.2.6 Convergence criteria

The modular architecture of GPAW allows the user to have precise control over
how the SCF loop decides that the electronic structure has converged to
sufficient precision. GPAW contains simple keywords for common convergence
criteria, such as “energy”, “forces”, (electron) “density”, and “work
function”, which are sufficient for the most common use cases.

Internally, all convergence criteria are instances of convergence classes, and
for each step through the SCF loop, each convergence class is called. When any
convergence class is called, it is passed a context that contains the current
state of the calculation, such as the wavefunctions and the Hamiltonian. The
criterion can thus pull relevant data from the calculation to decide if it is
converged. Because the convergence criterion itself is an object, it can store
information, such as previous values of the energy for comparison to the new
value. When all convergence criteria report that they are converged, the
calculation as a whole is considered to be converged and terminates.

This modular nature gives the user full control over how each convergence
criterion operates. For example, the user can easily ask the energy criterion
to check the differences in the last four values of the energy, rather than
the last three. If a convergence criterion itself is expensive to calculate,
it can make sense to not check it until the rest of the convergence criteria
are met. This can be accomplished by activating an internal “calculate_last”
flag within the convergence criterion.

Users can easily add their own custom convergence criteria to the SCF loop. If
a user would like to use a criterion not included by default with GPAW, it is
straightforward to write a new criterion as a Python class, and pass this to
the convergence dictionary of GPAW. For example, if one wanted to be sure the
bandgap of a semiconductor was converged, the criterion could check the
bandgap at each iteration and compare it to stored values from previous
iterations, and report that the calculation is converged when the peak-to-peak
variance among the last $n$ iterations is below a given threshold.

####  III.2.7 PAW data sets and pseudopotentials

GPAW can read PAW datasets from (possibly compressed) XML-files following the
PAW-XML specification [81]. Dataset files for most of the periodic table can
be downloaded from the GPAW webpage or installed with the gpaw install-data
command-line tool. The datasets are available for the LDA, PBE, revPBE, RPBE
and GLLBSC functionals. The electronic structure code Abinit [82] also reads
the PAW-XML format allowing GPAW and Abinit to share PAW dataset collections
such as the Jollet-Torrent-Holzwarth collection [83].

Specialized datasets can be generated with the gpaw dataset command-line tool.
This allows one to tweak the properties of a dataset. Some examples could be:
1) add more/less semi-core states; 2) increase/decrease the augmentation
sphere radius to make the pseudo wave functions more/less smooth; 3)
add/remove projector functions and corresponding pseudo and all-electron
partial waves; or 4) base the PAW dataset on a different XC-functional. These
changes will affect the accuracy and cost of the calculations.

GPAW is also able to use norm-conserving pseudopotentials (NCPP) such as HGH
[84] and pseudopotentials in the UPF format, such as SG15 [85]. Non-local
NCPPs can be considered an approximation to PAW: in the PAW description, the
non-local part of the Hamiltonian (the term containing $\Delta H_{ij}^{a}$ in
Eq. (6)) will adapt to the environment, whereas for a NCPP, $\Delta
H_{ij}^{a}$ will be diagonal and have a fixed value taken from a reference
atom. Because of the norm-conservation, NCPPs will have $\Delta S_{ij}^{a}=0$.

####  III.2.8 Parallelization

GPAW can parallelize over various degrees of freedoms depending on the type of
calculation, and implements multiple algorithms for achieving good parallel
performance and scalability. In calculations involving k-points,
parallelization over them is typically the most efficient, as little
communication is needed during the summation of wave functions to calculate
the density and any derived quantities. As the number of k-points is often
limited especially in large systems, parallelization is also possible over
real-space grids in the FD and LCAO modes, as well as over plane waves in the
PW mode. All modes also support parallelization over electronic bands, which
is particularly efficient for real-time TDDFT where the time-propagation of
each band can be carried out independently. Additional parallelization
possibilities exist depending on calculations, such as over electron–hole
pairs in linear-respose TDDFT calculations.

Parallelization is done mainly with MPI. In the FD and LCAO modes, it is
possible in addition to use also OpenMP within shared-memory nodes, which can
improve performance when the number of CPU cores per node is large. Dense
linear algebra, such as matrix diagonalization and Cholesky decomposition, can
be carried out with the parallel ScaLAPACK or ELPA libraries. This applies to
both the direct diagonalization in the LCAO mode as well as to subspace
diagonalizations in iterative Davidson and RMM-DIIS methods in the FD and PW
modes.

For ground-state calculations, GPAW will divide the cores into three MPI-
communicators: k-points, bands and domain. When parallelizing over k-points
and/or bands, all the cores of the k-point and/or band communicators will have
a copy of the density (possibly distributed over the domain communicator).
GPAW has the option to redistribute the density from the domain-communicator
to all cores so that operations such as evaluating the XC-energy and solving
the Poisson equation can be done more efficiently.

For each k-point, all the cores in the band and domain communicators will
cooperate on calculating matrix elements of the Hamiltonian and the overlap
operators. Dense-matrix linear-algebra operations on those matrices can be
done on a single core (most efficient for small systems) or with ScaLAPACK
where the matrices are distributed over either all or some of the cores from
the pool of cores in the band and domain communicators.

One drawback of Python is that in large parallel calculations its import
mechanism may incur a heavy load on the filesystem because all the parallel
processes are trying to read the same Python files. GPAW tries to alleviate
this by a special “broadcast imports” mechanism: during initial module imports
only a single process loads the modules, afterwards MPI is used to broadcast
the data to all the processes.

Parallel scalability depends strongly on the calculation mode and the system.
The FD mode offers the best scalability for high core counts as only nearest-
neighbour communication is needed across domains over a domain decomposition.
In PW mode, the limiting factor is all-to-all communication in parallelization
over plane waves. In LCAO mode, communications arise from multi-center
integrals of basis functions across domains. At best, GPAW scales to tens or
hundreds of nodes in supercomputers.

####  III.2.9 GPU implementation

The GPU implementation of GPAW works both on NVIDIA and AMD GPUs, targeting
either CUDA or HIP backends, respectively. GPAW uses a combination of manually
written GPU kernels, external GPU libraries (such as cuBLAS / hipBLAS), and
CuPy [86]. CuPy offers an easy-to-use Python interface for GPUs centered
around a NumPy-like GPU array and makes many hardware details completely
transparent to the end-user.

In the manually written GPU kernels, both GPU backends (CUDA and HIP) are
targeted using a header-only porting approach [87] in which generic GPU
identifiers are translated to vendor-specific identifiers at compile time. For
example, to allocate GPU memory, the identifier gpuMalloc is used in the code,
that is then translated either to cudaMalloc or hipMalloc depending on which
GPU backend is targeted. This allows us to avoid unnecessary code duplication
and still target multiple hardware platforms natively.

An earlier GPU implementation of GPAW [88, 89] served as the starting point
for the recent work on a new GPU code based on the rewritten ground-state
code. The objects that store quantities like $\tilde{\psi}_{n}(\mathbf{r})$,
$\tilde{p}_{i}^{a}(\mathbf{r})$,
$\langle\tilde{p}_{i}^{a}|\tilde{\psi}_{n}\rangle$, $D_{ii^{\prime}}^{a}$,
$\tilde{n}(\mathbf{r})$ and $\tilde{v}(\mathbf{r})$, use Numpy arrays for the
CPU code and CuPy arrays when running on a GPU. At the moment, GPUs can be
used for total-energy calculations with LDA/GGA in the PW mode.

Parallelization to multiple GPUs is done using MPI. Each MPI rank is assigned
a single GPU and communication between the GPUs is handled by MPI. Support for
GPU-aware MPI makes it possible to do direct GPU-to-GPU communication without
unnecessary memory copies between the GPU device and the host CPU.

###  III.3 XC-functionals

Exchange–correlation (XC) functionals provide a mapping between the
interacting and the non-interacting system of electrons. In Kohn–Sham DFT, the
density is built from a set of occupied non-interacting single-particle
orbitals $\psi_{i}$: $n=\sum_{i}f_{i}\left|\psi_{i}\right|^{2}$, with $f_{i}$
denoting the occupation numbers, leading to the same density as the
interacting system. The total energy in DFT is expressed as a sum of density
functionals for the different contributions:

$\displaystyle
E_{\mathrm{tot}}[n]=T_{\mathrm{S}}[n]+V_{\mathrm{ext}}[n]+U_{\mathrm{H}}[n]+E_{\mathrm{xc}}[n],$
(20)

where $T_{\mathrm{S}}[n]$ denotes the kinetic energy of the non-interacting
system, $V_{\mathrm{ext}}[n]$ the energy of the density in the external
potential, $U_{\mathrm{H}}[n]$ the classical Coulomb energy of the density
with itself, and $E_{\mathrm{xc}}[n]$ the so called exchange–correlation
energy, which collects all energy contributions missing in the prior terms and
therefore provides a mapping between the interacting and the non-interacting
system of electrons. While the first three terms can be calculated exactly,
even the form of $E_{\text{xc}}$ is unknown and although proven to exist and
be exact in principle, it has to be approximated in practice. A huge number of
approaches belonging to several families exists [1, 2, 90]. Several of these
approximations are available in GPAW, and an overview is given in the
following.

####  III.3.1 Libxc and Libvdwxc

The _libxc_ library [65] provides implementations of several (semi-)local
variants of the XC functional, given by the LDA, GGA, and MGGA families. These
are available in GPAW by a combination of their names from _libxc_ as, e.g.
“GGA_X_PBE+GGA_C_PBE” for PBE [3]. Additionally, GPAW provides its own
implementation of several (semi-)local functionals called by their short-
names, e.g. TPSS [91], PBE [3] and LDA, the latter with the correlation of
Perdew and Wang [92]. Several hybrids, see below, are implemented in GPAW with
the support of the _libxc_ library for their local parts.

For fully non-local van der Waals functionals, like the vdW-DF functional
[31], GPAW uses the efficient fast Fourier-transform convolution algorithm by
Roman-Perez and Soler [93] as implemented in the _libvdwxc_ library [66].

####  III.3.2 GLLB-sc

GPAW has an implementation of the solid-state modification of the
Gritsenko–van Leeuwen–van Lenthe–Baerends exchange–correlation model potential
(GLLB-sc) [94, 95]. This potential has been shown to improve the description
of the band gap [16] and the $d$-electron states in noble metals [52, 96, 97,
98].

####  III.3.3 Hubbard U

DFT+U calculations using a Hubbard-like correction can be performed in GPAW to
improve the treatment of Coulombic interactions of localized electrons. This
correction is most commonly applied to the valence orbitals of transition
metals to assist with obtaining experimental band gaps of oxides [99], which
may otherwise be underestimated [100, 101]. Also, formation energies and
magnetic states are often improved due to a more correct description of the
electronic structure. It may also be applied to main-group elements such as N,
O, P, or S, but this is less commonly done [102].

In the formalism chosen in GPAW [103], one uses a single $U^{\text{eff}}$
parameter rather than seperate $U$ and $J$ parameters for on-site Coulombic
and on-site exchange, respectively [104]. The correction influences the
calculation by applying an energy penalty to the system:

$U^{\text{eff}}=U-J,$ (21)
$E_{\text{DFT+U}}=E_{\text{DFT}}+\sum_{a,i}\frac{U^{\text{eff}}_{i}}{2}\textup{Tr}\left(\rho^{a}_{i}-\rho^{a}_{i}\rho^{a}_{i}\right),$
(22)

where the sum runs over the atoms $a$ and orbitals $i$ for which the
correction should be applied. $E_{\text{DFT}}$ is calculated by a standard
GPAW calculation and is corrected to $E_{\text{DFT+U}}$ by penalizing the
energy such that fully occupied or fully unoccupied orbitals are stabilized.
The magnitude of the correction depends on $U^{\text{eff}}_{i}$, and the
atomic orbital occupation matrix ($\rho^{a}_{i}$) controls which orbitals
contribute to the correction based on their occupation. In principle, any
orbital on any element that is partially occupied can be corrected.

GPAW supports normalization of the occupation matrix, accounting for the
truncated portion of the wave function outside the augmentation sphere. To
maintain consistency with other codes that do not support normalization, this
normalization can be disabled, but large disagreements are expected when
applied to $p$ orbitals. GPAW is one of the few codes that currently supports
multiple simultaneous corrections to orbitals on the same atom; this is useful
when two types of orbitals, such as $p$ and $d$ orbitals, are nearly
degenerate but both are partially occupied.

There is no $U^{\text{eff}}$ that is strictly correct, but methods such as RPA
[105] or linear response [106, 107] can allow for a value to be calculated
from first principles. More commonly, $U^{\text{eff}}$ is chosen semi-
empirically to fit experimental properties such as formation energy [108],
band gap [109], or more recently by machine learning predictions [110].

####  III.3.4 Hybrids

Hybrid functionals, especially range-separated functionals (see below),
correct problems present in calculations utilizing (semi-)local functionals
such as the wrong asymptotic behavior of the effective potential leading to an
improper description of Rydberg excitations [111], the improper calculation of
the total energy against fractional charges leading to the charge
delocalization error [112], and the wrong description of (long-range) charge
transfer excitations due to the locality of the exchange hole [113, 114, 115].

The exchange–correlation energy $E_{\mathrm{xc}}$, can be split into the
contributions from exchange, $E_{\mathrm{x}}$, and from correlation,
$E_{\mathrm{c}}$ [2, 116]. Hybrid functionals combine the exchange from the
(semi-)local functionals with the exchange from Hartree–Fock theory (HF).
Global hybrids such as PBE0 [117] mix the exchange from DFT with the exchange
from HF by a global fixed amount. Range-separated functionals (RSF) add a
separation function $\omega_{\mathrm{RSF}}$ [118, 119] to express the Coulomb
kernel in the exchange integrals as

$\displaystyle\frac{1}{r_{12}}$
$\displaystyle=\underbrace{\frac{1-\left[\alpha+\beta\left(1-\omega_{\mathrm{RSF}}\left(\gamma,r_{12}\right)\right)\right]}{r_{12}}}_{\mathrm{SR}}$
$\displaystyle+\underbrace{\frac{\left[\alpha+\beta\left(1-\omega_{\mathrm{RSF}}\left(\gamma,r_{12}\right)\right)\right]}{r_{12}}}_{\mathrm{LR}},$
(23)

where $r_{12}=|\vec{r}_{1}-\vec{r}_{2}|$. Here $\alpha$ and $\beta$ are mixing
parameters for global and range-separated mixing, respectively.
$\omega_{\mathrm{RSF}}$ is a soft function with values ranging from one for
$r_{12}=0$ to zero for $r_{12}\to\infty$, where the decay is controlled by the
parameter $\gamma$. Long-range-corrected RSFs such as LCY-PBE [120] use
(semi-)local approximations for short-range (SR) interaction and apply HF
exchange for long-range (LR) interaction. Short-range-corrected RSFs such as
HSE [5] reverse this approach. The parameter $\gamma$ is either fixed or can
be varied to match criteria of the ideal functional, e.g. that the energy of
the highest occupied molecular orbital matches the ionization potential [121,
116].

Details on the FD-mode implementation of long-range corrected RSF can be found
in Refs. [122, 123]. In general, the FD-mode implementation of hybrids is
limited to molecules, and forces have not been implemented. The PW-mode
implementation of hybrids handles $\mathbf{k}$-points, exploits symmetries and
can calculate forces.

####  III.3.5 SIC

A fully self-consistent and variational implementation of the Perdew-Zunger
self-interaction correction [124] (PZ-SIC) is available in GPAW. It corrects
the various problems with (semi-)local Kohn–Sham functionals mentioned above
in the context of hybrid functionals. Atomic forces are available in the GPAW
implementation with all three types of basis sets. A corrected KS functional
has the form

$E^{\mathrm{SIC}}[n_{1},\dots,n_{n}]=E^{\mathrm{KS}}[n]-\alpha\sum_{i}^{\mathrm{occ}}\left(U_{\mathrm{H}}[n_{i}]+E_{\mathrm{xc}}[n_{i},0]\right)$
(24)

where the self-Coulomb and self-XC of each occupied orbital density is
subtracted from the Kohn–Sham energy functional. Due to the explicit
dependence on the orbital densities, the corrected energy functional is not
invariant under a unitary transformation among occupied orbitals and thereby
not a KS functional. As a result, the minimization of $E^{\mathrm{SIC}}$
requires special direct minimization techniques (see section III.2.5) and
delivers a specific set of (typically localized) optimal orbitals. The
calculations should be carried out using complex orbitals [125, 126, 127]. The
full PZ-SIC has been shown to give an over-correction to binding energy as
well as band gaps and improved results for these properties are obtained by
scaling the SIC by $\alpha=1/2$ [128, 129], while the long-range form of the
effective potential, necessary for Rydberg state calculations, requires the
full correction [130].

PZ-SIC has been shown to give accurate results in cases where commonly used KS
functionals fail. This includes, for example, the Mn dimer where the PBE
functional gives qualitatively incorrect results while the corrected
functional gives close agreement with high-level quantum chemistry results as
well as experimental measurements [131]. Another example is the defect state
of a substitutional Al atom in $\alpha$-quartz [132]. Also, PZ-SIC has been
shown to improve values of the excitation energy of molecules obtained in
variational calculations of excited states [50, 133] (see section VIII.2).

####  III.3.6 BEEF

A great strength of Kohn–Sham DFT and its extensions is that a reasonably high
accuracy of physical, material, and chemical properties can be obtained at a
relatively moderate computational cost. DFT is thus often used to simulate
materials, reactions, and properties, where there de facto exists no “better”
alternative. Even though more accurate electronic structure methods in
principle might exist for a given application, the poor scaling of
systematically more accurate methods often makes these computationally
infeasible at the given system size that is being studied using DFT. One is
therefore often in the situation that the accuracy of a given DFT calculation
of some materials or chemical property cannot be verified against e.g. a more
accurate solution of the Schrödinger equation, even on the biggest available
supercomputers.

On the other hand, the wealth of available XC functionals naturally allows one
to look at how sensitive a DFT result is on the choice of functional, and
often accuracy is therefore being judged primarily by applying a small set of
different XC functions, especially if no accurate benchmark theoretical
simulation or experimental measurement is available. A challenge, however, is
that the different available functionals are often known to be particularly
good at simulating certain properties and poor at others. It is thus not at
all clear how much one should trust a given functional for a given simulated
property. The Bayesian Error Estimation (BEE) class of functionals [134]
attempts to develop a practical framework to establishing an error estimate
from a selected set of functional “ingredients”.

Assume that the XC model $M(a)$ is a function of a set of parameters, $a$,
that can be varied freely. If a benchmark data set, $D$, of highly accurate
properties established from experiments or higher-accuracy electronic
structure simulations is available, we can attempt to identify the ensemble of
models given by some distribution function, $P(M(a))$, such that the most
likely model in the ensemble, $M(a_{0})$, makes an accurate predictions for
the benchmark data set, while the spread of ensemble reproduces the spread
between the predictions of the most likely model and the benchmark data.

Bayes’ theorem provides a natural framework to search for the ensemble
distribution. If a joint distribution between $M(a)$ and $D$ exists, which we
would assume, then Bayes’ theorm gives

$P(M(a)|D)\propto P(D|M(a))P_{0}(M(a)),$ (25)

where $P_{0}(M(a))$ is the prior expectation to the distribution of models
before looking at the data, $D$, and $P(D|M(a))$ is a likelihood of seeing the
data given the model. To achieve a useful ensemble, much care has to be put
into finding a large enough set of varied and accurate data for different
materials and chemical properties, and much care has to be applied as to how
the ensemble is regularized to avoid overfitting [135, 136].

In all the BEE functionals, choices are made such that ultimately the XC
functional is linear in $a$, and such that the distribution of $a$ ends up
following a multidimensional Normal distribution given by a regularized
covariance matrix $\Gamma$:

$P(a)\propto\exp(-(a-a_{0})\Gamma(a-a_{0})),$ (26)

where $\Gamma$ has been scaled in such a way that the ensemble reproduce the
observed standard deviation between $M(a_{0})$ and $D$.

Figure 2:  Bayesian ensemble of XC functionals around BEEF–vdW. Orange solid
line is the BEEF–vdW exchange enhancement factor, while the blue lines depict
$F_{\text{x}}(s)$ for 50 samples of the randomly generated ensemble. Dotted
black lines mark the exchange model perturbations that yield DFT results ±1
standard deviation away from BEEF–vdW results.

Figure 2 shows an ensemble of exchange enhancement factors $F_{x}(s)$ from the
BEEF-vdW functional, where $s=|\nabla n|/(2k_{F}n)$ is the reduced density
gradient. The approach has given rise to several functionals BEEF-vdW [135],
mBEEF [136], and mBEEF-vdW [137], which all include error estimation and are
readily available in GPAW. Through the ASE interface, one can for example
utilize the error ensembles to establish error estimates on Python-implemented
models that are using DFT simulations for their parametrization. An example of
this is the application to error estimates on adsorption scaling relations,
microkinetics, and materials selection [138].

One risk to the approach of establishing error estimates from a small selected
“ensemble” of XC-functionals is clearly that if the simulated property in
question is poorly described by all functionals in the ensemble, then the
error estimate might be also become poor. This could for example be the case
if one tried to establish an error estimate for band gaps in oxides or van der
Waals bonding of adsorbates on surfaces based on an ensemble of GGA XC-
functionals, since no GGA functional may be accurate for simulating either
property.

##  IV Ion dynamics

GPAW can be employed as a ‘black-box’ calculator, supplying energies and
forces to other programs such as e.g. ASE, which then perform optimization of
ground state geometries and reaction paths or carry out molecular dynamics. In
fact, this is a key design principle behind GPAW. Methodological developments
and general implementations that are not directly dependent on fast access to
detailed electronic structure information should preferably be implemented
externally to GPAW. This lead to maximal simplicity of the GPAW code itself,
while also allowing for the external code to be utilized and tested with other
electronic structure codes and atomistic potentials. A key to the high
efficiency of GPAW simulations involving ionic displacements is the versatile
implementation of constraints in ASE. Here many types of constraints are
readily accessible, from the simple removal of degrees of freedom to more
exotic constraints allowing for rigid molecule dynamics [139], to harmonic
restoring forces [140], spacegroup preservation and combined ionic–unit cell
dynamics [141]. Many algorithms are available for various structure
optimization and molecular dynamics tasks.

###  IV.1 Structure relaxation

Local structure optimization in GPAW is typically achieved through the use of
an optimizer from ASE. Here a larger range of standard optimizers are
available such as quasi-Newton algorithms including BFGS [142, 143, 144, 145]
and limited-memory BFGS [146], or Newtonian dynamics-based algorithms such as
MDMin and FIRE [147]. That the optimizers have been implemented externally to
GPAW provides benefits in terms of simple procedures for restarting long
simulations and monitoring their convergence. Some optimizers from SciPy [68]
are also available through the open-source ASE package, which provides a
simple way to interface any optimizer in the SciPy format to GPAW.
Preconditioning is implemented in an accessible way [148], which often leads
to significant performance improvements.

Of the classical optimization methods, the quasi–Newton algorithms are often
highly competitive. Here one builds up information on the Hessian or inverse
Hessian matrix from calculated forces, ultimately leading to an accurate
harmonic model of the potential energy surface in the vicinity of the local
minimum. Such algorithms can, however, have problems both dealing with
anharmonicity in the potential energy surface and with any noise in the
electronic structure simulations. It often makes sense instead to fit a
Gaussian process to the calculated energies and forces, and minimize within
this model potential. This is implemented as the so-called GPmin method in ASE
and often converges on the order of three times faster than the best
quasi–Newton optimizers. [149]

###  IV.2 Reaction paths and barriers

Reliable calculations of energy barriers are of key importance for determining
the rate of atomistic processes. In many quantum-chemical codes utilizing
accurate atom-centered basis functions, this is achieved using analytical
second derivatives in the direct search for first-order saddle points. This
approach is less useful in the plane-wave based or grid-based modes of GPAW.
The Dimer method [150] is implemented in ASE and can be used with GPAW, but
often one would like to have an overview of the entire reaction path of an
atomic-scale process to verify that the correct energy barrier for a process
has been determined, and to obtain a full overview of the atomistic mechanism.
For this purpose, the Nudged Elastic Band (NEB) method is typically employed
through the ASE interface to GPAW. Both the original method [151] and a range
of later improvements are available [152, 153, 154, 155, 156, 157]. Special
care has to be taken in selecting the optimizer for a NEB algorithm, as this
choice can have drastic influence on the convergence rate. For optimization of
the reaction path, drastic performance improvements can be obtained by
carrying out the optimization in a surrogate machine-learning model fitted to
the potential energy surface [158, 159, 160]. GPAW has been used to drive NEB
calculations on both surface systems [161, 162, 163, 164] as well as in
molecules [164, 165].

###  IV.3 Global structure optimization

GPAW is integrated with various tools for global optimization of structures,
compositions, as well as materials morphologies. Some quite generally
applicable global optimization tools available are Basin Hopping [166], Minima
Hopping [167], and Genetic Algorithms [168]. Some of the most powerful global
optimization problems addressed using GPAW rely on machine-learning
accelerated global-optimization strategies. These strategies have for example
been applied to surfaces, clusters [169], and crystal structures in general
[170]. In other machine-learning accelerated global-optimization routines,
GPAW was used to generate initial databases of surface and bulk systems, and
for later model validation in a strategy that uses Gaussian processes to
generate surrogate potential-energy surfaces. These were then explored with
Bayesian-optimization techniques, achieving speed-ups over more conventional
methods of several orders of magnitude in finding optimal structures of the
systems under investigation [171, 172, 173]. The method has been augmented
through introducing extra (hyper)dimensions by interpolating between chemical
elements, which speeds up the efficiency of the global search [174]. GPAW has
been integrated with a Covariance Matrix Adaptation Evolution Strategy (CMA-
ES) framework, providing energies and forces, generating training data, and
evaluating CMA-ES candidate structures [175].

###  IV.4 Molecular dynamics and QM/MM

_Ab initio_ molecular dynamics (MD) can be performed and analyzed through the
ASE interface to GPAW. This include standard simulations such as constant-NVE,
constant-NVT, and constant-NPT simulations. There is access to e.g. Langevin-,
Andersen-, Nose-Hoover-, and Berendsen-dynamics. Due to the externalization of
the dynamics, the development of novel algorithms and analysis tools becomes
facile. Some examples of the use of GPAW span _ab initio_ MD studies that have
explored the liquid structure of water [176] and the water/Au(111)
electrochemical interface [177].

Furthermore, GPAW is capable of working with external electrostatic potential
terms from user-supplied point charge values and positions, enabling Quantum
Mechanical (QM) / Molecular Mechanical (MM) simulations. Since GPAW is from
the outset designed around highly efficient grid operations, the computational
overhead of evaluating this potential as well as the resulting forces on the
MM point charges from the QM density is kept low [178]. GPAW has been central
in a range studies of ion dynamics in solution, both in and out of
equilibrium. By using the molecular-dynamics functionality of ASE, researchers
have performed QM/MM MD simulations of excited-state bond formation in
photocatalyst model systems [179, 180, 181, 182], and electron transfer as
well as coupled solute-solvent conformational dynamics in photosensitizer
systems [183, 165]. Work on polarizable embedding QM/MM within GPAW is ongoing
[184, 185].

##  V Magnetism and spin

Many important technological applications utilize magnetic order or
manipulation of spin in materials. GPAW has a wide range of functionalities
that faciliate the analysis of magnetic properties. This includes calculations
with non-collinear spin, inclusion of external magnetic fields, spin–orbit
coupling and spin spiral calculations within the generalized Bloch theorem.
The implementaton of these features is described below while additional
methods for calculating magnetic excitations are described in section VI.4.

###  V.1 Spin–orbit coupling

The spin–orbit coupling is typically completely dominated by regions close to
the nuclei where the electrostatic potential becomes strong. Within the PAW
augmentation sphere of atom $a$, the spin-$\sigma$ component of orbital $n$ is
given by

$\psi_{n\sigma}(\mathbf{r})=\sum_{i}\langle\tilde{p}^{a}_{i}|\tilde{\psi}_{n\sigma}\rangle\phi_{i}^{a}(\mathbf{r}).$
(27)

Assuming a spherically symmetric potential, the spin–orbit Hamiltonian for
atom $a$ is then written as

$\hat{H}_{\mathrm{SOC}}^{a}=\frac{\mu_{\mathrm{B}}}{\hbar
m_{\mathrm{e}}ec^{2}}\frac{1}{r}\frac{\mathrm{d}V_{r}^{a}}{\mathrm{d}r}\hat{\mathbf{L}}\cdot\hat{\mathbf{S}},$
(28)

where $V_{r}^{a}$ is the radial electrostatic potential of atom $a$. We
evaluate $V_{r}^{a}$ as the spherical part of the XC- and Hartree-potential
from the local expansion of the density given by Eq. 11. Since the partial
waves $\phi^{a}_{i}$ are eigenstates of the scalar-relativistic Hamiltonian,
these are independent of spin and it is straightforward to evaluate the action
of $\hat{\mathbf{L}}$ on them.

The eigenenergies can be accurately calculated in a non-selfconsistent
treatment of the spin–orbit coupling and may be obtained by diagonalizing the
full Hamiltonian in a basis of scalar-relativistic orbitals [186]:

$H_{mn}=\varepsilon^{0}_{m}\delta_{mn}+\sum_{aij\sigma\sigma^{\prime}}\langle\tilde{\psi}_{m\sigma}^{0}|\tilde{p}^{a}_{i}\rangle\langle\phi_{i}^{a}|\hat{H}^{a}_{\mathrm{SOC,\sigma\sigma^{\prime}}}|\phi_{j}^{a}\rangle\langle\tilde{p}^{a}_{j}|\tilde{\psi}_{n\sigma^{\prime}}^{0}\rangle.$
(29)

Here $\varepsilon^{0}_{m}$ and $\psi_{m}^{0}$ represent the scalar-
relativistic eigenenergies and eigenstates, respectively. This constitutes a
fast post-processing step for any scalar-relativistic calculation and only
requires the projector overlaps
$\langle\tilde{p}^{a}_{j}|\tilde{\psi}_{n\sigma^{\prime}}^{0}\rangle$. It
should be noted that this approach in principle requires convergence with
respect to the number of scalar-relativistic states included in the basis, but
the eigenvalues typically converge rather rapidly with respect to the basis.
In Fig. 3, we show the band structure of a WS2 monolayer obtained from PBE
with non-self-consistent spin–orbit coupling. The W atoms introduce strong
spin–orbit coupling in this material and the valence band is split by 0.45 eV
at the K point. The spin degeneracy is retained along the $\Gamma$-M line,
which is left invariant by two non-commuting mirror symmetries.

Figure 3: Band structure of WS2 monolayer obtained from PBE with non-
selfconsistent spin–orbit coupling. The colors indicate the expectation value
of $S_{z}$ for each state. The grey lines show the band structure without
spin–orbit coupling.

For magnetic materials the non-selfconsistent treatment of spin–orbit coupling
is convenient for evaluating the magnetic anisotropy. The magnetic force
theorem [187] implies that rotating the magnetic moments away from the ground
state configuration yields a contribution to the energy, which is well
approximated by the change in Kohn–Sham eigenvalues. The change in energy for
a given orientation of the magnetization density can thus be obtained as

$\Delta
E^{\theta,\varphi}=\sum_{n}f_{n}^{\theta,\varphi}\varepsilon_{n}^{\theta,\varphi},$
(30)

where $\varepsilon_{n}^{\theta,\varphi}$ are the eigenvalues obtained from
diagonalizing Eq. (29) with the spins rotated to a direction defined by the
angles $(\theta,\varphi)$ and $f_{n}^{\theta,\varphi}$ are the associated
occupation numbers. This corresponds to rotating the xc-magnetic field
(defined below), which will lead to different eigenvalues when spin–orbit
coupling is included. In two-dimensional magnets, an easy-axis anisotropy is
decisive for magnetic order [188, 189, 41] and Eq. (30) is easily applied to
high throughput computations of magnetic properties [190, 191].

###  V.2 Self-consistent non-collinear magnetism

The Kohn–Sham framework for treating non-collinear magnetism was developed by
Barth and Hedin [192] and involves the spin density matrix

$\rho_{\sigma\sigma^{\prime}}(\mathbf{r})=\sum_{n}f_{n}\psi^{*}_{n\sigma}(\mathbf{r})\psi_{n\sigma^{\prime}}(\mathbf{r})$
(31)

as the basic variable. The electronic density and magnetization are then given
by $n=\mathrm{Tr}[\rho]$ and $\mathbf{m}=\mathrm{Tr}[\boldsymbol{\sigma}\rho]$
respectively. The XC-potential acquires four components, which are the
functional derivatives of the XC-energy with respect to the density matrix and
may be represented as a $2\times 2$ matrix acting on spinor Kohn–Sham states.
These can be expressed in terms of the density and magnetization, which lead
to the XC-part of the Kohn–Sham Hamiltonian

$H_{\mathrm{xc}}=v_{\mathrm{xc}}+\mathbf{B}_{\mathrm{xc}}\cdot\boldsymbol{\sigma},$
(32)

where the scalar potential and XC-magnetic field are given by

$v_{\mathrm{xc}}(\mathbf{r})=\frac{\delta E_{\mathrm{xc}}}{\delta
n(\mathbf{r})}\qquad\mathbf{B}_{\mathrm{xc}}(\mathbf{r})=\frac{\delta
E_{\mathrm{xc}}}{\delta\mathbf{m}(\mathbf{r})}.$ (33)

In GPAW, the self-consistent treatment of non-collinear spin is implemented
within LDA where $\mathbf{B}_{\mathrm{xc}}$ is approximated as

$\mathbf{B}_{\mathrm{xc}}^{\mathrm{LDA}}(\mathbf{r})=\frac{\delta
E_{\mathrm{xc}}^{\mathrm{LDA}}}{\delta
m(\mathbf{r})}\mathbf{\hat{m}}(\mathbf{r}).$ (34)

Here $m(\mathbf{r})=|\mathbf{m}(\mathbf{r})|$ is the magnitude and
$\mathbf{\hat{m}}(\mathbf{r})$ is the direction of the magnetization. The
generalization of Eq. 34 to GGAs are plagued by formal as well as numerical
problems [193, 194] and the self-consistent solution of the non-collinear
Kohn–Sham equations is presently restricted to LDA. It should be emphasized
that the PAW formalism allows for a fully non-collinear treatment that does
not rely on intra-atomic collinearity, which is often imposed in other
electronic structure packages.

Spin–orbit coupling may be included by adding Eq. 28 to the Kohn–Sham
Hamiltonian and this constitutes a fully self-consistent framework for
spin–orbit calculations within LDA.

###  V.3 Orbital magnetization

Current _ab initio_ methods for determining the orbital magnetization of a
material involve either the modern theory, i.e. a Berry-phase formula, or the
atom-centered approximation (ACA), where contributions to the expectation
value of the angular momentum operator are restricted within atom-centered
muffin-tin (MT) spheres with specified cutoff radii. The PAW formulation of
the wave functions allows for an approximation similar to the MT-ACA where
only the PAW expansion of the wave functions is assumed to contribute
significantly to the expectation value of the angular momentum. The orbital
magnetic moments can then be calculated through

$\mathbf{m}^{a}_{\mathrm{orb}}=-\frac{e}{2m_{e}}\sum_{a}\sum_{ii^{\prime}}\sum_{\sigma}D^{a}_{\sigma
ii^{\prime}}\left<\phi^{a}_{i}\middle|\widehat{\mathbf{L}}\middle|\phi^{a}_{i^{\prime}}\right>,$
(35)

where $D^{a}_{\sigma ii^{\prime}}$ are elements of the atomic density matrix
as defined in Eq. (9), $\phi^{a}_{i}(\mathbf{r})$ are the _bound_ all-electron
partial waves of atom $a$, and $\hat{\mathbf{L}}$ is the angular-momentum
operator.

We note that this expression does not entail cutoff radii for the atomic
contributions to the orbital magnetic moments; instead, the entire all-
electron partial waves are included although the all-electron atomic expansion
is only formally exact inside the PAW spheres. Additionally, this means that
the unbounded, i.e. non-normalisable, all-electron partial waves must be
excluded in Eq. (35)

A prerequisite for nonzero orbital magnetization is broken time-reversal
symmetry as represented by a complex Hamiltonian that is not unitarily
equivalent to a real counterpart. Practically this means that a finite orbital
magnetization requires magnetic order and either inclusion of spin–orbit
coupling or non-coplanar spin textures. In GPAW, the spin–orbit interaction
can be included either self-consistently in a non-collinear calculation or as
a non-self-consistent post-processing step following a collinear calculation.

Crystal | bcc-Fe | fcc-Ni | fcc-Co | hcp-Co
---|---|---|---|---
Easy axis | $[001]$ | $[111]$ | $[111]$ | $[001]$
PAW-ACA | $0.0611$ | $0.0546$ | $0.0845$ | $0.0886$
MT-ACA [195] | $0.0433$ | $0.0511$ | $0.0634$ | $0.0868$
Modern theory [195] | $0.0658$ | $0.0519$ | $0.0756$ | $0.0957$
Experiment [196] | $0.081$ | $0.053$ | $0.120$ | $0.133$
Table 2: Calculated and measured values for the orbital magnetization in units
of $\mu_{\mathrm{B}}$ per atom.

The orbital magnetization has been calculated using Eq. (35) for the simple
ferromagnets bcc-Fe, fcc-Ni, fcc-Co, and hcp-Co, where the spin–orbit
interaction is included non-self-consistently. The PAW-ACA results are
displayed in Table 2 along with MT-ACA results, modern theory results, and
measurements from experiment, demonstrating mainly that the PAW-ACA can be an
improvement over the MT-ACA and secondly that there is decent agreement
between the PAW-ACA and the modern theory for these systems.

###  V.4 Constant B-field

The coupling of the electronic spin magnetic moments to a constant external
magnetic field $\mathbf{B}$ can be included by the addition of a Zeeman term
in the Kohn–Sham Hamiltonian:

$\Delta H_{\mathrm{B}}=\mathbf{B}\cdot\boldsymbol{\sigma}.$ (36)

As an example, we can consider the spin-flop transition in Cr2O3. The ground
state is anti-ferromagnetic with a weak anisotropy that prefers alignment of
the spins with the $z$-axis, and if a magnetic field is applied along the
$z$-direction, it becomes favorable to align the spins along a perpendicular
direction (with a small ferromagnetic component) at the critical field where
the Zeeman energy overcomes the anisotropy. This is clearly a spin–orbit
effect and one has to carry out the calculations in a fully non-collinear
framework with self-consistent spin–orbit coupling. In Fig. 4 we show the
energy of the two alignments of spin as a function of external magnetic field
obtained with LDA+U ($U=2$ eV). The minimum energy configuration changes from
$S_{z}$ to $S_{x}$ alignment at 6 T, i.e., the spin flips at this critical
field value. This is in excellent agreement with the experimental value of
Ref. [197]. The critical field is, however, rather strongly dependent on the
chosen value of $U$ and increases by a factor of two in the absence of $U$.
The magnetic anisotropy (per unit cell) can be read off from the energy
difference at $B=0$ T.

Figure 4: Spin-flop transition in Cr2O3. The alignment of the spins with
respect to the magnetic field is sketched for small and large magnitudes of
the field. The canting (small ferromagnetic component after the spin flop) is
exaggerated for visualization. The actual canting is roughly $1^{\circ}$.

###  V.5 Spin spirals

The ground-state magnetic structure of frustrated and/or chiral magnets is
often non-collinear and may be incommensurate with the chemical unit cell. In
the classical isotropic Heisenberg model, the energy is always minimized by a
planar spin spiral [198]. The energy of a general planar spin spiral can be
evaluated efficiently within the chemical unit cell using the Generalized
Bloch’s Theorem (GBT) [199]. In the GBT implementation of GPAW [40], there is
no restriction on interatomic collinearity and thus it encodes a rotation of
the all-electron magnetization density by an angle
$\varphi=\mathbf{q}\cdot\mathbf{R}_{i}$ upon translations of a lattice vector
$\mathbf{R}_{i}$. The Kohn–Sham equations can then be solved self-consistently
for a fixed wave vector $\mathbf{q}$,

$\displaystyle\hat{H}_{\mathbf{q}}(\mathbf{k})\left|u_{\mathbf{q},n\mathbf{k}}\right\rangle=\varepsilon_{\mathbf{q},n\mathbf{k}}\left|u_{\mathbf{q},n\mathbf{k}}\right\rangle,$
(37)

using only the periodic part of the generalized Bloch orbitals,
$|\psi_{\mathbf{q},n\mathbf{k}}\rangle=U_{\mathbf{q}}^{\dagger}(\mathbf{r})e^{i\mathbf{k}\cdot\mathbf{r}}|u_{\mathbf{q},n\mathbf{k}}\rangle$,
where $U_{\mathbf{q}}$ denotes the generator of spin rotations about the
$z$-axis:

$\displaystyle
U_{\mathbf{q}}(\mathbf{r})=\begin{pmatrix}e^{i\mathbf{q}\cdot\mathbf{r}/2}&0\\\
0&e^{-i\mathbf{q}\cdot\mathbf{r}/2}\end{pmatrix}.$ (38)

The generalized Bloch Hamiltonian without spin–orbit coupling is given by

$\displaystyle\hat{H}_{\mathbf{q}}(\mathbf{k})$
$\displaystyle=e^{-i\mathbf{k}\cdot\mathbf{r}}U_{\mathbf{q}}(\mathbf{r})\hat{H}U_{\mathbf{q}}^{\dagger}(\mathbf{r})e^{i\mathbf{k}\cdot\mathbf{r}},$
(39)

and once Eq. (37) has been solved self-consistently for a given $\mathbf{q}$,
the corresponding spin-spiral energy is evaluated as usual,
$E_{\mathrm{SS}}(\mathbf{q})=E_{\mathbf{q}}^{\mathrm{DFT}}[n,\mathbf{m}]$.

Figure 5:  LSDA spin-spiral spectrum of the NiBr2 monolayer (structure taken
from the C2DB [15]). The ground state has an incommensurate wave vector
$\mathbf{Q}\simeq[0.1,0.1,0]$. The magnetic moment displays only weak
longitudinal fluctuations and the band gap remains finite for all wave vectors
$\mathbf{q}$, indicating that a Heisenberg model description of the material
would be appropriate. The inset shows the spin–orbit correction to the spin-
spiral energies (in meV) as a function of the normal vector $\mathbf{n}$ of
the planar spin spiral. $\mathbf{n}$ is depicted in terms of its stereographic
projection in the upper hemisphere above the monolayer plane. The spiral plane
is found to be orthogonal to $\mathbf{Q}$ and tilted slightly with respect to
the out-of-plane direction.

Using the GPAW functionality to compute the spin-spiral energy as a function
of $\mathbf{q}$, one can then search for the ground-state wave vector
$\mathbf{Q}$ that minimizes the energy. In Fig. 5, we show that the monolayer
NiBr2 has an incommensurate spin-spiral ground state, and that the local
magnetic moment on the Ni atom depends only weakly on the wave vector
$\mathbf{q}$.

Furthermore, the orientation of the ground state spin spiral can be obtained
by including spin–orbit coupling non-self-consistently in the projected
spin–orbit approximation [200] and search for the orientation that minimizes
the energy. The ordering vector $\mathbf{Q}$ and the normal vector
$\mathbf{n}$ of the spiral plane thus constitute a complete specification of
the magnetic ground state within the class of single-$\mathbf{q}$ states. The
normal vector is of particular interest since spin spirals may lead to
spontaneous breaking of the crystal symmetry, and the normal vector largely
determines the direction of magnetically induced spontaneous polarization
[40]. In Fig. 5, we show that $\mathbf{n}$ is perpendicular to the wavevector
$\mathbf{Q}$ in the monolayer NiBr2 ground state, corresponding to a cycloidal
spin spiral.

##  VI Response functions and excitations

Linear response functions are the bread and butter of condensed matter
physics. Their applications include the description of dielectric screening,
optical and electron energy loss spectra, many-body excitations, and ground
state correlation energies. In this section we describe the methods available
in GPAW for calculating electronic response functions as well as GW
quasiparticle band structures and optical excitations from the Bethe-Salpeter
Equation (BSE). In addition to electronic response function, GPAW can also
calculate the transverse magnetic susceptibility, which holds information
about the magnetic excitations, e.g. magnons, and can be used to derive
parameters for classical Heisenberg spin models and to estimate magnetic
transition temperatures. Finally, we present methods to calculate Raman
spectra of solids and quadratic optical response tensors for describing second
harmonics generation and the Pockels electro-optical effect.

###  VI.1 Linear-response TDDFT

To linear order, the change in electron density induced by a time-dependent
external (scalar) potential $V_{\mathrm{ext}}(\mathbf{r},t)$ is governed by
the electronic susceptibility, $\chi$:

$\delta n(\mathbf{r},t)=\int_{-\infty}^{\infty}dt^{\prime}\int
d\mathbf{r}^{\prime}\,\chi(\mathbf{r},\mathbf{r}^{\prime},t-t^{\prime})V_{\mathrm{ext}}(\mathbf{r}^{\prime},t^{\prime}).$
(40)

The susceptibility is itself given by the Kubo formula [201]

$\chi(\mathbf{r},\mathbf{r}^{\prime},t-t^{\prime})=-\frac{i}{\hbar}\theta(t-t^{\prime})\langle\left[\hat{n}_{0}(\mathbf{r},t),\hat{n}_{0}(\mathbf{r}^{\prime},t^{\prime})\right]\rangle_{0},$
(41)

where the expectation value is taken with respect to the ground state at zero
temperature and the density operators are cast in the interaction picture.

In the noninteracting Kohn–Sham system, the susceptibility can be evaluated
explicitly from the Kohn–Sham orbitals, eigenvalues and occupations:

$\displaystyle\chi_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)=$
$\displaystyle\lim_{\eta\rightarrow
0^{+}}2\,\sum_{n\mathbf{k}}\sum_{m\mathbf{k}^{\prime}}\left(f_{n\mathbf{k}}-f_{m\mathbf{k}^{\prime}}\right)$
$\displaystyle\times\frac{\psi^{*}_{n\mathbf{k}}(\mathbf{r})\psi_{m\mathbf{k}^{\prime}}(\mathbf{r})\psi^{*}_{m\mathbf{k}^{\prime}}(\mathbf{r}^{\prime})\psi_{n\mathbf{k}}(\mathbf{r}^{\prime})}{\hbar\omega-(\epsilon_{m\mathbf{k}^{\prime}}-\epsilon_{n\mathbf{k}})+i\hbar\eta}.$
(42)

Based on the Kohn–Sham susceptibility, $\chi_{0}$, the many-body
susceptibility can be calculated via a Dyson-like equation [202]:

$\displaystyle\chi(\mathbf{r},\mathbf{r}^{\prime},\omega)=\chi_{0}(\mathbf{r},\mathbf{r}^{\prime},\omega)+$
$\displaystyle\iint
d\mathbf{r}_{1}d\mathbf{r}_{2}\,\chi_{0}(\mathbf{r},\mathbf{r}_{1},\omega)$
$\displaystyle\times
K_{\mathrm{Hxc}}(\mathbf{r}_{1},\mathbf{r}_{2},\omega)\chi(\mathbf{r}_{2},\mathbf{r}^{\prime},\omega).$
(43)

Here, electronic interactions are accounted for via the Hartree-XC kernel,
which is defined in terms of the effective potentials of time-dependent DFT
(TDDFT) [203]:

$K_{\mathrm{Hxc}}(\mathbf{r},\mathbf{r}^{\prime},t-t^{\prime})=v_{c}(\mathbf{r}-\mathbf{r}^{\prime})+f_{\mathrm{xc}}(\mathbf{r},\mathbf{r}^{\prime},t-t^{\prime})$
(44)

where $v_{c}$ is the Coulomb interaction and $f_{\mathrm{xc}}$ is the XC-
kernel defined as

$f_{\mathrm{xc}}(\mathbf{r},\mathbf{r}^{\prime},t-t^{\prime})=\frac{\delta
v_{\mathrm{xc}}(\mathbf{r},t)}{\delta n(\mathbf{r}^{\prime},t^{\prime})}.$
(45)

In prototypical linear-response TDDFT (LR-TDDFT) calculations, the
exchange–correlation part of the kernel is either neglected (leading to the
random phase approximation (RPA)) or approximated via the adiabatic local
density approximation (ALDA),

$v_{\mathrm{xc}}(\mathbf{r},t)\simeq\left.\frac{\partial[\epsilon_{\mathrm{xc}}(n)n]}{\partial
n}\right|_{n(\mathbf{r},t)},$ (46)

where $\epsilon_{\mathrm{xc}}(n)$ denotes the XC-energy per electron of the
homogeneous electron gas of density $n$. The ALDA is in fact rather restricted
as it only ensures a correct description of the kernel for metals in the long
wave length limit. This implies that it cannot account for excitons in
extended systems and furthermore it leads to a divergent XC-hole. The latter
problem can be resolved by a simple renormalization of the ALDA, that
regulates the ontop XC-hole and drastically improves the description of local
correlations over the ALDA (and RPA) [37].

####  VI.1.1 Implementation for periodic systems

In crystalline systems, the susceptibility is periodic with respect to
translations on the Bravais lattice,
$\chi(\mathbf{r}+\mathbf{R},\mathbf{r}^{\prime}+\mathbf{R},\omega)=\chi(\mathbf{r},\mathbf{r}^{\prime},\omega)$.
Consequently, the susceptibility can be Fourier-transformed according to

$\chi_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)=\iint\frac{d\mathbf{r}d\mathbf{r}^{\prime}}{\Omega}e^{-i(\mathbf{G}+\mathbf{q})\cdot\mathbf{r}}\chi(\mathbf{r},\mathbf{r}^{\prime},\omega)e^{i(\mathbf{G}^{\prime}+\mathbf{q})\cdot\mathbf{r}^{\prime}},$
(47)

translating the Dyson equation (43) into a plane-wave matrix equation, which
is diagonal in the wave vector $\mathbf{q}$ and can be inverted numerically.

By using a plane-wave representation, the crux of the implementation then
becomes to calculate the reciprocal-space pair densities

$n_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\mathbf{G}+\mathbf{q})=\langle\psi_{n\mathbf{k}}|e^{-i(\mathbf{G}+\mathbf{q})\cdot\mathbf{r}}|\psi_{m\mathbf{k}+\mathbf{q}}\rangle,$
(48)

and to Fourier transform the XC-kernel

$f_{\mathrm{LDA}}(\mathbf{G})=\int_{\Omega_{\mathrm{cell}}}d\mathbf{r}\,e^{-i\mathbf{G}\cdot\mathbf{r}}\left.\frac{\partial^{2}[\epsilon_{\mathrm{xc}}(n)n]}{\partial
n^{2}}\right|_{n(\mathbf{r})}.$ (49)

Crucially, both can be evaluated by adding a PAW correction to the analogous
pseudo quantity, which itself can be evaluated by means of a fast Fourier
transform (FFT). Implementational details hereto are reported in Refs. [43]
and [36].

The Hartree contribution to the kernel (44) is simply given by the bare
Coulomb interaction (here given in atomic units)

$v_{\mathrm{c},\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q})=\frac{4\pi}{|\mathbf{G}+\mathbf{q}|^{2}}\delta_{\mathbf{G},\mathbf{G}^{\prime}},$
(50)

and can be evaluated analytically when $\mathbf{q}$ is finite. However, in the
optical limit, the $\mathbf{G}=\mathbf{0}$ component diverges and one needs to
be careful when inverting the Dyson equation (43). In GPAW, this is handled by
expanding the pair densities within $kp$-perturbation theory:

$\displaystyle\langle\psi_{n\mathbf{k}}|e^{-i(\mathbf{G}+\mathbf{q})\cdot\mathbf{r}}|\psi_{m\mathbf{k}+\mathbf{q}}\rangle_{\mathbf{q}\rightarrow\mathbf{0},\mathbf{G}=\mathbf{0}}=-i\mathbf{q}\cdot\frac{\langle\psi_{n\mathbf{k}}|\boldsymbol{\nabla}|\psi_{m\mathbf{k}}\rangle}{\epsilon_{m\mathbf{k}}-\epsilon_{n\mathbf{k}}}.$
(51)

In the expansion, the diverging Coulomb term is exactly cancelled by the
$\mathbf{q}$-dependence of the pair densities, so that the product
$\chi_{0}v_{\mathrm{c}}$ remains finite. For additional details, see Ref.
[36].

####  VI.1.2 Spectral representation

GPAW offers two different ways of dealing with the frequency dependence of the
Kohn–Sham susceptibility (42). One is to evaluate the expression explicitly
for the frequencies of interest. This is advantageous if one is interested in
a few specific frequencies. The other is to evaluate the associated spectral
function

$\displaystyle S_{0,\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)=$
$\displaystyle\frac{2}{\Omega}\sum_{\mathbf{k}}\sum_{n,m}(f_{n\mathbf{k}}-f_{m\mathbf{k}+\mathbf{q}})$
$\displaystyle\times
n_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\mathbf{G}+\mathbf{q})n_{m\mathbf{k}+\mathbf{q},n\mathbf{k}}(-\mathbf{G}^{\prime}-\mathbf{q})$
$\displaystyle\times\delta\left(\hbar\omega-\left[\epsilon_{m\mathbf{k}+\mathbf{q}}-\epsilon_{n\mathbf{k}}\right]\right)$
(52)

from which the susceptibility can be obtained via a Hilbert transform,

$\displaystyle\chi_{0,\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)=\lim_{\eta\rightarrow
0^{+}}\int_{-\infty}^{\infty}d\omega^{\prime}\,\frac{S_{0,\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega^{\prime})}{\omega^{\prime}-\omega+i\eta}.$
(53)

To converge the Hilbert transform, the spectral function is evaluated on a
nonlinear frequency grid spanning the range of eigenenergy differences
included in Eq. (52). Although the calculation becomes more memory-intensive
as a result, it is usually faster to compute $\chi_{0}$ via its spectral
function. Furthermore, since $S_{0}$ is linearly interpolatable, the
tetrahedron method can be employed to improve the convergence with respect to
$\mathbf{k}$-points.

###  VI.2 Dielectric function

The longitudinal part of the dielectric tensor is related to the
susceptibility $\chi$ as

$\displaystyle\varepsilon^{-1}_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)=\delta_{\mathbf{G}\mathbf{G}^{\prime}}+v_{c,\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q})\chi_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega).$
(54)

From the dielectric matrix, the macroscopic dielectric function, including
local field corrections, is given by

$\displaystyle\varepsilon_{M}(\mathbf{q},\omega)=\frac{1}{\varepsilon^{-1}_{00}(\mathbf{q},\omega)}$
(55)

from which it is straightforward to extract the optical absorption spectrum

$\displaystyle\mathrm{ABS}=\mathrm{Im}\>\varepsilon_{M}(\mathbf{q}\rightarrow
0,\omega)$ (56)

as well as the electron energy-loss spectrum

$\displaystyle\mathrm{EELS}=-\mathrm{Im}\>\frac{1}{\varepsilon_{M}(\mathbf{q},\omega)}.$
(57)

It is also possible to define a symmetrized version of the dielectric matrix
as

$\displaystyle\tilde{\varepsilon}_{\mathbf{GG}^{\prime}}(\mathbf{q},\omega)=v^{-1/2}_{\mathbf{G}}(\mathbf{q})\varepsilon_{\mathbf{GG}^{\prime}}(\omega,\mathbf{q})v^{1/2}_{\mathbf{G}^{\prime}}(\mathbf{q}).$
(58)

The susceptibility can be evaluated at the RPA level (by setting
$V_{\text{xc}}$ in Eq. 43 to zero) or using one of the XC-kernels implemented
in GPAW such as the the local ALDA, the non-local rALDA [37], or the bootstrap
kernel [204].

The evaluation $\chi_{0}$ is computationally demanding since it involves an
integration over the Brillouin zone (BZ) as well as a summation over occupied
and unoccupied states. The $\mathbf{k}$-point convergence can be increased
substantially with the tetrahedron method [205]. Contrary to the simple point
integration, where the $\delta$-function in Eq. 52 is replaced by a smeared-
out Lorentzian, the tetrahedron method utilizes linear interpolation of the
eigenvalues and matrix elements.

In Fig. 6 we compare the dielectric function computed using the two different
integration methods for two prototypical cases: (semiconducting) Si and
(metallic) Ag. Since Si is semiconducting there are no low-energy excitations
and consequently the imaginary part of the dielectric function is zero while
the real part is flat for low frequencies. The point integration and
tetrahedron integration yield the same value for $\omega\rightarrow 0$, but
the tetrahedron integration avoids the unphysical oscillations of
$\varepsilon$ at higher frequencies, exhibited by the point integration due to
the finite $\mathbf{k}$-point sampling. For metals, the $\mathbf{k}$-point
convergence is even slower and the difference between the point integration
and tetrahedron integration is thus even more pronounced for Ag. By increasing
the $\mathbf{k}$-point sampling, the point integration results will eventually
approach the results obtained with the tetrahedron method, but at a much
higher computational cost.

Figure 6: Real and imaginary parts of the dielectric function for Ag (solid)
and Si (dashed) using the tetrahedron and point integration methods.

####  VI.2.1 Screening in low-dimensional systems

In a three dimensional (3D) bulk crystal, the macroscopic dielectric function
is related to the macroscopic polarizability $\alpha_{M}$ of the material as

$\displaystyle\varepsilon_{M}(\mathbf{q},\omega)=1+4\pi\alpha_{M}(\mathbf{q},\omega).$
(59)

Since $\varepsilon_{M}$ is related to the macroscopic average of the induced
potential, it will depend on the unit cell in low-dimensional systems and tend
to unity when increasing the cell in the non-periodic direction. In contrast,
it is straighforward to generalize the definition of the polarizability to
such that it measures the induced dipole moment per length or area rather than
volume. The $d$-dimensional polarizability is thus defined as
$\alpha_{M}^{d}=\Omega^{d}\alpha_{M}$ where $\Omega^{d}$ is the cell volume
with the periodic directions projected out. For example, for a 2D material
$\Omega^{d}$ is simply the length of the unit cell in the non-periodic
direction. To improve convergence with respect to the size of the unit cell,
the Coulomb kernel is truncated using the reciprocal-space method introduced
in Ref. [206]. This enables efficient calculations of dielectric properties of
low-dimensional materials materials [207, 208]. This does, however, imply that
the polarizability cannot be evaluated directly from the dielectric constant
(which is just one), but one may obtain it directly from the susceptibility as

$\displaystyle\alpha_{M}=-\lim\limits_{\mathbf{q}\to
0}\frac{\Omega^{d}}{q^{2}}\chi_{00}(\mathbf{q}).$ (60)

###  VI.3 Adiabatic-connection fluctuation-dissipation theorem

The adiabatic-connection fluctuation-dissipation theorem (ACFDT) is a highly
promising approach for constructing accurate correlation functionals with non-
local characteristics. Unlike regular XC-functionals, the ACFDT correlation
functionals do not rely on error cancellation between exchange and
correlation. Instead, the ACFDT provides an exact theory for the electronic
correlation energy $E_{\text{c}}$ in terms of the interacting electronic
suseptibility, which can be combined with the exact exchange:

$E_{\text{c}}=-\int_{0}^{\infty}\frac{\mathrm{d}\omega}{2\pi}\int_{0}^{1}\mathrm{d}\lambda\mathrm{Tr}\left[v_{\text{c}}\chi^{\lambda}\left(w\right)-v_{\text{c}}\chi^{\text{KS}}(w)\right].$
(61)

The response function can be expressed in terms of the Kohn–Sham response
function and the exchange–correlation kernel $f_{\text{xc}}$

$\chi^{\lambda}(\omega)=\chi^{\text{KS}}(\omega)+\chi^{\text{KS}}(\omega)\left[\lambda
v_{c}+f_{\text{xc}}(\omega)\right]\chi^{\lambda}(\omega).$ (62)

The random phase approximation (RPA) can be derived from ACFDT if the XC-
kernel $f_{\text{xc}}$ is neglected. RPA has the strength of capturing non-
local correlation effects and provides high accuracy across different bonding
types including van der Waals interaction [209, 210, 29, 211].

Simple exchange–correlation kernels can also be incorporated into the response
function, such as the adiabatic LDA (ALDA) kernel. However, the locality of
adiabatic kernels leads to divergent characteristics of the pair-correlation
function [212]. This issue can be overcome by a renormalization (r) scheme
[37], which is implemented in GPAW as rALDA. This class of renormalized
kernels provides a significantly better description of the short-range
correlations and hence also yield highly accurate total correlation energies
[213, 214, 215, 216].

###  VI.4 Magnetic response

The LR-TDDFT framework described in Sec. VI.1 can be generalized to include
spin degrees of freedom [43, 217]. Using the four-component density as the
basic variable, $n^{\mu}(\mathbf{r})$ (where $\mu\in\\{0,x,y,z\\}$), one can
define the four-component susceptibility tensor:

$\chi^{\mu\nu}(\mathbf{r},\mathbf{r}^{\prime},t-t^{\prime})=-\frac{i}{\hbar}\theta(t-t^{\prime})\langle\left[\hat{n}^{\mu}_{0}(\mathbf{r},t),\hat{n}^{\nu}_{0}(\mathbf{r}^{\prime},t^{\prime})\right]\rangle_{0}.$
(63)

In a similar fashion to Eq. (43), the many-body $\chi^{\mu\nu}$ can be
calculated from the corresponding susceptibility tensor of the Kohn–Sham
system, $\chi^{\mu\nu}_{\mathrm{KS}}$. In the most general case, the Dyson
equation for $\chi^{\mu\nu}$ is a matrix equation in the $\mu$ and $\nu$
indices, explicitly coupling the charge and spin degrees of freedom. However,
for collinear magnetic systems, the transverse components of the
susceptibility decouple from the rest in the absence of spin–orbit coupling.

####  VI.4.1 Transverse magnetic susceptibility

Taking the spins to be polarized along the $z$-direction, the transverse
magnetic susceptibility of collinear nonrelativistic systems may be expressed
in terms of the spin-raising and spin-lowering density operators
$\hat{n}^{+}(\mathbf{r})=(\hat{n}^{x}(\mathbf{r})+i\hat{n}^{y}(\mathbf{r}))/2=\hat{\psi}^{\dagger}_{\uparrow}(\mathbf{r})\hat{\psi}_{\downarrow}(\mathbf{r})$
and
$\hat{n}^{-}(\mathbf{r})=(\hat{n}^{x}(\mathbf{r})-i\hat{n}^{y}(\mathbf{r}))/2=\hat{\psi}^{\dagger}_{\downarrow}(\mathbf{r})\hat{\psi}_{\uparrow}(\mathbf{r})$.
In the ALDA, the Dyson equation for $\chi^{+-}$ then becomes a scalar one
[43],

$\displaystyle\chi^{+-}(\mathbf{r},\mathbf{r}^{\prime},\omega)=$
$\displaystyle\chi^{+-}_{\mathrm{KS}}(\mathbf{r},\mathbf{r}^{\prime},\omega)+\int\mathrm{d}\mathbf{r}_{1}$
$\displaystyle\times\chi^{+-}_{\mathrm{KS}}(\mathbf{r},\mathbf{r}_{1},\omega)f^{-+}_{\mathrm{LDA}}(\mathbf{r}_{1})\chi^{+-}(\mathbf{r}_{1},\mathbf{r}^{\prime},\omega),$
(64)

where the transverse XC-kernel is given by
$f^{-+}_{\mathrm{LDA}}=2W_{\mathrm{LDA}}(\mathbf{r})/n^{z}(\mathbf{r})$. The
susceptibilities $\chi^{+-}$ and $\chi^{+-}_{\mathrm{KS}}$ are themselves
defined via the Kubo formula (63), where the latter can be evaluated
explicitly in complete analogy to Eq. (42):

$\displaystyle\chi^{+-}_{\mathrm{KS}}(\mathbf{r},\mathbf{r}^{\prime},\omega)=$
$\displaystyle\lim_{\eta\rightarrow
0^{+}}\sum_{n\mathbf{k}}\sum_{m\mathbf{k}^{\prime}}\left(f_{n\mathbf{k}\uparrow}-f_{m\mathbf{k}^{\prime}\downarrow}\right)$
$\displaystyle\times\frac{\psi^{*}_{n\mathbf{k}\uparrow}(\mathbf{r})\psi_{m\mathbf{k}^{\prime}\downarrow}(\mathbf{r})\psi^{*}_{m\mathbf{k}^{\prime}\downarrow}(\mathbf{r}^{\prime})\psi_{n\mathbf{k}\uparrow}(\mathbf{r}^{\prime})}{\hbar\omega-(\epsilon_{m\mathbf{k}^{\prime}\downarrow}-\epsilon_{n\mathbf{k}\uparrow})+i\hbar\eta}.$
(65)

In GPAW, the transverse magnetic susceptibility is calculated using a plane-
wave basis as described in Sec. VI.1.1, with the notable exception that no
special care is needed to treat the optical limit, since the Hartree kernel
plays no role in the Dyson equation (64). In terms of the temporal
representation, it is at the time of this writing only possible to do a
literal evaluation of Eq. (65) at the frequencies of interest. For metals,
this means that $\eta$ has to be left as a finite broadening parameter, which
has to be carefully chosen depending on the $\mathbf{k}$-point sampling, see
[43].

####  VI.4.2 The spectrum of transverse magnetic excitations

Based on the transverse magnetic susceptibility, one may calculate the
corresponding spectral function [43]

$\displaystyle S^{+-}(\mathbf{r},\mathbf{r}^{\prime},\omega)$
$\displaystyle=-\frac{1}{2\pi
i}\left[\chi^{+-}(\mathbf{r},\mathbf{r}^{\prime},\omega)-\chi^{-+}(\mathbf{r}^{\prime},\mathbf{r},-\omega)\right],$
(66)

which directly governs the energy dissipation in a collinear magnet relating
to induced changes in the spin projection along the $z$-axis $S_{z}$. In
particular, one can decompose the spectrum into contributions from majority
and minority spin excitations,
$S^{+-}(\mathbf{r},\mathbf{r}^{\prime},\omega)=A^{+-}(\mathbf{r},\mathbf{r}^{\prime},\omega)-A^{-+}(\mathbf{r}^{\prime},\mathbf{r},-\omega)$,
that is, into spectral functions for the excited states where $S_{z}$ has been
lowered or raised by one unit of spin angular momentum respectively:

$\displaystyle A^{+-}(\mathbf{r},\mathbf{r}^{\prime},\omega)=\sum_{\alpha>0}$
$\displaystyle\langle
0|\hat{n}^{+}(\mathbf{r})|\alpha\rangle\langle\alpha|\hat{n}^{-}(\mathbf{r}^{\prime})|0\rangle$
$\displaystyle\times\delta\left(\hbar\omega-\left[E_{\alpha}-E_{0}\right]\right).$
(67)

Here, $\alpha$ iterates the system eigenstates with $\alpha=0$ denoting the
ground state.

Figure 7:  Spectrum of transverse magnetic excitations for ferromagnetic hcp-
Co evaluated at $\mathbf{q}=5\mathbf{q}_{\mathrm{M}}/6$. The spectrum was
calculated using 8 empty-shell bands per atom, a plane-wave cutoff of 800 eV,
a $(60,60,36)$ $\mathbf{k}$-point mesh and a spectral broadening of $\eta=50$
meV. In the upper panel, the spectrum diagonal is depicted for the 1st and 2nd
Brillouin zone out of the hexagonal plane, from which the acoustic and optical
magnon frequencies can be respectively extracted. In the lower panel, the
spectrum of majority excitations is shown. The full spectral weight
$A(\omega)$ is calculated as the sum of all positive eigenvalues of $S^{+-}$,
the acoustic and optical mode lineshapes $a_{0}(\omega)$ and $a_{1}(\omega)$
are obtained via the two largest eigenvalues (which are significantly larger
than the rest) and the Stoner spectrum is extracted as the difference
$A_{\mathrm{S}}(\omega)=A(\omega)-a_{0}(\omega)-a_{1}(\omega)$.

For collinear magnetic systems, the spectrum $S^{+-}$ is composed of two types
of excitations: collective spin-wave excitations (referred to as magnons) and
excitations in the Stoner-pair continuum (electron–hole pairs of opposite
spin). Since GPAW employs a plane-wave representation of the spectrum,
$S^{+-}_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)$, one can directly
compare the calculational output to the inelastic neutron scattering cross-
section measured in experiments [218]. In particular, one can extract the
magnon dispersion directly by identifying the position of the peaks in the
spectrum diagonal, see Fig. 7. In this way, GPAW allows the user to study
various magnon phenomena in magnetic systems of arbitrary collinear order,
such as nonanalytic dispersion effects in itinerant ferromagnets [43],
correlation-driven magnetic phase transitions [219] and emergence of distinct
collective modes inside the Stoner continuum of an antiferromagnet [220].

Additionally, one can analyze the spectrum more intricately by extracting the
majority and minority eigenmodes from $S^{+-}$ as the positive and negative
eigenvalues respectively. Contrary to analysis of the plane-wave diagonal,
this makes it possible to completely separate the analysis of each individual
magnon lineshape as well as the many-body Stoner continuum, see Fig. 7.

####  VI.4.3 Liechtenstein MFT

Not only can the transverse magnetic susceptibility be used to study magnetic
excitations in a literal sense, but one can also use it to map the spin-
degrees of freedom to a classical Heisenberg model,

$E_{\mathrm{H}}=-\frac{1}{2}\sum_{i,j}\sum_{a,b}J_{ij}^{ab}\,\mathbf{u}_{ia}\cdot\mathbf{u}_{jb},$
(68)

where $i,j$ and $a,b$ denote the indices of the Bravais lattice and the
magnetic sublattice, respectively, $\mathbf{u}_{ia}$ being the direction of
spin-polarization of the given magnetic site. Based on the magnetic force
theorem (MFT), the LSDA Heisenberg exchange parameters $J_{ij}^{ab}$ can be
calculated from the reactive part of the static Kohn–Sham susceptibility [42]
and the effective magnetic field $B_{\mathrm{xc}}(\mathbf{r})=\delta
E_{\mathrm{xc}}/\delta m(\mathbf{r})$, using the well-known Liechtenstein MFT
formula [187]:

$J_{ij}^{ab}=-2\int_{\Omega_{ia}}\mathrm{d}\mathbf{r}\int_{\Omega_{jb}}\mathrm{d}\mathbf{r}^{\prime}\,B_{\mathrm{xc}}(\mathbf{r})\chi^{\prime+-}_{\mathrm{KS}}(\mathbf{r},\mathbf{r}^{\prime})B_{\mathrm{xc}}(\mathbf{r}^{\prime}).$
(69)

Here $\Omega_{ia}$ denotes the site volume, which effectively defines the
Heisenberg model (68).

Using GPAW’s plane-wave representation of $\chi^{\prime+-}_{\mathrm{KS}}$, one
can directly compute the lattice Fourier-transformed exchange parameters [42],

$\displaystyle\bar{J}^{ab}(\mathbf{q})$
$\displaystyle=\sum_{i}J_{0i}^{ab}e^{i\mathbf{q}\cdot\mathbf{R}_{i}}$ (70a)
$\displaystyle=-\frac{2}{\Omega_{\mathrm{cell}}}B_{\mathrm{xc}}^{\dagger}K^{a\dagger}(\mathbf{q})\chi^{\prime+-}_{\mathrm{KS}}(\mathbf{q})K^{b}(\mathbf{q})B_{\mathrm{xc}},$
(70b)

where the right-hand side of the second equality is written in a plane-wave
basis and $K_{a}$ denotes the sublattice site-kernel:

$K^{a}_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q})=\frac{1}{\Omega_{\mathrm{cell}}}\int_{\Omega_{0a}}\mathrm{}\mathrm{d}\mathbf{r}\,e^{-i(\mathbf{G}-\mathbf{G}^{\prime}+\mathbf{q})\cdot\mathbf{r}}.$
(71)

Since an a priori definition for the magnetic site volumes does not exist,
GPAW supplies functionality to calculate exchange parameters based on
spherical, cylindrical and/or parallelepipedic site configurations of variable
size.

Figure 8:  Magnon spectrum of ferromagnetic hcp-Co calculated using ALDA LR-
TDDFT (shown as a heat map), compared to the spin-wave dispersion of
Liechtenstein MFT. The acoustic magnon mode $a_{0}(\omega)$ is shown to the
left of the A-point, while the optical magnon mode $a_{1}(\omega)$ is shown to
the right. Please note that the two modes are degenerate in both the
$\mathrm{A}$ and $\mathrm{K}$ high-symmetry points. The calculations were
carried out using 8 empty-shell bands per atom, a plane-wave cutoff of 800 eV
and a $(36,36,24)$ $\mathbf{k}$-point mesh. A finite value of $\eta=100$ meV
was used to broaden the ALDA spectrum. For the MFT calculations, the
dispersion was calculated using linear spin-wave theory based on a Heisenberg
model of closed-packed spherical sites centered at each of the Co atoms.

Upon calculation of the exchange parameters $\bar{J}^{ab}(\mathbf{q})$, it is
straight-forward to compute the magnon dispersion within the classical
Heisenberg model using linear spin-wave theory and to estimate thermal
quantities such as the Curie temperature, see e.g. [42]. In Fig. 8 the MFT
magnon dispersion of hcp-Co is compared to majority magnon spectrum calculated
within LR-TDDFT. For Co, the two are in excellent agreement except for the
dispersion of the optical magnon mode along the $\mathrm{K}-\mathrm{M}-\Gamma$
high-symmetry path, where MFT underestimates the magnon frequency and neglects
the fine structure of the spectrum. The fine structure of the ALDA spectrum
appears due to the overlap between the magnon mode and the Stoner continuum.
This gives rise to so-called Kohn anomalies (nonanalytical points in the
magnon dispersion), which is a trademark of itinerant electron magnetism.
Since the itinerancy is largely ignored in a localized spin model such as
(68), one cannot generally expect to capture such effects.

###  VI.5 GW approximation

GPAW supports standard G0W0 quasiparticle (QP) calculations based on a first-
order perturbative treatment of the linearized QP equation [9]

$\displaystyle
E_{n\mathbf{k}\sigma}^{\text{QP}}=\varepsilon_{n\mathbf{k}\sigma}+Z_{n\mathbf{k}\sigma}\langle\psi_{n\mathbf{k}\sigma}|\Sigma_{\mathrm{GW}}(\varepsilon_{n\mathbf{k}\sigma})+v_{\text{x}}-v_{\text{xc}}|\psi_{n\mathbf{k}\sigma}\rangle,$
(72)

where $\varepsilon_{n\mathbf{k}\sigma}$ and $\psi_{n\mathbf{k}\sigma}$ are
Kohn–Sham eigenvalues and wave functions, and $v_{\text{xc}}$ and
$v_{\text{x}}$ are the local XC-potential and nonlocal exchange potential,
respectively. $\Sigma_{\mathrm{GW}}$ is the (dynamical part of) the GW self-
energy whose frequency-dependence is accounted to first order by the
renormalization factor
$Z_{n\mathbf{k}\sigma}=(1-\mathrm{Re}\>\Sigma_{\mathrm{GW}}^{\prime}(\varepsilon_{n\mathbf{k}\sigma}))^{-1}$.
As indicated by the $\sigma$-index, spin-polarized G0W0 calculations are
supported.

The GW self-energy is calculated in a plane-wave basis using full frequency
integration along the real axis to evaluate the convolution between $G$ and
$W$ [38]. Compared to alternative schemes employing contour deformation
techniques or analytical continuation [221, 222, 223], this approach is time-
consuming but numerically accurate and can provide the full spectral function.
A highly efficient and accurate evaluation of the self-energy based on a
multipole expansion of the screened interaction, $W$ [224], is currently being
implemented, and a GPU version of the full GW code is under development.

An important technical issue concerns the treatment of the head and wings of
$W$ ($\mathbf{G}=0$ and/or $\mathbf{G}^{\prime}=0$, respectively) in the
$\mathbf{q}\rightarrow 0$ limit. The divergence of the Coulomb interaction
appears both in the evaluation of the inverse dielectric matrix
$\varepsilon^{-1}_{\mathbf{GG}^{\prime}}$ and in the subsequent evaluation of
the screened interaction

$\displaystyle
W_{\mathbf{GG}^{\prime}}(\mathbf{q},\omega)=\sum_{\mathbf{G}^{\prime\prime}}\varepsilon^{-1}_{\mathbf{GG}^{\prime\prime}}(\mathbf{q},\omega)v_{c,\mathbf{G}^{\prime\prime}\mathbf{G}^{\prime}}(\mathbf{q}).$

For 3D bulk crystals, GPAW obtains these components by evaluating
$\varepsilon^{-1}_{\mathbf{GG}^{\prime}}$ on a dense $\mathbf{k}$-grid
centered at $\mathbf{k}=0$ while $v_{c,\mathbf{G}\mathbf{G}^{\prime}}$ can be
integrated numerically or analytically around the $\Gamma$-point.

For low-dimensional structures, in particular atomically thin 2D
semiconductors, GPAW can use a truncated Coulomb interaction to avoid
interactions between periodical images when evaluating $W$ [206]. It has been
shown that the use of a truncated Coulomb kernel leads to slower
$\mathbf{k}$-point convergence [207]. To mitigate this drawback, a special 2D
treatment of $W(q)$ at $q=0$, can be applied to significantly improve the
$\mathbf{k}$-point convergence [225]. A detailed account of the GW
implementation in GPAW can be found in Ref. [38].

Figure 9:  The real (Re) and imaginary (Im) parts of the self-energy matrix
elements at the Gamma point for valence (top) and conduction (bottom) bands
evaluated with Yambo and GPAW. Both codes are using full frequency integration
with a broadening of 0.1 eV. Yambo is using norm conserving pseudo potentials,
and GPAW its standard PAW setup. The k-point grid was $12\times 12\times 12$,
the plane wave cutoff was 200 eV, and the number of bands 200 for both codes.
The results are virtually indistinguishable.

Figure 9 shows two matrix elements of the dynamical GW self-energy for the
valence and conduction band states at the Gamma point. As can be seen, the
agreement with the corresponding quantities obtained with the Yambo GW code
[226] is striking.

###  VI.6 Bethe–Salpeter Equation (BSE)

In addition to the LR-TDDFT discussed in Sec. VI.1, the interacting response
function may be approximated by solving the Bethe–Salpeter equation (BSE)
[11]. In particular, for a certain wave vector $\mathbf{q}$, one may obtain
the two-particle excitations by diagonalizing the Hamiltonian

$\displaystyle H_{\begin{subarray}{c}\mathbf{k}_{2}m_{1}m_{2}\\\
\mathbf{k}_{4}m_{3}m_{4}\end{subarray}}(\mathbf{q})$
$\displaystyle=\delta_{m_{1}m_{3}}\delta_{m_{2}m_{4}}\delta_{\mathbf{k}_{2}\mathbf{k}_{4}}(\varepsilon_{m_{1}\mathbf{k}_{2}+\mathbf{q}}-\varepsilon_{m_{2}\mathbf{k}_{2}})$
(73)
$\displaystyle-(f_{m_{1}\mathbf{k}_{2}+\mathbf{q}}-f_{m_{2}\mathbf{k}_{2}})K_{m_{2}\mathbf{k}_{2}m_{1}\mathbf{k}_{2}+\mathbf{q}m_{3}\mathbf{k}_{4}+\mathbf{q}m_{4}\mathbf{k}_{4}},$

where $\varepsilon_{\mathbf{k}m}$ are the Kohn–Sham eigenvalues and
$f_{\mathbf{k}m}$ the associated occupation numbers. The kernel is defined by
$K=v_{\mathrm{c}}-W$ with

$\displaystyle v_{c}=$ $\displaystyle\langle
m_{2}\mathbf{k}_{2},m_{3}\mathbf{k}_{4}+\mathbf{q}|\hat{v}_{\mathrm{c}}|m_{1}\mathbf{k}_{2}+\mathbf{q},m_{4}\mathbf{k}_{4}\rangle$
(74) $\displaystyle W=$ $\displaystyle\langle
m_{2}\mathbf{k}_{2},m_{3}\mathbf{k}_{4}+\mathbf{q}|\hat{W}|m_{4}\mathbf{k}_{4},m_{1}\mathbf{k}_{2}+\mathbf{q}\rangle,$
(75)

where $\hat{W}=\varepsilon^{-1}v_{\mathrm{c}}$ is the static screened Coulomb
interaction. The matrix elements of the kernel are evaluated in a plane-wave
basis where they are easily expressed in terms of the pair densities (48) and
the reciprocal-space representation of the dielectric matrix (54).

In the Tamm–Dancoff approximation, one only includes states with
$\varepsilon_{m_{1}\mathbf{k}_{2}+\mathbf{q}}-\varepsilon_{m_{2}\mathbf{k}_{2}}>0$
in Eq. (73) and the BSE Hamiltonian becomes Hermitian [11]. The interacting
retarded response function may then be written as

$\displaystyle\chi^{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)=\frac{1}{\Omega}\sum_{\lambda}\frac{B_{\lambda}(\mathbf{q},\mathbf{G})C_{\lambda}^{*}(\mathbf{q},\mathbf{G}^{\prime})}{\omega-
E_{\lambda}(\mathbf{q})+i\eta},$ (76)

where

$\displaystyle
B_{\lambda}(\mathbf{q},\mathbf{G})=\sum_{\mathbf{k}m_{1}m_{2}}n_{m_{2}\mathbf{k}m_{1}\mathbf{k}+\mathbf{q}}(\mathbf{G}+\mathbf{q})A^{\lambda}_{m_{1}m_{2}\mathbf{k}}(\mathbf{q}),$
(77) $\displaystyle
C_{\lambda}(\mathbf{q},\mathbf{G})=\sum_{\mathbf{k}m_{1}m_{2}}n_{m_{2}\mathbf{k}m_{1}\mathbf{k}+\mathbf{q}}(\mathbf{G}+\mathbf{q})A^{\lambda}_{m_{1}m_{2}\mathbf{k}}(\mathbf{q})$
$\displaystyle\qquad\qquad\qquad\qquad\times(f_{m_{1}\mathbf{k+q}}-f_{m_{2}\mathbf{k}}),$
(78)

and $E_{\lambda}(\mathbf{q})$ denotes the eigenvalue of the Hamiltonian (73)
corresponding to the eigenvector
$A^{\lambda}_{m_{1}m_{2}\mathbf{k}}(\mathbf{q})$.

In GPAW, the construction of the BSE Hamiltonian (73) proceeds in two steps
[39]. First, the static screened interaction is calculated at all inequivalent
$\mathbf{q}$-points in a plane-wave basis and the Kernel is then subsequently
expressed in a basis of two-particle KS states. The first step is efficiently
parallelized over either states or $\mathbf{k}$-points and the second step is
parallelized over pair densities. The Hamiltonian elements are thus
distributed over all CPUs and the diagonalization is carried out using
ScaLAPACK such that the full Hamiltonian is never collected on a single CPU.
The dimension of the BSE Hamiltonian and memory requirements are, therefore,
only limited by the number of CPUs used for the calculation. We note that the
implementation is not limited to the Tamm–Dancoff approximation, but
calculations become more demanding without it. The response function may be
calculated for spin-paired as well as spin-polarized systems and spin–orbit
coupling can be included non-selfconsistently [227, 45]. In low-dimensional
systems it is important to eliminate the spurious screening from periodic
images of the structure, which is accomplished with the truncated Coulomb
interaction of Ref. [206].

The most important application of BSE is arguably the calculation of optical
absorption spectra of solids where BSE provides an accurate account of the
excitonic effects that are not captured by semi-local approximations for
$K_{\mathrm{Hxc}}$ (44). In 2D systems, the excitonic effects are particularly
pronounced due to inefficient screening [207, 227] and in Fig. 10 we show the
2D polarizability of WS2 calculated from BSE with $\mathbf{q}=\mathbf{0}$.
Comparing with Fig. 3, it is observed that the absorption edge is expected to
be located at the $K$ point, where spin–orbit coupling splits the highest
valence band by 0.45 eV. This splitting is seen as two excitonic peaks below
the band gap, which is interpreted as distinct excitons originating from the
highest and next-highest valence bands (the splitting of the lowest conduction
band is negligible in this regard). For comparison, we also show the RPA
polarizabilty obtained with the BSE module by neglecting the screened
interaction in the kernel and this shows the expected absorption edge at the
band gap. This yields identical results to the Dyson equation approach of Sec.
VI.1, but has the advantage that the eigenvalues and weights of Eq. (76) are
obtained directly such that the artificial broadening $\eta$ may be varied
without additional computational cost. The eigenstate decomposition also
allows one to access “dark states” and the BSE calculation reveals two (one
for each valley) triplet-like excitons that are situated 70 meV below the
lowest bright exciton in Fig. 10.

Figure 10: 2D polarizability of WS2 calculated from the BSE and the RPA. For
this calculation, we included spin–orbit coupling and used the 2D Coulomb
truncation to eliminate screening from periodic images. A gamma-centered
uniform $\mathbf{k}$-point grid of $48\times 48$ was applied and 8 valence
states and 8 conduction states (shifted by 1 eV to match the GW band gap [45])
was included in the Tamm–Dancoff approximation. This yields a BSE Hamiltonian
of size $N\times N$ with $N=147456$, which is easily diagonalized with
ScaLAPACK on 240 CPUs.

In addition to optical properties, GPAW allows for solving the BSE at finite
$\mathbf{q}$, which can used to obtain plasmon dispersion relations from the
EELS (57) and magnon dispersions from the transverse magnetic susceptibility
of Sec. VI.4 [228].

###  VI.7 Electron–phonon coupling

The electron-phonon coupling is the origin of several important materials
properties, ranging from electrical and thermal conductivity to
superconductivity. In addition it provides access to the deformation
potential, which can be used to obtain transport properties for electrons and
holes in solids [229].

The first-order electron–phonon coupling matrix
$g_{mn}^{\nu}(\mathbf{k},\mathbf{q})$ measures the strength of the coupling
between a phonon branch $\nu$ with wave vector $\mathbf{q}$ and frequency
$\omega_{\nu}$ and the electronic states $m(\mathbf{k}+\mathbf{q})$ and
$n(\mathbf{k})$ [230, 231]:

$g_{mn}^{\nu}(\mathbf{k},\mathbf{q})=\sqrt{\frac{\hbar}{2m_{0}\omega_{\nu}}}M_{mn}^{\nu}(\mathbf{k},\mathbf{q})$
(79)

with

$M_{mn}^{\nu}(\mathbf{k},\mathbf{q})=\langle\psi_{m\mathbf{k}+\mathbf{q}}|\nabla_{u}v^{\text{KS}}\cdot\mathbf{e}_{\nu}|\psi_{n\mathbf{k}}\rangle.$
(80)

Here $m_{0}$ is the sum of the masses of all the atoms in the unit cell,
$\nabla_{u}$ denotes the gradient with respect to atomic displacements and
$\mathbf{e}_{\nu}$ projects the gradient onto the direction of the phonon
displacements. In the case of the three translational modes at
$|\mathbf{q}|=0$, the matrix elements $g_{mn}^{\nu}$ vanishes as a consequence
of the acoustic sum rule [231].

In GPAW, $\nabla_{u}v^{\text{KS}}(\mathbf{r})$ is determined using a finite-
difference method with a super cell description of the system. This step can
be performed in any of the wave function representations available in GPAW.
The derivative is then projected onto a set of atomic orbitals from an LCAO
basis $\phi_{NM}$, where $N$ denotes the cell index and $M$ the orbital index:

$\displaystyle\mathbf{g}_{\begin{subarray}{c}NM\\\
N^{\prime}M^{\prime}\end{subarray}}^{sc}=FT\left[\langle\phi_{NM}(\mathbf{k})|\nabla_{u}v^{\text{KS}}|\phi_{N^{\prime}M^{\prime}}(\mathbf{k})\rangle\right].$
(81)

The Fourier transform from the Bloch to the real space representation makes it
possible to compute $M_{mn}^{\nu}$ for arbitrary $\mathbf{q}$. Finally, the
electron–phonon coupling matrix is obtained by projecting the matrix
corresponding to the supercell into the primitive unit cell bands $m,n$ and
phonon modes $\nu$:

$M_{mn}^{\nu}(\mathbf{k},\mathbf{q})=\sum_{\begin{subarray}{c}NM\\\
N^{\prime}M^{\prime}\end{subarray}}C_{mM}^{*}C_{nM^{\prime}}\mathbf{g}_{\begin{subarray}{c}NM\\\