#  koopmans: an open-source package for accurately and efficiently predicting
spectral properties with Koopmans functionals

Edward Linscott<EMAIL_ADDRESS>[  Nicola Colonna  Laboratory for
Neutron Scattering and Imaging, Paul Scherrer Institut, 5232 Villigen,
Switzerland  Riccardo De Gennaro  [  Ngoc Linh Nguyen  Faculty of Materials
Science and Engineering, Phenikaa University, Hanoi 12116, Vietnam  Giovanni
Borghi  [  Andrea Ferretti  Centro S3, CNR–Istituto Nanoscienze, 41125 Modena,
Italy  Ismaila Dabo  Department of Materials Science and Engineering,
Materials Research Institute, and Institutes of Energy and the Environment,
The Pennsylvania State University, University Park, Pennsylvania 16802, USA
Nicola Marzari<EMAIL_ADDRESS>[

###### Abstract

Over the past decade we have developed Koopmans functionals, a computationally
efficient approach for predicting spectral properties with an orbital-density-
dependent functional framework. These functionals impose a generalized
piecewise linearity condition to the entire electronic manifold, ensuring that
orbital energies match the corresponding electron removal/addition energy
differences (in contrast to semi-local DFT, where a mismatch between the two
lies at the heart of the band gap problem and, more generally, the
unreliability of Kohn-Sham orbital energies). This strategy has proven to be
very powerful, yielding molecular orbital energies and solid-state band
structures with comparable accuracy to many-body perturbation theory but at
greatly reduced computational cost while preserving a functional formulation.
This paper reviews the theory of Koopmans functionals, discusses the
algorithms necessary for their implementation, and introduces koopmans, an
open-source package that contains all of the code and workflows needed to
perform Koopmans functional calculations and obtain reliable spectral
properties of molecules and materials.

THEOS] Theory and Simulation of Materials (THEOS), École Polytechnique
Fédérale de Lausanne, 1015 Lausanne, Switzerland Contributed equally to this
work Contributed equally to this work [MARVEL] National Centre for
Computational Design and Discovery of Novel Materials (MARVEL), École
Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland THEOS] Theory
and Simulation of Materials (THEOS), École Polytechnique Fédérale de Lausanne,
1015 Lausanne, Switzerland Phenikaa Research and Technology Institute (PRATI),
A&A Green Phoenix Group JSC, No. 167 Hoang Ngan, Trung Hoa, Cau Giay, Hanoi
11313, Vietnam THEOS] Theory and Simulation of Materials (THEOS), École
Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland Now at Liceo
Manfredo Fanti, 41012 Carpi, Italy THEOS] Theory and Simulation of Materials
(THEOS), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland
[MARVEL] National Centre for Computational Design and Discovery of Novel
Materials (MARVEL), École Polytechnique Fédérale de Lausanne, 1015 Lausanne,
Switzerland Laboratory for Materials Simulations, Paul Scherrer Institut, 5232
Villigen, Switzerland

##  1 Introduction

How can one accurately and efficiently predict spectral properties of
molecules and materials _ab initio_? Currently, the most accurate and popular
approaches to compute charged excitation energies are Green’s functions
methods such as many-body perturbation theory (GW) 1, 2 or wavefunction
methods such as quantum Monte Carlo 3 and equation-of-motion coupled cluster 4
— although for the latter, calculations for the solid state (rather than for
molecules) are far from routine. Of these approaches, GW is computationally
the least expensive, scaling as $\mathcal{O}(N^{4})$, where $N$ is the number
of electrons in the system.

Despite ongoing progress in the field of GW5, performing these calculations is
not straightforward. The aforementioned scaling of $\mathcal{O}(N^{4})$ can
still be an obstacle, and the calculations themselves can be challenging: they
converge slowly with respect to the number of empty states included (which
increases the importance of constructing transferable pseudopotentials that
avoid ghost states6), and there is a strong interdependence of the results on
different calculation parameters, which makes achieving convergence
challenging at best. This hampers routine applications of GW (especially in a
high-throughput context, where the calculations must be unsupervised) 7.
Finally, while in principle GW and many-body perturbation theory are
systematically improvable — that is to say, by increasing the number of
diagrams included in the calculations, the results should progressively
converge to the correct answer (with GW outperforming GW0 in turn
outperforming G0W0) — in practice this does not appear to hold8.

Alternatively, one could try and calculate the energies of electronic
excitations with density-functional theory (DFT) 9, 10. DFT has proven to be a
remarkably successful theory for predicting the ground-state properties of
solids, surfaces, nanoparticles, and molecules 11, 12. It is typically
inexpensive, and these days such calculations are generally robust and can be
treated as a “black box”. However, DFT is a theory of total energies, and
while the Kohn-Sham auxiliary system is a powerful construct, the Kohn-Sham
eigenvalues are not necessarily related to the energies of charged excitations
(with the exception of the highest occupied molecular orbital, or HOMO, which
is related to the exponential decay of the density 13). Nevertheless, these
eigenvalues can bear qualitative or even quantitative resemblance to
experimental quasiparticle energies, and it is common practice to interpret
them as such, motivated by the fact that the Kohn-Sham potential is the best
local and static approximation to the electronic self-energy 14. Asides from
this theoretical disconnect, problems also arise from the additional
approximations inherent in exchange-correlation functionals. In the case of
local and semi-local functionals, a key qualitative failure arises from the
erroneous convex curvature in the total energy as a function of the total
number of electrons in the system, which should instead be piecewise-linear
15. This curvature explains in part the disagreement between first ionization
potentials as calculated via total energy differences compared to Kohn-Sham
eigenvalues.

Many strategies have emerged that attempt to restore the piecewise linearity
of the energy functional — the hope being that the resulting Kohn-Sham
eigenvalues will yield accurate excitation energies. For example, DFT+_U_
imposes a penalty functional to a localized subspace that restores linearity
in the energy with respect to the occupation of this subspace16, 17.
Similarly, hybrid functionals mix semi-local functionals with Hartree-Fock
exchange (which happens to exhibit a concave curvature), which means that for
a specific mixing fraction of the two functionals there will be an overall
error cancellation 18, 19, 20, 21. Recent state-of-the-art approaches that
employ curvature corrections to yield reliable quasiparticle energies include
screened, range-separated, and dielectric-dependent hybrid functionals with
tuned mixing or range-separation parameters22, 23, 24, 25, 26, 27, as well as
the Koopmans-Wannier method of Wang and coworkers 28 and the localized orbital
scaling correction (LOSC) of Yang and coworkers 29, 30, 31. Piecewise
linearity is also central to ensemble density functional theory 32, 33. Even
DM21, the recent machine-learned exchange-correlation functional created by
Google DeepMind, was constructed around the idea of restoring piecewise
linearity 34.

Starting in 2009, we have introduced and developed the concept of Koopmans
functionals 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47. By imposing a
generalized piecewise linearity condition and relating quasiparticle energies
to total energy differences, these functionals address the above issues, and
as a consequence they yield spectroscopic properties (such as molecular
ionization potentials, electron affinities, solid-state band structures, and
band-edge alignments) with comparable accuracy to state-of-the-art GW
approaches, but at greatly reduced computational cost while preserving a
functional formulation. This has all been implemented in koopmans, an open-
source package that allows non-experts to perform their own Koopmans
functional calculations, and which is built upon the popular Quantum ESPRESSO
distribution. This paper provides an overview of the theory of Koopmans
functionals (Section 2), describes the algorithms that enable their
implementation in the koopmans package (Section 3), and demonstrates how these
tools can be deployed to predict spectral properties using the examples of
ozone, silicon, and zinc oxide (Section 4).

##  2 Koopmans functionals

###  2.1 Fundamental concepts

For a spectral theory, the orbital energies $\varepsilon_{i}$ should match the
total energy differences corresponding to electron removal $E(N)-E_{i}(N-1)$
and addition $E_{i}(N+1)-E(N)$. This is trivially true for the exact Green’s
function, whose poles correspond directly to these total energy differences,
but there is no such connection in Kohn-Sham DFT. The only exception to this
is the HOMO, but even there the violation of piecewise linearity in density
functional approximations leads to a mismatch between the HOMO eigenvalue with
the corresponding total energy difference (i.e. the negative of the ionization
potential).

Koopmans functionals restore this correspondence, by imposing the condition
that the orbital energies
$\varepsilon_{i}=\left\langle\varphi_{i}\right|\hat{H}\left|\varphi_{i}\right\rangle=\frac{dE}{df_{i}}$
of orbitals $\varphi_{i}$ should be independent of that orbital’s occupation
$f_{i}$:

$\varepsilon_{i}=\text{constant with respect to}\ f_{i}$ (1)

It follows from Janak’s theorem that this is equivalent to a “generalized”
piecewise linearity condition where the total energy of the system is
piecewise linear with respect to the change of occupation of _any_ orbital.
This is a sufficient but not a necessary condition to fulfil the much more
well-known piecewise linearity condition48, which states that the total energy
of the system is piecewise linear with respect to its _total_ number of
electrons. In passing, we mention that eq. 1 is reminiscent of a photoemission
experiment, where an electron is removed from a Dyson orbital.

Imposing this condition will require a beyond-DFT approach, and is not simply
a matter of correcting density functional approximations within a DFT
framework. We can see that this must be the case by considering the exact
density functional, for which the Koopmans corrections must be non-vanishing.
(This is because while the negative of the HOMO energy for the exact density
functional matches the ionization potential, there is no such guarantee for
the other eigenenergies13.)

The generalized piecewise linearity condition of eq. 1 is imposed on a “base”
functional (here, approximate or exact DFT) by removing, orbital-by-orbital,
the non-linear dependence of the energy $E$ on the orbital occupation $f_{i}$
and replacing it with a term that is linear in $f_{i}$:

$E^{\text{Koopmans}}={E^{\text{DFT}}}+\sum_{i}\left[-\left(E^{\mathrm{DFT}}-\left.E^{\mathrm{DFT}}\right|_{f_{i}=0}\right)+f_{i}\eta_{i}\right]$
(2)

where $\left.E^{\mathrm{DFT}}\right|_{f_{i}=f}$ corresponds to the DFT energy
of the $(N-1+f)$-electron system, with the occupancy of orbital $i$
constrained to be $f$. The first term in the square brackets removes the
dependence of the total energy on $f$, and the second term replaces it with a
term explicitly linear in $f$. This construction is reminiscent of the SIC
functional of Ref. 49, but here the correction is generalized to the entire
electronic manifold.

Here, one must choose a suitable slope $\eta_{i}$ for this linear term; one
option is to use the energy difference between fully-occupied and empty
orbitals

$\eta^{\mathrm{KI}}_{i}=\left.E^{\mathrm{DFT}}\right|_{f_{i}=1}-\left.E^{\mathrm{DFT}}\right|_{f_{i}=0}$
(3)

giving rise to the Koopmans integer (KI) functional. Note that this
formulation provides Koopmans functionals with meaningful eigenvalues, because
they now correspond to total energy differences, which in the scope of DFT are
formally meaningful and much more reliable than Kohn-Sham eigenvalues. It can
be seen from eqs. 2 and 3 that the KI functional gives, at integer
occupations, the same total energy as the base functional, but has different
derivatives and hence yields different spectral properties. (This will be
discussed further in Section 2.4.1).

Equations 2 and 3 are difficult to evaluate unless we only consider the
explicit dependence of the DFT energy on the orbital occupancies, neglecting
the implicit dependence of the orbitals $\varphi_{i}(\mathbf{r})$ on their own
occupation $f_{i}$, in which case

$\left.E^{\mathrm{DFT}}\right|_{f_{i}=f}=E^{\mathrm{DFT}}[\rho-\rho_{i}+fn_{i}]$
(4)

where $n_{i}(\mathbf{r})=|\varphi_{i}(\mathbf{r})|^{2}$ is the density of
orbital $i$ and
$\rho_{i}(\mathbf{r})=f_{i}|\varphi_{i}(\mathbf{r})|^{2}=f_{i}n_{i}(\mathbf{r})$
is the occupancy-weighted density of orbital $i$. Orbital relaxation — or,
equivalently, screening — is instead accounted for _post hoc_ by scaling the
unscreened correction by a scalar coefficient $\alpha_{i}$. Crucially, these
coefficients can be calculated _ab initio_ at the level of DFT via linear
response or total energy differences 43, 44. This brings us, finally, to the
Koopmans energy functional:

$E^{\text{Koopmans}}[\rho,{\\{\rho_{i}\\}}]={E^{\text{DFT}}[\rho]}+\sum_{i}\alpha_{i}\left[-\left(E^{\mathrm{DFT}}[\rho]-E^{\mathrm{DFT}}[\rho-\rho_{i}]\right)+f_{i}\eta_{i}\right]$
(5)

Figure 1: The total energy $E$ of a CH4 molecule as a function of the
occupancy of the $1a_{1}$ molecular orbital (crosses) and one of the $1t_{2}$
molecular orbitals (pluses). The absolute total energy is shown in (a) while
the deviation from piecewise linearity is shown in (b). For both orbitals, PBE
gives a total energy that is erroneously convex, while the KI correction
successfully linearizes the total energy. Screening is key to this success; in
its absence, the KI correction overcorrects the PBE base functional and yields
a concave energy curve. Note that each orbital $\varphi_{i}$ is obtained from
the charge-neutral system ($f_{i}=1$) and is frozen throughout (while all
others are relaxed). If that orbital was not frozen, then as $f_{i}\rightarrow
0$ the orbital would always morph into the LUMO of the $N-1$-electron system
and both sets of curves would be identical.

In Figure 1 we show the efficacy of this linearizing correction when applied
to two orbitals in methane. The full derivation of eq. 5 can be found in
Supporting Information S1. This functional is actually very different from
semi-local DFT functionals; this will be elaborated upon in the following
sections.

###  2.2 Orbital-density dependence

The one important distinction that is worth making immediately is that
Koopmans functionals are not density functionals, but _orbital-density-
dependent_ (ODD) functionals. This is because they are dependent on the
individual orbital densities $\\{\rho_{i}\\}$ and not just the total
electronic density $\rho$. A direct consequence of this is that Koopmans
functionals — much like other ODD functionals such as the Perdew-Zunger self-
interaction correction (PZSIC) — are no more invariant under unitary
transformation of the occupied manifold, and their minimization requires extra
care. The variation of $E^{\rm Koopmans}$ in eq. 5 with respect to an
arbitrary change of each orbital $\varphi_{i}$ (density $\rho_{i}$) leads to
the Euler-Lagrange equations

$h^{\rm DFT}|\varphi_{i}\rangle+v_{i}^{\rm
ODD}|\varphi_{i}\rangle=\sum_{j}\Lambda_{ji}|\varphi_{j}\rangle$ (6)

where $h^{\rm DFT}(\mathbf{r})=\frac{\delta E^{\rm
DFT}}{\delta\rho(\mathbf{r})}$ is the Hamiltonian of the underlying DFT energy
functional, $v_{i}^{\rm ODD}(\mathbf{r})$ is the orbital-density-dependent
potential associated with the orbital $\varphi_{i}$, and $\Lambda_{ji}$ is the
matrix of Lagrangian multipliers enforcing orthonormality constraints. Because
of the ODD contribution, within the space spanned by the orbitals
$\\{\varphi_{i}\\}$) the energy is representation-dependent and a proper
minimization of the functional requires its variation with respect to
infinitesimal unitary transformations among the occupied orbitals to vanish
50, 51, 40, leading to the Pederson condition 50

$\langle\varphi_{i}|h_{i}|\varphi_{j}\rangle=\langle\varphi_{i}|h_{j}|\varphi_{j}\rangle.$
(7)

The self-consistent solution of eqs. 6 and 7 define the proper minimum of the
Koopmans functionals, and the minimizing orbitals are known as the
_variational_ orbitals. The implementation of this minimization procedure will
be discussed later in Section 3.1.

At the minimum, as a consequence of eq. 7, the $\Lambda$ matrix becomes
Hermitian and can be diagonalized allowing us to define a set of _canonical_
orbitals and energies. This mirrors the definition of canonical orbitals and
energy in Hartree-Fock theory where, among all the equivalent sets of orbitals
(those related by a unitary transformations) that minimize the functional, the
canonical orbitals are recognized as those that also make the energy
functional stationary when a fraction of electron is added to or removed from
the system, thus qualifying these as electron addition/removal energies. This
also applies to ODD functionals, as discussed in detail in Ref. 52 for the
case of PZSIC. Moreover, canonical orbitals typically display the symmetry of
the Hamiltonian operator (e.g. are Bloch states in periodic systems 45 as
shown in Fig. 2(a)) and, in analogy to exact DFT, the energy of the highest
occupied canonical orbitals has been numerically shown to determine the
asymptotic decay of the ground-state charge density 53. For all these reasons,
the canonical orbitals and the corresponding eigenvalues are usually
interpreted as Dyson orbitals and quasiparticle energies. Nevertheless, it is
important to stress that the reliability of canonical energies (and their
correspondence with total energy differences) is not directly imposed by the
Koopmans correction, but instead is inherited via the variational orbitals.
That is to say: the Koopmans corrections are applied to the variational
orbitals, and thus the Koopmans functional is linear with respect to the
occupancy of variational orbitals. The canonical orbitals are composed of some
linear combination of variational orbitals, and their energies (i.e. the
quasiparticle energies) are subject to a weighted combination of corrective
potentials arising from their constituent variational orbitals.

Given their central role in the theory, it is important to discuss the key
features of variational orbitals. In contrast to canonical orbitals,
variational orbitals are typically very localized in space (see Fig. 2(b)). As
was recognized long ago50, eq. 7 is a localization condition that, once
satisfied, leads to orbitals that resemble Boys orbitals in molecules or,
equivalently, maximally localized Wannier functions in periodic systems54. The
localization of the variational orbitals is a common feature of ODD
functionals and a key property for Koopmans functionals, in particular when it
comes to dealing with periodic systems. By applying Koopmans corrections to a
set of localized orbitals, the corrections are well-defined and non-vanishing
for both small molecules, infinite bulk systems, and everything in between,
preserving size-consistency 43. Contrast this to if we were to apply the
corrections to the canonical orbitals, in which case they would become ill-
defined in the bulk limit. In order to understand why this is the case, it is
useful to return to the connection between the Koopmans construction and the
$\Delta$SCF approach. In a nutshell, the ultimate effect of the Koopmans
correction is to revert the wrong eigenvalue from the underlying (approximate)
density functional into a total energy difference ($\Delta$SCF) between the
neutral system and the system with plus or minus one electron evaluated using
the same density functional. This means that the success of the approach
relies on the quality of the $\Delta$SCF value at the approximate DFT level.
It is well known that evaluating this total energy difference when removing an
electron from a completely delocalized state reduces to the derivative of the
total energy with respect to the particle number 55, 56, 57, which, for a
local or semilocal density-functional approximation, is the negative of the
KS-DFT eigenvalue. This means that for a standard density functional in the
thermodynamic limit there is no difference between the $\Delta$SCF and the KS
eigenvalues and as a consequence the Koopmans corrections vanish. To overcome
this issue, two routes are possible: either improving the base functional in
such a way to have improved $\Delta$SCF energies in the most general case, or
retaining the simplicity of local and semi-local density-functionals and
working in a localized representation of the orbitals 58, 28. Indeed, the
total energy differences of approximate density functionals also become
accurate when computed on localized orbitals (e.g. typically, semi-local
$\Delta$SCF calculations accurately predict localized defect levels relative
to the average electrostatic potential 59). Thus, by applying the Koopmans
corrections to the variational orbitals (and not the canonical orbitals), the
Koopmans corrections are well-defined and non-vanishing also in the bulk
limit, and yield accurate band structures compared to experiment. See Ref. 43
for more details.

(a) a canonical orbital

(b) two variational orbitals

Figure 2: Example canonical and variational orbitals of polyethylene. Figure
adapted from Ref. 43.

Moving from a DFT framework to an ODDFT framework may appear like an
unnecessary complication. This is not the case: ODDFTs are a very natural way
to generalize a static functional theory like DFT to predict spectral
information. Ultimately, the spectral properties of a many-body electronic
system are exactly described by its non-local and dynamic self-energy. The
exact Kohn-Sham potential is the best local and approximation to this self-
energy 14. If we instead consider local but dynamic approximations, one enters
into the domain of spectral functional theories, where the exact spectral
functional predicts exactly the spectral density $\rho(\mathbf{r},\omega)$ 60.
ODDFTs can be interpreted as energy-discretized spectral functional theories
39, so as such an ODDFT framework is a sensible choice when attempting to
predict spectral properties.

###  2.3 Accounting for screening effects

As discussed earlier in Section 2.1, we account for orbital relaxation _post
hoc_ via screening parameters $\\{\alpha_{i}\\}$ and we can calculate these
parameters _ab initio_. But how?

The crucial point is that we would like the total energy to be piecewise
linear: that is, we would like orbital energies (specifically, the expectation
value of the Hamiltonian on a given variational orbital) to match the
corresponding total energy differences when adding/removing an electron from
this orbital, without the frozen-orbital assumption that we made earlier.
Specifically, we would like $\lambda_{ii}(\alpha,f)=\Delta
E^{\text{Koopmans}}_{i}$, where

$\lambda_{ii}(\alpha,f)=\left.\left\langle\varphi_{i}\right|\hat{h}^{\mathrm{DFT}}+\alpha\hat{v}^{\mathrm{Koopmans}}\left|\varphi_{i}\right\rangle\right|_{f_{i}=f}=\left.\frac{\partial
E^{\text{Koopmans}}}{\partial f_{i}}\right|_{f_{i}=f}$ (8)

is the expectation value of the Hamiltonian for a given variational orbital
$\varphi_{i}$, and

$\Delta
E^{\text{Koopmans}}_{i}=\begin{cases}E^{\text{Koopmans}}(N)-E^{\text{Koopmans}}_{i}(N-1)&\text{for
occupied orbitals}\\\
E^{\text{Koopmans}}_{i}(N+1)-E^{\text{Koopmans}}(N)&\text{for empty
orbitals}\end{cases}$ (9)

where $E^{\text{Koopmans}}_{i}(N\pm 1)$ is the total energy of the system
where we add/remove an electron from variational orbital $i$ and allow the
rest of the system to relax, with all the other orbitals remaining orthogonal
to $\left|\varphi_{i}\right\rangle$.

We use this condition to determine the screening parameters _ab initio_.
Specifically, given a starting guess $\\{\alpha^{0}_{i}\\}$ for the screening
parameters, an improved guess for the screening parameters can be obtained via

$\alpha^{n+1}_{i}=\alpha^{n}_{i}\frac{\Delta
E^{\text{Koopmans}}_{i}-\lambda_{ii}(0,1)}{\lambda_{ii}(\alpha^{n}_{i},1)-\lambda_{ii}(0,1)};\qquad\Delta
E^{\text{Koopmans}}_{i}=E^{\text{Koopmans}}(N)-E^{\text{Koopmans}}_{i}(N-1)$
(10)

for occupied orbitals and

$\alpha^{n+1}_{i}=\alpha^{n}_{i}\frac{\Delta
E^{\text{Koopmans}}_{i}-\lambda_{ii}(0,0)}{\lambda_{ii}(\alpha^{n}_{i},0)-\lambda_{ii}(0,0)};\qquad\Delta
E^{\text{Koopmans}}_{i}=E^{\text{Koopmans}}_{i}(N+1)-E^{\text{Koopmans}}(N)$
(11)

for empty orbitals, where $E^{\text{Koopmans}}_{i}(N\pm 1)$ is the total
energy of the $N\pm 1$ electron system where we take the $N$-electron system,
take this variational orbital $i$ and fill/empty it, and then hold it frozen
while the rest of the system is allowed to relax (while remaining orthogonal).
These equations yield the screening parameters that satisfy
$\lambda_{ii}(\alpha,f)=\Delta E_{i}^{\text{Koopmans}}$ if we assume a linear
dependence of $\lambda_{ii}$ on $\alpha_{i}$ and approximate the total energy
as a function of $f_{i}$ to second order. By iterating to self-consistency we
lift these approximations and guarantee that $\lambda_{ii}(\alpha,f)=\Delta
E_{i}^{\text{Koopmans}}$ is satisfied. Typically, only a few iterations are
required in order to reach self-consistency, especially if one starts from a
physically-motiviated initial guess (such as the static limit of the inverse
dielectric function $\varepsilon^{-1}$ in the case of bulk systems). All of
these ingredients for calculating $\alpha^{n+1}_{i}$ are obtained from
constrained Koopmans and DFT calculations. Specifically, a $N$-electron
Koopmans calculation yields $E^{\text{Koopmans}}(N)$ and
$\lambda_{ii}(\alpha,f)$ (for both $\alpha=\alpha_{i}^{n}$ and 0, and $f=1$
for filled orbitals and $0$ for empty). Meanwhile, a constrained $N\pm
1$-electron calculation yields $E^{\text{Koopmans}}_{i}(N\pm 1)$.

For a periodic system, this method for determining the screening parameters
requires a supercell treatment. This is because the $N\pm 1$-electron systems
contain a charged defect (because we have filled/emptied a localized orbital)
and a supercell is required in order to remove the spurious interactions
between periodic images 43, 45. Section 3.3 will discuss an efficient linear-
response reformulation of this problem that avoids a supercell treatment (and
can also be used for molecules).

###  2.4 Koopmans variants

As we saw previously in Section 2.1, there is some freedom in how one defines
a Koopmans functional. Namely, one must choose values for $\eta_{i}$, the
gradient of the energy as a function of the the occupancy of orbital $i$, for
each value of $i$ (modulo the corresponding screening term). In that section,
we briefly introduced the Koopmans integer (KI) approach (eq. 3), but that is
just one of several different ways one can define these gradient terms, and it
is possible to define several variants.

####  2.4.1 KI

In the KI approach, $\eta_{i}$ is chosen as the total energy difference of two
adjacent electronic configurations with integer occupations as given by the
base DFT functional:

$\eta_{i}^{\rm KI}=\left.E^{\rm DFT}\right|_{f_{i}=1}-\left.E^{\rm
DFT}\right|_{f_{i}=0}=\int_{0}^{1}\langle\varphi_{i}|\hat{h}^{\rm
DFT}(f)|\varphi_{i}\rangle df$ (12)

where $\hat{h}^{\mathrm{DFT}}(f)$ is the DFT Hamiltonian with the occupancy of
orbital $i$ constrained to $f$. In this case, the explicit expression for the
unscreened KI Koopmans’ correction to orbital $i$, which we denote as
$\Pi_{i}^{\mathrm{KI}}$, becomes

$\Pi_{i}^{\mathrm{KI}}=E_{\rm Hxc}[\rho-\rho_{i}]-E_{\rm
Hxc}[\rho]+f_{i}\left(E_{\rm Hxc}[\rho-\rho_{i}+n_{i}]-E_{\rm
Hxc}[\rho-\rho_{i}]\right)$ (13)

where $\rho_{i}(\mathbf{r})=f_{i}|\varphi_{i}(\mathbf{r})|^{2}$ and
$n_{i}(\mathbf{r})=|\varphi_{i}(\mathbf{r})|^{2}$. $E_{\rm Hxc}$ denotes the
Hartree and exchange-correlation energy corresponding to the underlying base
functional.

It can be seen that at integer occupations the KI energy correction vanishes;
that is, $\Pi^{\rm KI}_{i}=0$. In other words, for integer occupations the KI
functional preserves the potential energy surface of the base functional! But
while the energy correction is vanishing, the potential is non-vanishing — for
example, the KI potential correction to an occupied variational orbital is

$\displaystyle\frac{\delta\Pi^{\mathrm{KI}}_{i\sigma}}{\delta\rho_{j\sigma^{\prime}}(\mathbf{r})}=\left[-E_{\mathrm{Hxc}}[\rho-
n_{i\sigma}]+E_{\mathrm{Hxc}}[\rho]-\int d\mathbf{r^{\prime}}\
v^{\sigma}_{\mathrm{Hxc}}[\rho](\mathbf{r}^{\prime})\
n_{i\sigma}(\mathbf{r^{\prime}})\right]\delta_{ij}\delta_{\sigma\sigma^{\prime}}$
(14)

(here the spin index $\sigma$ has been decoupled from the orbital index).
Unlike the energy correction in eq. 13, this term is non-zero, which means
that the KI correction will affect the spectral properties of the system while
leaving the total energy unchanged.

####  2.4.2 KIPZ

In the KIPZ approach the slope $\eta_{i}$ is also chosen as the total energy
difference of two adjacent electronic configurations with integer occupations,
but this time using the Perdew-Zunger (PZ) one-electron-self-interaction
corrected (SIC) functional applied to the approximate DFT base functional,

$\eta_{i}^{\rm KIPZ}=\left.E^{\rm PZ}\right|_{f_{i}=1}-\left.E^{\rm
PZ}\right|_{f_{i}=0}=\int_{0}^{1}\langle\varphi_{i}|\hat{h}_{i}^{\rm
PZ}(s)|\varphi_{i}\rangle ds,$ (15)

In this instance, the explicit expression for the unscreened energy correction
corresponding to orbital $i$ (denoted $\Pi_{i}^{\rm KIPZ}$) becomes

$\Pi_{i}^{\rm KIPZ}=-\int_{0}^{f_{i}}\langle\varphi_{i}|\hat{h}^{\rm
DFT}(f)|\varphi_{i}\rangle df+f_{i}\int_{0}^{1}\langle\varphi_{i}|\hat{h}^{\rm
PZ}_{i}(f)|\varphi_{i}\rangle df,$ (16)

where

$\hat{h}_{i}^{\rm PZ}(f)=\hat{h}^{\rm DFT}(f)-\hat{v}^{\rm DFT}_{\rm
Hxc}\left[f|\varphi_{i}(\mathbf{r})|^{2}\right]$ (17)

is the PZ self-interaction correction applied to the $i^{\rm th}$ variational
orbital with constrained occupation $f$, which removes the Hartree-plus-
exchange-correlation potential for that orbital. The KIPZ correction can be
rewritten as

$\Pi^{\rm KIPZ}_{i}=\Pi^{\rm KI}_{i}-f_{i}E_{\rm Hxc}[n_{i}]$ (18)

which makes the physics of this correction clear: it is nothing less than the
KI correction with the addition of a (screened) Perdew-Zunger self-interaction
correction. This added correction removes one-electron self-interaction and
makes the KIPZ functional exact for one-electron systems. In the many-electron
case, it provides different (and typically improved) total energies and forces
than the base functional41, albeit with a screening coefficient for the
Perdew-Zunger correction that is inherited from a spectral condition. More
details are provided in Supporting Information S2.

####  2.4.3 Comparing KI and KIPZ

The KIPZ correction is more computationally expensive than the KI approach,
for the following reasons: we have already mentioned that the KI energy
correction vanishes for integer orbital occupations. Furthermore, for occupied
orbitals, the KI corrective potential is scalar (i.e. it does not have a
spatial dependence) and therefore the total energy is invariant with respect
to unitary rotations of the variational orbitals. Consequently, once the
variational orbitals (and, by extension, the total density) are initialized
they remain unchanged during the entire energy minimization procedure. This
also implies that the screening parameters of occupied variational orbitals
converge instantly (in eq. 10, $\Delta E_{i}$ and $\lambda_{ii}(0,1)$ are
independent of $\alpha^{n}_{i}$ and $\lambda_{ii}(\alpha^{n}_{i},1)$ is linear
in $\alpha^{n}_{i}$). Contrast this with KIPZ: the KIPZ energy does not match
that of the base functional, nor is it invariant with respect to the unitary
rotations of occupied orbitals. This means we must directly minimize the
energy with respect to the shape of the variational orbitals, greatly
increasing the computational cost of these calculations. Furthermore, the KIPZ
ground-state density and variational orbitals are a function of the screening
parameters, which means that the screening parameters must be calculated self-
consistently, further increasing the computational cost.

Despite its additional computational cost, KIPZ has some desirable advantages
over KI: for instance, it is one-electron-self-interaction-free. For this
reason, we also have introduced the “perturbative KIPZ” (pKIPZ) method, where
the KIPZ Hamiltonian is applied non-self-consistently to the KI density and
variational orbitals, as a way of approximating the KIPZ result at reduced
computational cost without significantly compromising the accuracy 46.

It is important to note that the KI functional’s invariance with respect to
unitary rotations of the occupied variational orbitals introduces an ambiguity
in its definition: the variational orbitals are no longer well-defined. This
ambiguity is resolved by formally defining the KI functional as the
$\gamma\rightarrow 0$ limit of the “KI$\gamma$PZ” functional, which is the
KIPZ functional with the PZ contribution to the correction scaled by a
prefactor $\gamma$. This is discussed further in the Supporting Information
S3.1.1.

Finally, we note that the original formulations of Koopmans functionals also
introduced the K and the KPZ functionals 35, 36, 38. These are similar to the
KI and KIPZ functionals, except that the slope $\eta_{i}$ is evaluated at
half-occupation rather than as the total energy difference between integer
occupations. These formulations provide almost identical results but more
cumbersome than their integer counterparts.

####  2.4.4 Total energies and forces with different Koopmans variants

The design of Koopmans functionals focuses on predicting spectral properties.
However, it is worthwhile pausing to consider how accurately these functionals
will predict structural properties (namely, total energies and forces). The KI
functional, as we have already discussed, yields the same total energy — and
by extension, the same forces — as its base functional. The KIPZ functional,
on the other hand, gives total energies and forces that correspond to its base
functional augmented with a screened PZ correction.

There are instances where these two approaches yield significantly different
results. For example, in a study of the geometry of adenine, thymine, and
uracil, the KIPZ@PBE functional predicted bond lengths with a relative mean
absolute error compared to experiment of $0.65\%$, which was slightly better
than PBE0 ($0.76\%$) and PZ@PBE ($0.83\%$), and was markedly better than PBE
($1.63\%$) — and, by extension, KI@PBE 42. That same study showed that the
KIPZ@PBE functional captured the tilt of the amino groups of nucleobases with
respect to their aromatic rings, whereas PBE wrongly predicts a near-planar
structure. However, the addition of a PZ correction does not necessarily
improve structural properties across the board. Ref. 38 compared structural
properties for the reference G2-1 set of molecules, and found that KIPZ@PBE
predicted bond angles less accurately (with a mean relative error of 2.2% for
KIPZ@PBE compared to 1.4% for PBE) despite predicting bond lengths slightly
better (1.5% for KIPZ@PBE compared to 2.3% for PBE).

We stress that these considerations regarding structural properties are
somewhat orthogonal to the Koopmans functional formalism. One should not use
the KI functional to calculate structural properties alone (because the ODD
formalism comes at increased computational cost but provides no change in the
structural properties). If desired, improved geometrical properties and
accurate spectral properties can be simultaneously obtained by combining the
KI correction with a more advanced base functional that predicts structural
properties more reliably.

###  2.5 Important caveats

Before concluding this section, there are a few further important points that
must be made.

####  2.5.1 Restriction to systems with a non-zero band gap

First, the Koopmans formulation is only well-defined for systems with a non-
zero band gap. This is because the Koopmans correction (eq. 8) is defined in
terms of the diagonal elements of the occupation matrix. A band gap (however
small) means that the occupancy matrix is block-diagonal, and can always be
chosen to be the identity for the occupied manifold and zero for the
unoccupied manifold. In the absence of a band gap, the occupancy matrix is not
block-diagonal and a well-defined Koopmans functional would require some
(currently unknown) corrections for the off-diagonal components. While it
would be desirable to derive an off-diagonal correction and to lift this
restriction, the current theory remains powerful — after all, it is in
insulating and semi-conducting systems where DFT exhibits one of its most
striking failures in the underestimation of the band gap.

However, we note that we often we rely on semi-local DFT as the base
functional to define or initialize the variational orbitals. If the base
functional also predicts a non-zero band gap, then the valence and conduction
manifold can be disentangled 61, the occupancy matrix will be block-diagonal,
and the Koopmans correction can immediately be applied. However, if the base
functional wrongly predicts a metallic state, then the valence and conduction
manifolds are not so easily disentangled. In these cases, one might be able to
first employ other base functionals to open a gap (such as DFT +  _U_) or
deploy novel projectability disentanglement methods to separate the valence
and conduction manifolds 62.

The occupancies of variational orbitals $f_{i}$ have been a central quantity
in constructing the Koopmans formalism. This restriction to systems with a
band gap means that these variational orbital occupancies will always be
either 0 or 1, and consequently some terms in the formalism vanish (for
example, the KI correction to the energy; eq. 13) but others do not (for
example, the KI correction to the potential; eq. 14).

####  2.5.2 Empty state localization in the bulk limit

While minimizing the Koopmans energy functional for bulk systems leads to
well-localized occupied orbitals, the same process does not lead to well-
localized empty orbitals. This is because (a) low-lying conduction bands are
often entangled with highly-delocalized nearly-free-electron bands, and (b)
the Koopmans correction to empty states contains a leading Hartree term that
incentivizes delocalization (see Ref. 38). However, the Koopmans correction
ought to be applied to localized orbitals, and vanishes in the limit of
infinitely delocalized states (as discussed in Section 2.2). In light of this,
we typically apply the Koopmans correction non-self-consistently on a
maximally localized Wannier function representation of the empty manifold.
This approach is heuristic but effective, as demonstrated by previous works
43, 45.

####  2.5.3 Symmetries

Because a Koopmans potential $v^{\mathrm{Koopmans}}[\rho,\rho_{i}]$ is
constructed via a variational orbital density, these potentials can break the
translational symmetry of periodic systems. However, the variational orbitals
crucially possess the translational properties of Wannier functions; that is,
for each variational orbital $\varphi_{\mathbf{R}}$ there exists a periodic
replica $\varphi_{\mathbf{R+R^{\prime}}}$ such that

$\varphi_{\mathbf{R}}(\mathbf{r}-\mathbf{R^{\prime}})=\varphi_{\mathbf{R}+\mathbf{R^{\prime}}}(\mathbf{r})$
(19)

where $\mathbf{R}$ and $\mathbf{R}^{\prime}$ can be any pair of Bravais
lattice vectors. Thanks to this property, the collective potential
$\sum_{i}v^{\mathrm{Koopmans}}[\rho,\rho_{i}]|\varphi_{i}\rangle\langle\varphi_{i}|$
inherits the translational symmetry of the overall system and thus it remains
possible to describe the system’s electronic structure with a band-structure
picture. For more details, refer to Ref. 45.

More generally, the orbital-density dependence of Koopmans functionals might
unphysically break the crystal point group symmetry. This is a common feature
of non-rotationally invariant methods that are based on localized orbitals 53,
63. Here, the symmetry of the localized representation plays an important
role, especially in small systems and in the atomic limit. Possible solutions
to this issue have been recently suggested64, and this point is worthy of
further investigation.

##  3 Algorithms and implementation

The formulation of Koopmans functionals, as outlined in the previous section,
is inherently more complex than a “standard” semi-local DFT calculation, and
requires non-standard algorithms and bespoke implementation within electronic-
structure codes. This section describes these algorithms and describes how
Koopmans functionals have been implemented in Quantum ESPRESSO and the open-
source package koopmans.

###  3.1 Orbital optimization

In order to work with Koopmans functionals, we must be able to minimize an
orbital-density-dependent functional. In other words, we must optimize a set
of orbital densities $\\{\rho_{i}\\}$ such that the Koopmans energy functional
(eq. 5) is minimized. This orbital optimization is performed separately for
the occupied and then the empty manifold using an optimization algorithm
similar to that employed in the ensemble DFT approach65: the orbital densities
are parameterized via a set of wavefunctions $\phi_{i}$ and a unitary rotation
matrix $U$, such that $\rho_{i}=|(U\phi)_{i}|^{2}$, and then the energy is
then minimized via the nested loop:

$E=\min_{\\{\phi_{i}\\}}\left(\min_{U}E^{\mathrm{Koopmans}}[\\{|(U\phi)_{i}|^{2}\\}]\right)$
(20)

where in the inner loop the unitary rotation matrix $U$ is optimized (which
leaves the total density unchanged), and in the outer loop the wave functions
are optimized. Both steps are performed using the conjugate-gradient
algorithm. The optimization is performed separately for the occupied and empty
manifolds to ensure that the occupation matrix remains block-diagonal (as
discussed in Section 2.5.1).

One important ingredient in ODD energy minimization is the use of complex
orbitals. Because the ODD energy is not invariant with respect to unitary
rotations of the variational orbitals, it can no longer be assumed (as in the
case for DFT) that the variational orbitals are real, and thus the
aforementioned wave functions $\phi_{i}$ must be complex in order to find the
true minimum of the ODD functional.38, 66, 67, 68, 63

In addition to the generic orbital minimization procedure, we must also
perform constrained minimization calculations (as required by the finite-
difference method for calculating screening parameters; Section 2.3). Here,
the total ODD energy is minimized while removing/adding one electron to a
particular variational orbital. (This gives us $E_{i}(N\pm 1)$ from eq. 9).
This orbital must be frozen during the minimization, otherwise it would morph
into the valence band maximum/conduction band minimum, and one must also
impose the standard orthogonality condition with all other orbitals belonging
to the same spin channel. Image correction methods such as Martina-Tuckerman
or Gygi-Baldereschi 69, 70 must be used to avoid spurious interaction between
charged periodic images. For periodic systems this also means that these
calculations must be performed in a supercell. These charged defect
calculations also require special care in low-dimensional materials 71.
Further details regarding the orbital minimization procedure are presented in
Ref. 40.

###  3.2 The kcp.x code

These orbital minimization algorithms are implemented in the code kcp.x. In
other words, kcp.x can be used to obtain the ground-state energy and the
minimizing set of variational orbitals of an arbitrary system for a given
orbital-density-dependent functional (PZ, KI, or KIPZ).

kcp.x can be used to calculate screening parameters via the finite-difference
approach, and is applicable to both periodic and aperiodic systems. By design,
it does not use $\mathbf{k}$-point sampling for periodic systems, because the
finite-difference approach mandates the use of a supercell (as discussed
above), rendering $\mathbf{k}$-point sampling of the Brillouin zone
superfluous. Instead, the dimensions of the supercell can be used to
effectively sample $\mathbf{k}$-space for bulk systems, and the band structure
for the equivalent primitive cell can be reconstructed at the end of the
calculation using an unfolding procedure 45. Despite the absence of
$\mathbf{k}$-space sampling (which is embarrassingly parallel), the kcp.x code
still uses MPI parallelism: it is parallelized over the plane wave basis. This
allows for the distribution of linear algebra operations and Fourier
transforms across processors.

Because Koopmans functionals are a correction applied on top of a local or
semilocal functional, and these functionals are computationally inexpensive
compared to their ODD counterparts, before commencing orbital minimization
with kcp.x it is efficient to initialize the variational orbitals as Kohn-Sham
orbitals or maximally localized Wannier functions 54. To support the use of
Wannier functions for periodic systems, we have implemented an interface that
takes set of $\mathbf{k}$-indexed Wannier functions from a Wannier90
calculation and maps it to an enlarged set of $\Gamma$-only Wannier functions
defined on the corresponding supercell. Given that kcp.x implements the full
minimization of the ODD functional, in principle one could use the output of
kcp.x to perform geometry optimizations, calculate phonons via the frozen-
phonon method, calculate electron-phonon coupling, model excitons, and so on.

For historical reasons, kcp.x is implemented on top of cp.x, the code within
Quantum ESPRESSO usually responsible for performing Car-Parrinello molecular
dynamics (hence the name “kcp.x”), which already contained algorithms similar
to the direct functional minimization required by Koopmans functionals. It is
important to note that kcp.x is not meant to perform molecular dynamics like
cp.x. The implementation is built on top of version 4.1 of Quantum ESPRESSO.
The modifications made to implement Koopmans functionals are (a) extensive and
(b) of no relevance to the standard functioning of the cp.x code, so these
modifications have not yet been incorporated within the official Quantum
ESPRESSO repository, nor was the private version of the code kept aligned with
subsequent Quantum ESPRESSO releases. Fast-forward to today, and kcp.x has
effectively become a standalone code.

###  3.3 Screening parameters via linear response calculations in reciprocal
space

While the finite-difference approach of kcp.x can provide us with all of the
ingredients to calculate the screening parameters, it is somewhat cumbersome,
since one must perform several constrained DFT and Koopmans calculations, and
for periodic systems these must be performed in a supercell. An alternative to
this approach is to compute the screening coefficients via density-functional
perturbation theory (DFPT) 72.

In this approach, one first approximates the energy as a quadratic function of
the occupation number (which is typically a very good approximation), and the
expression for the screening coefficients reduces to

$\alpha_{i}=\frac{d^{2}E_{\rm DFT}/df_{i}^{2}}{\partial^{2}E_{\rm
DFT}/\partial f_{i}^{2}}=\frac{\langle n_{i}|\epsilon^{-1}f_{\rm
Hxc}|n_{i}\rangle}{\langle n_{i}|f_{\rm Hxc}|n_{i}\rangle},$ (21)

where $\frac{d}{df_{i}}$ ($\frac{\partial}{\partial f_{i}}$) represents
variations that do (do not) account for orbital relaxation,
$\epsilon(\mathbf{r},\mathbf{r}^{\prime})$ is the microscopic dielectric
function of the material, $f_{\rm
Hxc}(\mathbf{r},\mathbf{r}^{\prime})=\delta^{2}E_{Hxc}/\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})$
is the Hartree-plus-exchange-and-correlation kernel, and
$n_{i}(\mathbf{r})=|\varphi_{i}(\mathbf{r})|^{2}$ is the orbital density at
integer occupation44. This can be evaluated by considering the density
response $\Delta^{i}n(\mathbf{r})$ induced in the system by the perturbing
potential $v^{i}_{\rm pert}(\mathbf{r})=\int d\mathbf{r}^{\prime}f_{\rm
Hxc}(\mathbf{r},\mathbf{r}^{\prime})n_{i}(\mathbf{r}^{\prime})$. This
perturbation is the Hartree-plus-exchange-and-correlation potential generated
when adding/removing an infinitesimal fraction of an electron to/from orbital
$i$. One determines $\Delta^{i}n$ self-consistently via DFPT73, and then the
screening parameters are given by

$\alpha_{i}=1+\frac{\langle v^{i}_{\rm pert}|\Delta^{i}n\rangle}{\langle
n_{i}|v^{i}_{\rm pert}\rangle}.$ (22)

Evaluating the screening coefficients within this linear-response approach
only requires quantities available from a $N$-electron calculation, which
means that in the case of periodic solids there is no need for a supercell.
Instead, we can reduce the cost of these calculations by taking advantage of
the translational symmetry of the system 45 and recasting the supercell
problem in a basis of Wannier functions. These Wannier functions take the form
$w^{\mathbf{R}i}(\mathbf{r})$, where the orbital label explicitly denotes the
lattice vector $\mathbf{R}$ of the home cell inside the supercell. In this
basis, the DFPT expression for the screening coefficients (eq. 21) can be
decomposed into a set of independent problems (monochromatic perturbations),
one for each $\mathbf{q}$ point sampling the Brillouin zone of the primitive
cell.73 The now $\mathbf{q}$-dependent charge density variation $\Delta
n^{\mathbf{0}i}_{\mathbf{q}}(\mathbf{r})$ induced by the perturbing potential
$v^{\mathbf{0}i}_{\rm pert,\mathbf{q}}$ is obtained self-consistently via DFPT
(eqs. 15-17 of Ref. 73), and then the screening coefficients are obtained by
summing over $\mathbf{q}$:

$\alpha_{\mathbf{0}i}=1+\frac{\sum_{\mathbf{q}}\langle v^{\mathbf{0}i}_{\rm
pert,\mathbf{q}}|\Delta^{\mathbf{0}i}_{\mathbf{q}}n\rangle}{\sum_{\mathbf{q}}\langle
n^{\mathbf{0}i}_{\mathbf{q}}|v^{\mathbf{0}i}_{\rm pert,\mathbf{q}}\rangle}.$
(23)

The KI Hamiltonian at a particular $\mathbf{k}$ point is then given to second
order by

$H_{ij}^{\mathrm{KI}(2)}(\mathbf{k})=H_{ij}^{\mathrm{DFT}}(\mathbf{k})+\alpha_{\mathbf{0}j}\Delta
H_{ij}^{\mathrm{KI}(2)}(\mathbf{k})$ (24)

where the second-order KI contribution to the Hamiltonian is

$\Delta H_{vv^{\prime}}^{\mathrm{KI}(2)}(\mathbf{k})=-\frac{1}{2}\left\langle
n_{\mathbf{q}}^{\mathbf{0}v}\middle|v_{\rm
pert,\mathbf{q}}^{\mathbf{0}v}\right\rangle\delta_{vv^{\prime}}$ (25)

for valence bands and

$\Delta H_{cc^{\prime}}^{\mathrm{KI}(2)}(\mathbf{k})=-\frac{1}{2}\left\langle
n_{\mathbf{q}}^{\mathbf{0}c}\middle|v_{\rm
pert,\mathbf{q}}^{\mathbf{0}c}\right\rangle\delta_{cc^{\prime}}+\frac{1}{N_{\mathbf{q}}}\sum_{\mathbf{q}}\left\langle
v_{\mathrm{pert},\mathbf{q}}^{\mathbf{0}c^{\prime}}\middle|n_{\mathbf{k},\mathbf{k}+\mathbf{q}}^{cc^{\prime}}\right\rangle$
(26)

for conduction bands, where
$n^{cc^{\prime}}_{\mathbf{k},\mathbf{k}+\mathbf{q}}(\mathbf{r})=(w^{c}_{\mathbf{k}}(\mathbf{r}))^{*}w^{c^{\prime}}_{\mathbf{k}+\mathbf{q}}(\mathbf{r})$;
$w_{\mathbf{k}}^{c}(r)$ the periodic part of the electronic state in the
Wannier gauge. As expected, the KI contribution to the valence bands is
$\mathbf{k}$-independent. The total Hamiltonian is then diagonalized in order
to obtain the canonical eigenstates and energies. Given the fact that the
Hamiltonian is written in a basis of Wannier functions, it is also possible to
employ standard interpolation techniques to obtain the KI eigenvalues at any
arbitrary $\mathbf{k}$-point54.

However, the DFPT approach does come with some limitations. The principal
limitation is that the energy is approximated to second order in the
perturbing potential. In most cases this is very accurate, correctly capturing
the quadratic Hartree contribution and only missing the non-quadratic, higher-
order exchange-correlation contributions.

###  3.4 The kcw.x code

The calculation of screening parameters via DFPT and the subsequent
construction of the Koopmans Hamiltonian and band structure, as described
above, has been implemented in the code kcw.x. Because all of these
calculations are performed in a basis of Wannier functions, this code obtains
Wannier functions via an interface with Wannier90. (The “w” in kcw stands for
“Wannier”.) Because all of these equations are formulated in terms of a
primitive cell with $\mathbf{k}$-point sampling, kcw.x uses MPI to parallelize
over $\mathbf{k}$-points. It also parallelizes over plane-wave orbitals (as
already introduced in the context of kcp.x).

While much of the above applies to periodic systems, kcw.x can still be used
to perform calculations on aperiodic systems. The Wannier function basis still
remains valid, but we no longer have multiple $\mathbf{k}$-points.

kcw.x is part of the official Quantum ESPRESSO distribution (from version 7.1
onward).

###  3.5 Comparing kcp.x and kcw.x

kcp.x and kcw.x implement different Koopmans strategies and, as such, they
have different use-cases, largely defined by their computational scaling. The
two codes scale differently largely due to the fact that kcw.x operates in a
primitive cell while kcp.x operates in a supercell. Calculating one screening
parameter using kcp.x requires multiple SCF calculations, each of which takes
a computational time $T^{\rm SC}$ that roughly scales as
$\mathcal{O}\left(\left({N_{\rm orb}^{\rm SC}}\right)^{3}\right)$, where
$N_{\rm orb}^{\rm SC}$ is the number of orbitals in the supercell. Meanwhile,
calculating one screening parameter using the kcw.x DFPT approach scales as
$T^{\rm PC}\propto N_{\mathbf{q}}N_{\mathbf{k}}{N_{\rm orb}^{\rm PC}}^{3}$.
This is the typical computational time for the SCF cycle
$N_{\mathbf{k}}{N_{\rm orb}^{\rm PC}}^{3}$ times the number of independent
monochromatic perturbations $N_{\mathbf{q}}$. Using the relation $N_{\rm
orb}^{\rm SC}=N_{\mathbf{k}}N_{\rm orb}^{\rm PC}$ and the fact that
$N_{\mathbf{q}}\lesssim N_{\mathbf{k}}$, the ratio between the supercell and
primitive computational times is roughly proportional to $N_{\mathbf{q}}$.
Thus, as the supercell size (or equivalently the number of $\mathbf{q}$-points
in the primitive cell) increases, the kcw.x DFPT approach becomes more and
more computationally efficient.73 For aperiodic systems, $N_{\mathbf{q}}=1$
and the two approaches scale similarly, but with different prefactors.

Note that these scaling relations pertain to the calculation of a single
screening parameter, whereas a full Koopmans workflow requires the calculation
of one screening parameter per unique variational orbital in the system. Here,
the word “unique” is very important; orbitals that are related by symmetry
will share the same screening parameter and therefore the screening does not
need to be recalculated for each orbital. This means that in the worst-case
scenario, where none of the variational orbitals are related by symmetry, the
overall scaling of the workflow has an additional $N_{\rm orb}$ prefactor, but
for many systems (and for periodic systems in particular) the number of unique
variational orbitals in the system can be many times smaller than the total
number of orbitals. Furthermore, (a) the calculation of screening parameters
for separate orbitals is embarrassingly parallelizable, and (b) it is possible
to predict the screening parameters via machine learning, avoiding the need to
repetitively calculate screening parameters altogether 74.

The superior scaling of kcw.x comes at a cost, as it makes two approximations
that kcp.x does not: the DFPT approach expands the total energy only to second
order when computing screening parameters (see Section 3.3), and it does not
optimize the variational orbitals. These are instead defined via Wannier
functions, which often closely resemble the minimizing orbitals of the
Koopmans energy functional. This also means that kcw.x only implements the KI
functional. Without orbital minimization one cannot perform KIPZ calculations,
and pKIPZ would require the PZ kernel (i.e. the second derivative of the PZ
energy with respect to the density), and this is not implemented in common
electronic-structure codes.

###  3.6 Workflow management

Running a Koopmans calculation with either kcp.x or kcw.x requires a few
additional steps compared to a standard semi-local DFT calculation. In this
section, we will focus on the workflows that one needs to perform in order to
complete a Koopmans functional calculation, and how these are publicly
disseminated in open-source form.

Typically, these workflows can be divided into three steps:

  1. 1.

an initialization step, where the variational orbitals are initialized

  2. 2.

the calculation of screening parameters

  3. 3.

a final calculation using the final screening parameters

Depending on the method used for calculating screening parameters (that is,
either finite differences with kcp.x or DFPT with kcw.x), the resulting
workflows look very different. Differences also emerge between calculations on
molecules and solids. For the latter (and for large molecular systems), we
have already seen that maximally localized Wannier functions are typically
used as the variational orbitals (for KI) or as a starting guess for the
variational orbitals (for KIPZ). This necessitates an additional
Wannierization procedure 54 and an interface between Wannier90 and
kcp.x/kcw.x. Meanwhile, for calculating the screening parameters via finite
differences, we must perform a combination of different constrained orbital
minimizations. In all cases, the workflows typically comprise of several if
not dozens of calculations, often involving different electronic structure
codes that must handshake with one another. This can greatly benefit from
automation.

###  3.7 The koopmans package

These workflows are all implemented within the koopmans package. Users
exclusively interact with koopmans, rather than the electronic structure codes
directly (which can include, in addition to kcp.x and kcw.x, pre-existing
codes such as pw.x, pw2wannier90.x, and wannier90.x 75, 76, 77).

Typically, a user provides koopmans with a single input JSON file (some
examples are provided in Supporting Information S4). Based on the settings
provided in this input file, koopmans proceeds through the requested workflow.
Whenever an electronic structure calculation needs to be performed, it
generates the corresponding input file, calls the relevant code, waits for it
to complete, and then parses the output file. Between successive calculations,
it computes intermediate variables, moves and modifies files, _etc._. In other
words, the workflow runner takes care of the banal aspects of performing a
Koopmans calculation, allowing users to concern themselves with scientific
matters (e.g. “what functional do I want to use?”) rather than getting bogged
down in practical details (e.g. “are the Wannier function files in the correct
format for the next calculation to be able to read?”)

The koopmans package is shipped with versions of Quantum ESPRESSO that contain
kcp.x and kcw.x, meaning that it contains everything that is required to
perform Koopmans functional calculations from start to finish.

Further details on the koopmans package can be found in Supporting Information
S5. A step-by-step explanation of the workflows themselves can be found in
Supporting Information S3.

##  4 Example calculations

This Koopmans functional formalism has already proven to be very powerful. In
Ref. 46, Koopmans functionals were found to predict the ionization potentials
of a set of 100 small molecules with comparable/superior accuracy to state-of-
the-art GW approaches. Importantly, Koopmans functionals do not only correct
the ionization potential (i.e. the charged excitation where the most weakly
bound electron is removed) but _any_ single-particle charged excitation. This
was shown for a large set of molecules relevant for photovoltaic applications
41, with Koopmans functionals yielding ultraviolet photoemission spectra that
agree quantitatively with experiment. One can see similar accuracy in the
prediction of band gaps and band structures of periodic systems43, 45, 73; in
a study of prototypical semiconductors and insulators, Koopmans functionals
were found to yield band gaps with a mean absolute error of 0.22 eV, compared
to 0.18 eV when using self-consistent GW with vertex corrections 43.
Importantly, alignment between the valence band edge and the vacuum level was
also very good: across six semiconductors the mean absolute error was 0.19 eV,
compared to 0.39 eV for G0W0 and 0.49 eV for self-consistent GW with vertex
corrections. Finally, Koopmans functionals can accurately describe the
spectral properties of liquids, with the KIPZ functional predicting the
electronic density of states of liquid water with comparable accuracy to self-
consistent GW with vertex corrections 78.

However, all of these calculations were performed by individuals with expert
knowledge of Koopmans functionals and with specific expertise on the codes
that implement them. This final section demonstrates the capabilities of the
koopmans package by way of several examples. All of the following calculations
are possible using a very minimalist input file (see Supporting Information
S4). Note that the following calculations use slightly underconverged
parameters (specifically, the energy cutoff, cell size, and/or the size of the
$k$-point grid). Our focus here is to provide example calculations that can be
reproduced easily by readers, rather than providing high-quality reference
results.

###  4.1 The ionization potential and electron affinity of ozone

First, we present the calculation of the ionization potential and electron
affinity of ozone using koopmans.

This calculation is run with the simple command koopmans ozone.json; the input
and output files for which can be found in Supporting Information S4.1. In
short, this command prompts the full sequence of Quantum ESPRESSO calculations
necessary to initialize the density and variational orbitals, calculate the
screening parameters, and run a final KI calculation. The Quantum ESPRESSO
input and output files for these calculations are all stored in various
subdirectories of the current working directory. In principle one can then
simply parse the quantities of interest from the output files (but there are
easier ways, as explained in Supporting Information S5.2). Refer to Supporting
Information S3.1 for a detailed step-by-step description of this workflow.

| IP | EA |
---|---|---|---
PBE | $7.95$ | $6.17$ | This work
G0W0 | $11.80\pm 0.25$ | $2.34\pm 0.25$ | Ref. 79
scGW0@PBE | $12.57$ |  | Ref. 80
scGW0@HF | $13.16$ |  | Ref. 80
scGW | $12.54$ |  | Ref. 80
qsGW | $13.21$ |  | Ref. 80
CCSD(T) | $12.55$ |  | Ref. 81
KI@[PBE,KS] | $12.52$ | $1.82$ | This work
KI@[PBE,KS] | $12.91$ |  | Ref. 46
experiment | $12.73$ | $2.10$ | Refs. 82, 83, 84, 85
Table 1: The vertical ionization potential (IP) and electron affinity (EA) of
ozone, as calculated using functional, perturbative, and quantum chemistry
methods, as well as experiment. KI@[PBE,KS] denotes the KI Koopmans correction
on top of the PBE base functional, with Kohn-Sham orbitals defining the
variational orbitals, showing excellent performance compared to state-of-the-
art methods. The uncertainties in the G0W0 correspond to the standard
deviation of values reported in Ref. 79, which presents calculations using a
range of codes and basis sets.

The ionization potential (IP) and electron affinity (EA) of ozone, as given by
this calculation, are listed in Table 1, showing the excellent performance of
the KI functional compared to state-of-the-art methods.

###  4.2 The band structure of silicon

koopmans can also perform calculations on bulk systems. Here one typically
performs a Wannierization procedure in order to generate maximally localized
Wannier functions to use as variational orbitals. Running this calculation
gives rise to a similar output to the previous case, with the notable
exception that the initialization procedure now involves Wannierization (see
Supporting Information S4.2).

The band structure that one obtains from this calculation is shown in Figure
3, the band gap is displayed in Table 2, alongside energy differences between
particular symmetry points in the band structure.

The experimental band gap is reproduced with accuracy comparable to self-
consistent GW with vertex corrections, and the energy differences between
symmetry points are reproduced with comparable accuracy to G0W0 (the only
perturbative method for which these data were available). Note that the PBE
and KI valence-to-valence energy differences match. This occurs because the
occupied manifold is comprised of four identical Wannier functions, and thus
the KI correction to these bands amounts to a rigid shift. Contrast this with
the valence-to-conduction energy differences, which are markedly better for
the KI functional.

Figure 3: Band structure of bulk silicon, calculated using the KI functional
with the PBE base functional, MLWFs as variational orbitals, and screening
parameters calculated via finite differences. The PBE band structure is also
plotted for comparison.

| PBE _a_ | G0W0 _b_ | scG$\tilde{\mathrm{W}}$_c_ | KI@[PBE,MLWFs] | KIPZ@PBE _d_ | exp _e_
---|---|---|---|---|---|---
|  |  |  | this work | Ref. 45
$E_{g}$ | $0.49$ | $1.06$ | $1.14$ | $1.16$ | $1.12$ | $1.15$ | 1.17
$\Gamma_{1v}$ $\rightarrow$ | $\Gamma_{25^{\prime}v}$ | $11.97$ | $12.04$ |  | $11.97$ | $11.96$ | $12.09$ |  $12.5$± | $0.6$
$X_{1v}$ $\rightarrow$ | $\Gamma_{25^{\prime}v}$ | $7.82$ |  |  | $7.82$ |  |  | 7.75
$X_{4v}$ $\rightarrow$ | $\Gamma_{25^{\prime}v}$ | $2.85$ | $2.99$ |  | $2.85$ | $2.84$ | $2.86$ | 2.90
$L_{2^{\prime}v}$ $\rightarrow$ | $\Gamma_{25^{\prime}v}$ | $9.63$ | $9.79$ |  | $9.63$ | $9.63$ | $9.74$ |  $9.3$± | $0.4$
$L_{1v}$ $\rightarrow$ | $\Gamma_{25^{\prime}v}$ | $6.98$ | $7.18$ |  | $6.98$ | $6.96$ | $7.04$ |  $6.8$± | $0.2$
$L_{3^{\prime}v}$ $\rightarrow$ | $\Gamma_{25^{\prime}v}$ | $1.19$ | $1.27$ |  | $1.19$ |  |  |  $1.2$± | $0.2$
$\Gamma_{25^{\prime}v}$ $\rightarrow$ | $\Gamma_{15c}$ | $2.48$ | $3.29$ |  | $3.17$ | $3.18$ | $3.20$ |  $3.35$± | $0.01$
$\Gamma_{25^{\prime}v}$ $\rightarrow$ | $\Gamma_{2^{\prime}c}$ | $3.28$ | $4.02$ |  | $3.95$ | $3.94$ | $3.95$ |  $4.15$± | $0.05$
$\Gamma_{25^{\prime}v}$ $\rightarrow$ | $X_{1c}$ | $0.62$ | $1.38$ |  | $1.28$ | $1.30$ | $1.31$ | 1.13
$\Gamma_{25^{\prime}v}$ $\rightarrow$ | $L_{1c}$ | $1.45$ | $2.21$ |  | $2.12$ | $2.12$ | $2.13$ |  $2.04$± | $0.06$
$\Gamma_{25^{\prime}v}$ $\rightarrow$ | $L_{3c}$ | $3.24$ | $4.18$ |  | $3.91$ | $3.93$ | $3.94$ |  $3.9$± | $0.1$
MSE | $0.35$ | $0.02$ |  | $0.01$ | $0.00$ | $0.03$ |  |
MAE | $0.44$ | $0.21$ |  | $0.14$ | $0.16$ | $0.17$ |  |

_a_ this work; _b_ Ref. 86 for $E_{g}$ and Ref. 87 for the transitions; _c_
Ref. 88; _d_ Ref. 45; _e_ Ref. 89

Table 2: The band gap $E_{g}$ and energy differences between symmetry points
in the band structure of bulk silicon (in eV), calculated with various
functional and perturbative approaches. These are compared against
experimental values via the mean signed and mean absolute errors (MSE and MAE
respectively). The calculated values for the band gap and the valence-to-
conduction transitions have been shifted by $-$0.06 eV to account for zero-
point renormalization 90.

###  4.3 The band structure of zinc oxide

The previous example used the finite-difference approach for calculating the
screening parameters. In this final example, we will instead use the DFPT
approach to calculate the band structure of zinc oxide. Refer to Supporting
Information S3.2 for a step-by-step description of what this entails. The
calculated band structure is shown in Figure 4 and the band gaps are listed in
Table 3. The corresponding input and output files are provided in Supporting
Information S4.3. In this instance, the band gap is predicted with better
accuracy than state-of-the-art self-consistent GW with vertex corrections.
This is also true of the average $d$-band energy, although these bands remain
slightly too high in energy relative to experiment. Finally, the bandwidth of
the oxygen $2p$ bands (the six highest-energy occupied bands) is much improved
going from LDA to KI. Note that this is a major departure from the earlier
calculations on silicon, where the KI correction to the occupied bands amounts
to a rigid shift and thus such a bandwidth would not change. In this instance,
these bands comprise of variational orbitals of multiple different characters,
each of which is subject to its own potential shift, and thus the overall band
shape can (and does) change.

Figure 4: Band structure of zinc oxide, calculated using the KI functional
with the LDA base functional, MLWFs as variational orbitals, and screening
parameters calculated via DFPT. The LDA band structure is also plotted for
comparison.

| LDA | G0W0 _a_ | GW0 _b_ | scG$\tilde{\mathrm{W}}$_b_ | KI@[LDA,MLWFs] | exp
---|---|---|---|---|---|---
| this work | Ref. 73 |  |  |  | this work | Ref. 73
$E_{g}$ | $0.53$ | $0.63$ | $1.96$ | $2.84$ | $3.04$ | $3.35$ | $3.52$ | 3.44 _c_
$\langle\varepsilon_{d}\rangle$ | $-5.42$ | $-5.14$ | $-6.1$ | $-6.4$ | $-6.7$ | $-7.06$ | $-6.93$ |  $-$7.5 to $-$8.81 _d_
$\Delta$ | $4.43$ | $4.15$ |  |  |  | $4.93$ | $4.99$ | 5.3 _e_

_a_ Ref. 86; _b_ Ref. 88; _c_ Ref. 91; _d_ Ref. 92; _e_ Ref. 93

Table 3: The band gap $E_{g}$, average $d$-band energy
$\langle\varepsilon_{d}\rangle$, and bandwidth $\Delta$ of bulk zinc oxide
(all in eV), as given by various functional and perturbative approaches, and
compared to experiment. All of the computational values for the band gap have
been shifted by -0.16 eV to account for zero-point renormalization 94, 90. In
contrast to the earlier calculation on bulk silicon, here we used the LDA base
functional to align with Ref. 73. That work used finer parameters than this
work (most notably, a finer $k$-point grid throughout).

##  5 Conclusions

Koopmans functionals are a powerful computational tool for predicting the
spectral properties of atoms, molecules, liquids, and crystalline and
amorphous solids from first principles with a functional approach. This has
already been demonstrated in their ability to calculate the ionization
potentials and electron affinities of small molecules,41, 46 the photoemission
spectra of large molecules41, 42, the electronic structure of liquid water78,
and the band structures and ionization potentials of prototypical
semiconductors and insulators43, 45, 73, all at a level of accuracy comparable
to state-of-the-art many-body perturbation methods.

The newly released koopmans package now makes it possible, for the first time,
for non-experts to use these functionals in their own research. Experts will
also benefit from their calculations becoming much more robust and
reproducible. For more information, we refer the reader to the website
koopmans-functionals.org.

The koopmans package will continue to be maintained and developed. In
particular, Koopmans calculations on periodic systems require the user to
perform a Wannierization of the electronic states, and correctly configuring
this calculation can be onerous. In the near future we will add support for
automated Wannierization 62, 61.

The second focus of ongoing development will be parallelization. Large swathes
of the Koopmans workflow (for example, the calculation of screening
parameters) are embarrassingly parallel. For example, one could calculate a
revised value of the screening parameter for orbital $i$ entirely
independently of the calculation of the screening parameter for orbital $j$.
(This is true for both the finite difference and DFPT schemes.) However,
koopmans performs each calculation in the workflow serially i.e. multiple
calculations are not run simultaneously. (N.B. We are not saying that
individual calculations must be run on a single core; all the codes support
MPI parallelization.) Integration of the workflows within a workflow engine
such as AiiDA would allow us to massively reduce the workflows’ walltimes 95.
Integration within AiiDA would come with the added benefits of AiiDA’s
provenance tracking and error detection/recovery. Combined with the automated
Wannierization and efficient parallelism, high-throughput studies with
Koopmans functionals are just around the corner.

We gratefully acknowledge financial support from the Swiss National Science
Foundation (SNSF – project numbers 179138 and 213082). This research was also
supported by the NCCR MARVEL, a National Centre of Competence in Research,
funded by the Swiss National Science Foundation (grant number 205602). We also
acknowledge support from MaX – MAterials design at the eXascale – a European
Centre of Excellence funded by the European Union’s program HORIZON-EUROHPC-
JU-2021-COE-01 (grant number 101093374). NLN acknowledges support by the
Vietnam National Foundation for Science and Technology Development (NAFOSTED)
under grant number 103.02-2021.95. ID acknowledges support from the Division
of Materials Research of the National Science Foundation through CAREER Award
number DMR-1654625. We thank Susi Lehtola for his helpful comments on the
manuscript.

The Supporting Information contains a derivation of the functional form of
Koopmans functionals, further details on the KIPZ functional, a detailed
description of the Koopmans workflows, example input and output koopmans
files, and additional details regarding the koopmans package. All of the input
and output files related to this paper can also be found on Materials Cloud at
10.24435/materialscloud:9w-sp.

Supporting Information

##  S1 Derivation of the functional form of Koopmans functionals

The brief derivation of the functional form of Koopmans functionals is as
follows: let us assume a functional of the form

$E^{\text{Koopmans}}=E^{\mathrm{DFT}}+\sum_{i}\Pi_{i}$ (S1)

If we take the derivative with respect to the occupancy of the $j^{\rm th}$
variational orbital then we have

$\eta_{j}=\left.\frac{dE^{DFT}}{df_{j}}\right\rvert_{f_{j}=f}+\left.\frac{d\Pi_{j}}{df_{j}}\right\rvert_{f_{j}=f}=\langle\varphi_{j}\rvert\hat{h}^{\text{DFT}}(f)\lvert\varphi_{j}\rangle+\left.\frac{d\Pi_{j}}{df_{j}}\right\rvert_{f_{j}=f}$
(S2)

where we assumed that the cross-term derivatives $d\Pi_{i}/df_{j}$ vanish, and
because $E^{\text{Koopmans}}$ ought to be linear in $f_{j}$, we replaced its
derivative with some yet-to-be deterimined constant $\eta_{j}$. For the second
equality we invoked Janak’s theorem, and $f$ is some number between 0 and 1.

Assuming that the energy correction $\Pi_{j}$ is zero at integer occupancies,
is independent of $f_{i}$ for $i\neq j$, and neglecting for the moment any
orbital relaxation as the orbital occupancies change, it follows that

$\Pi^{u}_{j}=-\int_{0}^{f_{j}}\langle\varphi_{j}\rvert\hat{h}^{\text{DFT}}(f)\lvert\varphi_{j}\rangle
df+f_{j}\eta_{j}=-\left(E^{\mathrm{DFT}}[\rho]-E^{\mathrm{DFT}}[\rho-\rho_{i}]\right)+f_{j}\eta_{j}$
(S3)

where the $u$ superscript denotes the fact that we neglected orbital
relaxation, and thus this term is “unscreened”. To account for this screening
we must introduce some screening parameters $\\{\alpha_{i}\\}$ such that
$\Pi_{j}=\alpha_{j}\Pi^{u}_{j}$. Having done this, we arrive at eq. 5, the
final result.

##  S2 KIPZ details

In previous works, KIPZ has been presented in slightly different ways. In eq.
27 of Ref. 38, KIPZ was introduced as

$\displaystyle\Pi_{i}^{\mathrm{KIPZ}}=$
$\displaystyle-\int_{0}^{f_{i}}\langle\varphi_{i}|\hat{H}^{\mathrm{PZ}}_{i}(s)|\varphi_{i}\rangle
ds+f_{i}\int_{0}^{1}\langle\varphi_{i}|\hat{H}^{\mathrm{PZ}}_{i}(s)|\varphi_{i}\rangle
ds-E_{\mathrm{Hxc}}[\rho_{i}]$ (S4)

where
$\hat{H}_{i}^{\mathrm{PZ}}(s)=\hat{H}^{\mathrm{DFT}}(s)-\hat{v}_{\mathrm{Hxc}}^{\mathrm{DFT}}[s|\varphi_{i}(\mathbf{r})|^{2}]$.
(Ref. 38 included an erroneous sum over $i$ in the definition of this
Hamiltonian.) Meanwhile, eq. 6 of Ref 43 defined KIPZ as

$\displaystyle\Pi_{i}^{\mathrm{KIPZ}}=$
$\displaystyle-\int_{0}^{f_{i}}\langle\varphi_{i}|\hat{H}^{\mathrm{DFT}}(s)|\varphi_{i}\rangle
ds+f_{i}\int_{0}^{1}\langle\varphi_{i}|\hat{H}_{i}^{\mathrm{PZ}}(s)|\varphi_{i}\rangle
ds$ (S5)

and in that same paper it was also stated that

$\displaystyle\Pi_{i}^{\mathrm{KIPZ}}=$
$\displaystyle\Pi_{i}^{\mathrm{KI}}-f_{i}E_{\mathrm{Hxc}}[n_{i}]$ (S6)

One can prove that these three definitions are equivalent via the identity

$\displaystyle\int_{0}^{f}\langle\varphi_{i}|v_{\mathrm{Hxc}}[sn_{i}]|\varphi_{i}\rangle
ds=E_{\mathrm{Hxc}}[fn_{i}],$ (S7)

from which it follows that

$\displaystyle\int_{0}^{f_{i}}\langle\varphi_{i}|\hat{H}^{\mathrm{PZ}}_{i}(s)|\varphi_{i}\rangle
ds$ $\displaystyle=$
$\displaystyle\int_{0}^{f_{i}}\langle\varphi_{i}|\hat{H}^{\mathrm{DFT}}(s)|\varphi_{i}\rangle
ds-\int_{0}^{f_{i}}\langle\varphi_{i}|v_{\mathrm{Hxc}}[sn_{i}]|\varphi_{i}\rangle
ds$ $\displaystyle=$
$\displaystyle\int_{0}^{f_{i}}\langle\varphi_{i}|\hat{H}^{\mathrm{DFT}}(s)|\varphi_{i}\rangle
ds-E_{\mathrm{Hxc}}[f_{i}n_{i}]$ (S8)

which proves that eqs. S4 and S5 are equivalent. Furthermore,

$\displaystyle-\int_{0}^{f_{i}}\langle\varphi_{i}|\hat{H}^{\mathrm{DFT}}(s)|\varphi_{i}\rangle
ds+f_{i}\int_{0}^{1}\langle\varphi_{i}|\hat{H}_{i}^{\mathrm{PZ}}(s)|\varphi_{i}\rangle
ds$ $\displaystyle=$
$\displaystyle-\int_{0}^{f_{i}}\langle\varphi_{i}|\hat{H}^{\mathrm{DFT}}(s)|\varphi_{i}\rangle
ds+f_{i}\int_{0}^{1}\langle\varphi_{i}|\hat{H}_{i}^{\mathrm{DFT}}(s)|\varphi_{i}\rangle
ds-
f_{i}\int_{0}^{1}\langle\varphi_{i}|v_{\mathrm{Hxc}}[sn_{i}]|\varphi_{i}\rangle
ds$ $\displaystyle=$
$\displaystyle\Pi_{i}^{\mathrm{KI}}-f_{i}E_{\mathrm{Hxc}}[n_{i}]$ (S9)

and thus eqs. S5 and S6 are equivalent.

In the unscreened case, the KIPZ functional as defined above is equivalent to
the KI correction applied to an unscreened PZ base functional (i.e. KI@PZ).
However, in the general case when screening is accounted for, KIPZ and KI@PZ
are not equivalent. Instead, the PZ corrections incorporated within the KIPZ
functional each inherit their own screening coefficient from the generalized
PWL condition. This is desirable, because scaling down the PZ correction has
been shown to improve energetics and thermochemistry 96, 97, 98, but it would
be interesting to explore alternative prescriptions for the scaling of the PZ
correction that are decoupled from the generalized PWL condition.

##  S3 Details of the Koopmans workflows

This appendix contains a detailed breakdown of the two key Koopmans workflows:
one for calculations where the screening parameters are calculated via finite
differences (Figure S1), and the other via DFPT (Figure S2)

###  S3.1 The finite-difference workflow

In this workflow, we calculate screening parameters via the method described
in Section 2.3 of the main text.

Initialisation (molecules) Initialisation (solids) Calculating screening
parameters Final calculation Postprocessing (solids) kcp pw w90 bespoke code
quantity of interest DFT PZ unitary rotation or PZ DFT (primitive cell)
wannierize fold to supercell $\alpha^{n}$KI/$\alpha^{n}$KIPZ
DFT/$\alpha^{n}$KIPZ ($N-1$) DFT/$\alpha^{n}$KIPZ ($N+1$) $E_{i}(N-1)$
$\lambda_{ii}(0,1)$ $\lambda_{ii}({\alpha^{n}},1)$ $E(N)$
$\lambda_{ii}(\alpha^{n},0)$ $\lambda_{ii}(0,0)$ $E_{i}(N+1)$
$\alpha^{n+1}_{i\in\text{filled}}$ $\alpha^{n+1}_{i\in\text{empty}}$
$\\{\alpha_{i}\\}$ not converged $\alpha$KI/$\alpha$KIPZ $\\{\alpha_{i}\\}$
converged unfold to primitive one per filled orbital (index $i$) one per empty
orbital (index $i$)abcdefghijklmnopqrst

Figure S1: The finite-difference workflow. Individual nodes represent a
calculation performed with Quantum ESPRESSO, except for the red nodes, which
represent key quantities of interest that we extract from a preceding
calculation. Individual nodes are explained more fully in the text.

####  S3.1.1 Initialization

The first step in this workflow is the initialization of the density and the
variational orbitals. Depending on the system and functional in question, this
can look quite different. For molecules, one can start with a DFT calculation
to obtain the ground-state Kohn-Sham eigenvalues (node a in Figure S1). This
is typical of KI calculations, which share the same ground-state density as
the base DFT functional, in which case the final density has already been
determined by this very first calculation. Note, however, that a unitary
rotation of the occupied Kohn-Sham orbital densities leaves the total density
(and therefore the total energy) unchanged, which means that the variational
orbitals from a DFT calculation are not uniquely defined. In order to resolve
this issue, one then performs a unitary rotation of the occupied Kohn-Sham
variational orbitals to minimize the PZ energy (node b). This gives us a
unique set of variational orbitals, while leaving the total density unchanged,
and fulfils the definition of the KI functional as the $\gamma\rightarrow 0$
limit of the “KI$\gamma$PZ” functional (as introduced in Section 2.4.3) in the
main text. If one is using the KIPZ functional, it is better to initialize the
density and the variational orbitals by performing a full PZ calculation (node
c). In contrast to the previous KI approach where the DFT and KI ground-state
densities match, the KIPZ and PZ ground-state densities and variational
orbitals are similar but not identical, so the PZ solution serves as a
suitable initial guess for a KIPZ calculation. Note that all of the above
calculations are performed with the $\Gamma$-point-only kcp.x code.

For solids, the approach for initializing the density and variational orbitals
is very different. Here, we take advantage of the periodicity of the lattice
by initializing the variational orbitals using maximally localized Wannier
functions (MLWFs). This approach is justified by the Wannier-like character of
the true minimizing orbitals 45. Practically, the Wannierization procedure
involves a DFT calculation with the pw.x code (node d), followed by a
Wannierization procedure using the Wannier90 and pw2wannier90.x codes (node
e).

There are two important points when it comes to the Wannierization. The first
is that the occupied and the empty manifolds must be Wannierized separately.
This guarantees that the occupancy matrix is diagonal in the basis of
variational orbitals, as required by Koopmans functionals (see Section 2.5 of
the main text). The second important point is that mixing bands that are far
apart in energy-space is generally detrimental to the Koopmans results. To
avoid this, each block of bands that are well-separated in energy-space are
Wannierized separately, preventing inter-block mixing during the
Wannierization procedure. This is a similar but cruder approach to the so-
called dually-localized Wannier functions, where the Wannier functions
minimize a localization criteria that is a mix of spatial and energy
localization 99.

The one final task is to map the Wannier functions in the $k$-sampled
primitive cell to the equivalent $\Gamma$-point-only supercell in a format
readable by the kcp.x code that will handle the subsequent calculation of the
screening parameters (node f). In this procedure the supercell dimensions
match those of the $k$-grid used during the initialization.

####  S3.1.2 Calculating the screening parameters

Having initialized the density and the variational orbitals, the next task to
perform is the calculation of the screening parameters. To this end, let us
restate equations 10 and 11 from the main text, with which we calculate these
parameters:

$\alpha^{n+1}_{i}=\alpha^{n}_{i}\frac{\Delta
E^{\text{Koopmans}}_{i}-\lambda_{ii}(0,1)}{\lambda_{ii}(\alpha^{n}_{i},1)-\lambda_{ii}(0,1)};\qquad\Delta
E^{\text{Koopmans}}_{i}=E^{\text{Koopmans}}(N)-E^{\text{Koopmans}}_{i}(N-1)$
(S10)

for occupied orbitals and

$\alpha^{n+1}_{i}=\alpha^{n}_{i}\frac{\Delta
E^{\text{Koopmans}}_{i}-\lambda_{ii}(0,0)}{\lambda_{ii}(\alpha^{n}_{i},0)-\lambda_{ii}(0,0)};\qquad\Delta
E^{\text{Koopmans}}_{i}=E^{\text{Koopmans}}_{i}(N+1)-E^{\text{Koopmans}}(N)$
(S11)

for empty orbitals. In order to calculate these screening parameters, we
therefore require three calculations. The first calculation is a KI or KIPZ
calculation with using a trial screening parameter $\alpha^{0}$ (node g). This
gives us access to the energy of the $N$-electron system $E(N)$ (node m) as
well as all of the requisite expectation values Koopmans Hamiltonian on the
variational orbitals $\lambda_{ii}(\alpha,f)$ (nodes k, l, n, and o). The
second and third calculations (nodes h and i) are calculations on the $N\pm
1$-electron systems, where orbital $i$ is frozen and its occupancy is fixed to
0 (in the case of occupied orbitals) or 1 (empty orbitals). These calculations
yield the total energies $E_{i}(N\pm 1)$ (nodes j and p). For these
calculations in particular, ensuring that there is no spurious interactions
between images is crucial (because now we have a charged defect in our
system). This requires the use of both a sufficiently large supercell and a
correction scheme such as Gygi-Baldereschi 70. Note that since KI yields the
same total energies as the base functional, these two calculations can be
performed at the DFT level when performing the KI workflow.

Having performed these three calculations (nodes g-i) and extracted all of the
requisite information (nodes j-p), we can then calculate the screening
parameters (nodes q and r) according to the above equations. If the screening
parameters are converged, we can then proceed to the final calculation; if
not, the process is repeated.

We note that this iterative procedure almost universally converges very
quickly. Indeed, for the KI functional and with occupied orbitals, it is
guaranteed to converge instantly. This is because, as mentioned earlier,
$\Delta E_{i}$ is independent of $\alpha$, as are occupied variational
orbitals, and consequently $\lambda^{\mathrm{KI}}_{ii}(\alpha,1)$ is linear in
$\alpha$. This is not the case for empty orbitals for the KI functional (for
which $\lambda^{\mathrm{KI}}_{ii}(\alpha,0)$ is not strictly linear) or for
the KIPZ functional (where additionally $\Delta E_{i}$ is dependent on
$\alpha$). Even for these functionals, the screening parameters tend to
converge in a few iterations.

####  S3.1.3 The final calculation and postprocessing for solids

We now perform a KI or KIPZ calculation with the finalized set of screening
parameters (node s). For a molecular system we are now done: the KI/KIPZ
calculation yields a Hamiltonian in the basis of variational orbitals which we
diagonalize to extract the quasiparticle energies.

However, for a calculation on a periodic system one final step is required.
This is because all of the preceding kcp.x calculations have been performed in
a $\Gamma$-point-only supercell. In order to extract the band structure, we
must now unfold the band structure by taking advantage of the MLWF basis, as
described in Ref. 45. This step is performed within python by the koopmans
workflow manager itself. One trick that we can perform at this stage is
“smooth interpolation”. In the supercell, our Hamiltonian in the basis of
Wannier functions is given by

$h^{\mathrm{DFT}}_{mn}(\mathbf{R})+v^{\mathrm{Koopmans}}_{mn}(\mathbf{R})$
(S12)

The Koopmans potential is very smooth and slowly-varying in
$\mathbf{k}$-space, applying an almost-constant shift to the Kohn-Sham DFT
bands. Consequently, the dominant contribution to the dispersion of the bands
comes from the DFT Hamiltonian, and it makes sense to construct the
$\mathbf{k}$-indexed Hamitonian as

$h_{mn}(\mathbf{k})=\sum_{\mathbf{R}^{\prime}}e^{i\mathbf{k}\cdot\mathbf{R}^{\prime}}h_{mn}^{\mathrm{DFT}}(\mathbf{R^{\prime}})+\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}v_{mn}^{\mathrm{Koopmans}}(\mathbf{R})$
(S13)

where now $\\{\mathbf{R}^{\prime}\\}$ corresponds to a much larger supercell
or, equivalently, a much denser $\mathbf{k}$-point grid . The advantage of
this strategy is that it improves the interpolation of the band structure at
very little computational cost. Suppose we perform a smooth interpolation with
a grid twice as fine as the default grid. The only additional computational
cost in this instance is having to generate the DFT Hamiltonian in the Wannier
basis for this finer grid (i.e. we repeat nodes d and e). This only represents
a small fraction of the total workflow, and includes only DFT and not ODDFT
calculations, so it only fractionally increases the total computational cost.
Contrast this with the alternative, where one could perform the entire
calculation with a grid twice as fine. This would require us to perform (among
other things) calculations on a supercell containing eight times as many
atoms, drastically increasing the computational cost of the workflow as a
whole. For more details on the smooth interpolation procedure, refer to Ref.
45.

###  S3.2 The DFPT workflow

The DFPT workflow is depicted in Figure S2. It is simpler than the finite-
difference procedure, because orbital relaxation is not implemented, and
instead Wannier functions are used as approximations to the true variational
orbitals. This means that only the KI and pKIPZ functionals can be used in
this scheme, means that the screening parameters do not need to be calculated
self-consistently, and makes the Koopmans functional effectively a post-
processing step on top of a DFT calculation.

Initialisation (molecules) Initialisation (solids) Calculating screening
parameters Final calculation Postprocessing (solids) Initialisation
Calculating screening parameters Final calculation pw0 w90 kcw quantity of
interest DFT wannierize Wannier function conversion DFPT calculation
$\\{\alpha_{i}\\}$ $\alpha_{0}$KI/$\alpha_{0}$KIPZ uvwxyz

Figure S2: The DFPT workflow. Individual nodes are explained in the text.

####  S3.2.1 Initialization

The initialization procedure for the DFPT workflow is very similar to that for
solids in the finite-difference workflow. A primitive cell calculation pw.x
calculation (node u) is followed by a Wannierization procedure in order to
define the density and the variational orbitals (node v). The one difference
is that now, instead of using kcw.x to map these Wannier functions to a
supercell that is readable by kcp.x, we use kcw.x to convert the Wannier
functions to more convenient format for subsequent calculations (node w). Note
that kcw.x does not map to a supercell because this workflow operates entirely
within the primitive cell with $k$-point sampling.

####  S3.2.2 Calculating the screening parameters

DFPT calculations evaluating equation 21 (from the main text) are then
performed by a single kcw.x run (node x). These calculations yield the
screening parameters (node y).

####  S3.2.3 The final calculation and postprocessing for solids

Having calculated the screening parameters, the Koopmans Hamiltonian is
constructed in the basis of Wannier functions and then diagonalized, resulting
in the full band structure of the system at hand (node z). Like the previous
calculations, this is performed using kcw.x.

##  S4 Example input and output files

Below are example koopmans input and output files for ozone, silicon, and zinc
oxide. All of these files can also be found on Materials Cloud at
10.24435/materialscloud:9w-sp.

###  S4.1 Ozone

An example minimal input file for ozone is as follows

⬇

{

”workflow”: {

”functional”: ”ki”,

”method”: ”dscf”,

”init_orbitals”: ”kohn-sham”,

”n_max_sc_steps”: 5,

”keep_tmpdirs”: false,

”pseudo_library”: ”sg15”

},

”atoms”: {

”cell_parameters”: {

”periodic”: false,

”vectors”: [[14.1738, 0.0, 0.0],

[0.0, 12.0, 0.0],

[0.0, 0.0, 12.66]],

”units”: ”angstrom”

},

”atomic_positions”: {

”positions”: [[”O”, 7.0869, 6.0, 5.89],

[”O”, 8.1738, 6.0, 6.55],

[”O”, 6.0, 6.0, 6.55]],

”units”: ”angstrom”

}

},

”calculator_parameters”: {

”ecutwfc”: 65.0,

”nbnd”: 10

}

}

This input file contains several blocks. The workflow block allows the user to
specify the details of the workflow. Here we can see we are performing a KI
calculation (line 3) calculating the screening parameters via the finite-
difference procedure (line 4), and using the Kohn-Sham orbitals to initialize
our variational orbitals (line 5; this is common practice for molecules). The
atoms block (lines 10-23) contains standard keywords specifying the system
configuration, such as the cell_parameters and atomic_positions. These mirror
the equivalent blocks in Quantum ESPRESSO input files (albeit in JSON format).
Finally, the calculator_parameters block allows the user to specify settings
specific to a particular code (e.g. a w90 subblock for specifying Wannier90
settings). In this instance we are providing a particular energy cutoff (line
26) and specifying the total number of orbitals to compute (line 27).

The output of koopmans ozone.json, which prompts a sequence of Quantum
ESPRESSO calculations necessary to initialize the density and variational
orbitals (lines 16-22), calculate the screening parameters (lines 24-68), and
run a final KI calculation (lines 85-87).

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in ozone.bib in work involving this calculation

Initialization of density and variational orbitals

==================================================

Running init/dft_init_nspin1… done

Running init/dft_init_nspin2_dummy… done

Running init/dft_init_nspin2… done

Overwriting the variational orbitals with Kohn-Sham orbitals

Copying the spin-up variational orbitals over to the spin-down channel

Calculating screening parameters

================================

SC iteration 1

————–

Running calc_alpha/iteration_1/ki… done

Orbital 1

———

Running calc_alpha/iteration_1/orbital_1/dft_n-1… done

Orbital 2

———

Running calc_alpha/iteration_1/orbital_2/dft_n-1… done

Orbital 3

———

Running calc_alpha/iteration_1/orbital_3/dft_n-1… done

Orbital 4

———

Running calc_alpha/iteration_1/orbital_4/dft_n-1… done

Orbital 5

———

Running calc_alpha/iteration_1/orbital_5/dft_n-1… done

Orbital 6

———

Running calc_alpha/iteration_1/orbital_6/dft_n-1… done

Orbital 7

———

Running calc_alpha/iteration_1/orbital_7/dft_n-1… done

Orbital 8

———

Running calc_alpha/iteration_1/orbital_8/dft_n-1… done

Orbital 9

———

Running calc_alpha/iteration_1/orbital_9/dft_n-1… done

Orbital 10

———-

Running calc_alpha/iteration_1/orbital_10/pz_print… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[2 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 3 … 8 9 10

0 -0.477344 -0.920497 -1.124088 … -0.888608 -1.050159 0.711209

[1 rows x 10 columns]

SC iteration 2

————–

Running calc_alpha/iteration_2/ki_nspin1_dummy… done

Running calc_alpha/iteration_2/ki_nspin1… done

Running calc_alpha/iteration_2/ki_nspin2_dummy… done

Running calc_alpha/iteration_2/ki_nspin2… done

Orbital 10

———-

Running calc_alpha/iteration_2/orbital_10/pz_print… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

2 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[3 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 … 9 10

0 -4.773441e-01 -9.204971e-01 … -1.050159e+00 7.112091e-01

1 2.131628e-14 -2.415845e-13 … -4.440892e-14 1.776357e-14

[2 rows x 10 columns]

Screening parameters have been converged

Final KI calculation

====================

Running final/ki_final… done

Workflow complete

...

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in ozone.bib in work involving this calculation

Initialization of density and variational orbitals

==================================================

Running init/dft_init_nspin1… done

Running init/dft_init_nspin2_dummy… done

Running init/dft_init_nspin2… done

Overwriting the variational orbitals with Kohn-Sham orbitals

Copying the spin-up variational orbitals over to the spin-down channel

Calculating screening parameters

================================

SC iteration 1

————–

Running calc_alpha/iteration_1/ki… done

Orbital 1

———

Running calc_alpha/iteration_1/orbital_1/dft_n-1… done

Orbital 2

———

Running calc_alpha/iteration_1/orbital_2/dft_n-1… done

Orbital 3

———

Running calc_alpha/iteration_1/orbital_3/dft_n-1… done

Orbital 4

———

Running calc_alpha/iteration_1/orbital_4/dft_n-1… done

Orbital 5

———

Running calc_alpha/iteration_1/orbital_5/dft_n-1… done

Orbital 6

———

Running calc_alpha/iteration_1/orbital_6/dft_n-1… done

Orbital 7

———

Running calc_alpha/iteration_1/orbital_7/dft_n-1… done

Orbital 8

———

Running calc_alpha/iteration_1/orbital_8/dft_n-1… done

Orbital 9

———

Running calc_alpha/iteration_1/orbital_9/dft_n-1… done

Orbital 10

———-

Running calc_alpha/iteration_1/orbital_10/pz_print… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[2 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 3 … 8 9 10

0 -0.477344 -0.920497 -1.124088 … -0.888608 -1.050159 0.711209

[1 rows x 10 columns]

SC iteration 2

————–

Running calc_alpha/iteration_2/ki_nspin1_dummy… done

Running calc_alpha/iteration_2/ki_nspin1… done

Running calc_alpha/iteration_2/ki_nspin2_dummy… done

Running calc_alpha/iteration_2/ki_nspin2… done

Orbital 10

———-

Running calc_alpha/iteration_2/orbital_10/pz_print… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

2 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[3 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 … 9 10

0 -4.773441e-01 -9.204971e-01 … -1.050159e+00 7.112091e-01

1 2.131628e-14 -2.415845e-13 … -4.440892e-14 1.776357e-14

[2 rows x 10 columns]

Screening parameters have been converged

Final KI calculation

====================

Running final/ki_final… done

Workflow complete

...

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in ozone.bib in work involving this calculation

Initialization of density and variational orbitals

==================================================

Running init/dft_init_nspin1… done

Running init/dft_init_nspin2_dummy… done

Running init/dft_init_nspin2… done

Overwriting the variational orbitals with Kohn-Sham orbitals

Copying the spin-up variational orbitals over to the spin-down channel

Calculating screening parameters

================================

SC iteration 1

————–

Running calc_alpha/iteration_1/ki… done

Orbital 1

———

Running calc_alpha/iteration_1/orbital_1/dft_n-1… done

Orbital 2

———

Running calc_alpha/iteration_1/orbital_2/dft_n-1… done

Orbital 3

———

Running calc_alpha/iteration_1/orbital_3/dft_n-1… done

Orbital 4

———

Running calc_alpha/iteration_1/orbital_4/dft_n-1… done

Orbital 5

———

Running calc_alpha/iteration_1/orbital_5/dft_n-1… done

Orbital 6

———

Running calc_alpha/iteration_1/orbital_6/dft_n-1… done

Orbital 7

———

Running calc_alpha/iteration_1/orbital_7/dft_n-1… done

Orbital 8

———

Running calc_alpha/iteration_1/orbital_8/dft_n-1… done

Orbital 9

———

Running calc_alpha/iteration_1/orbital_9/dft_n-1… done

Orbital 10

———-

Running calc_alpha/iteration_1/orbital_10/pz_print… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[2 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 3 … 8 9 10

0 -0.477344 -0.920497 -1.124088 … -0.888608 -1.050159 0.711209

[1 rows x 10 columns]

SC iteration 2

————–

Running calc_alpha/iteration_2/ki_nspin1_dummy… done

Running calc_alpha/iteration_2/ki_nspin1… done

Running calc_alpha/iteration_2/ki_nspin2_dummy… done

Running calc_alpha/iteration_2/ki_nspin2… done

Orbital 10

———-

Running calc_alpha/iteration_2/orbital_10/pz_print… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

2 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[3 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 … 9 10

0 -4.773441e-01 -9.204971e-01 … -1.050159e+00 7.112091e-01

1 2.131628e-14 -2.415845e-13 … -4.440892e-14 1.776357e-14

[2 rows x 10 columns]

Screening parameters have been converged

Final KI calculation

====================

Running final/ki_final… done

Workflow complete

...

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in ozone.bib in work involving this calculation

Initialization of density and variational orbitals

==================================================

Running init/dft_init_nspin1… done

Running init/dft_init_nspin2_dummy… done

Running init/dft_init_nspin2… done

Overwriting the variational orbitals with Kohn-Sham orbitals

Copying the spin-up variational orbitals over to the spin-down channel

Calculating screening parameters

================================

SC iteration 1

————–

Running calc_alpha/iteration_1/ki… done

Orbital 1

———

Running calc_alpha/iteration_1/orbital_1/dft_n-1… done

Orbital 2

———

Running calc_alpha/iteration_1/orbital_2/dft_n-1… done

Orbital 3

———

Running calc_alpha/iteration_1/orbital_3/dft_n-1… done

Orbital 4

———

Running calc_alpha/iteration_1/orbital_4/dft_n-1… done

Orbital 5

———

Running calc_alpha/iteration_1/orbital_5/dft_n-1… done

Orbital 6

———

Running calc_alpha/iteration_1/orbital_6/dft_n-1… done

Orbital 7

———

Running calc_alpha/iteration_1/orbital_7/dft_n-1… done

Orbital 8

———

Running calc_alpha/iteration_1/orbital_8/dft_n-1… done

Orbital 9

———

Running calc_alpha/iteration_1/orbital_9/dft_n-1… done

Orbital 10

———-

Running calc_alpha/iteration_1/orbital_10/pz_print… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[2 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 3 … 8 9 10

0 -0.477344 -0.920497 -1.124088 … -0.888608 -1.050159 0.711209

[1 rows x 10 columns]

SC iteration 2

————–

Running calc_alpha/iteration_2/ki_nspin1_dummy… done

Running calc_alpha/iteration_2/ki_nspin1… done

Running calc_alpha/iteration_2/ki_nspin2_dummy… done

Running calc_alpha/iteration_2/ki_nspin2… done

Orbital 10

———-

Running calc_alpha/iteration_2/orbital_10/pz_print… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

2 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[3 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 … 9 10

0 -4.773441e-01 -9.204971e-01 … -1.050159e+00 7.112091e-01

1 2.131628e-14 -2.415845e-13 … -4.440892e-14 1.776357e-14

[2 rows x 10 columns]

Screening parameters have been converged

Final KI calculation

====================

Running final/ki_final… done

Workflow complete

...

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in ozone.bib in work involving this calculation

Initialization of density and variational orbitals

==================================================

Running init/dft_init_nspin1… done

Running init/dft_init_nspin2_dummy… done

Running init/dft_init_nspin2… done

Overwriting the variational orbitals with Kohn-Sham orbitals

Copying the spin-up variational orbitals over to the spin-down channel

Calculating screening parameters

================================

SC iteration 1

————–

Running calc_alpha/iteration_1/ki… done

Orbital 1

———

Running calc_alpha/iteration_1/orbital_1/dft_n-1… done

Orbital 2

———

Running calc_alpha/iteration_1/orbital_2/dft_n-1… done

Orbital 3

———

Running calc_alpha/iteration_1/orbital_3/dft_n-1… done

Orbital 4

———

Running calc_alpha/iteration_1/orbital_4/dft_n-1… done

Orbital 5

———

Running calc_alpha/iteration_1/orbital_5/dft_n-1… done

Orbital 6

———

Running calc_alpha/iteration_1/orbital_6/dft_n-1… done

Orbital 7

———

Running calc_alpha/iteration_1/orbital_7/dft_n-1… done

Orbital 8

———

Running calc_alpha/iteration_1/orbital_8/dft_n-1… done

Orbital 9

———

Running calc_alpha/iteration_1/orbital_9/dft_n-1… done

Orbital 10

———-

Running calc_alpha/iteration_1/orbital_10/pz_print… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_1/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[2 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 3 … 8 9 10

0 -0.477344 -0.920497 -1.124088 … -0.888608 -1.050159 0.711209

[1 rows x 10 columns]

SC iteration 2

————–

Running calc_alpha/iteration_2/ki_nspin1_dummy… done

Running calc_alpha/iteration_2/ki_nspin1… done

Running calc_alpha/iteration_2/ki_nspin2_dummy… done

Running calc_alpha/iteration_2/ki_nspin2… done

Orbital 10

———-

Running calc_alpha/iteration_2/orbital_10/pz_print… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1_dummy… done

Running calc_alpha/iteration_2/orbital_10/dft_n+1… done

alpha

1 2 3 4 … 7 8 9 10

0 0.600000 0.600000 0.600000 0.60000 … 0.600000 0.6000 0.600000 0.600000

1 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

2 0.655689 0.727571 0.783859 0.66386 … 0.729888 0.7419 0.779264 0.717389

[3 rows x 10 columns]

Delta E_i - epsilon_i (eV)

1 2 … 9 10

0 -4.773441e-01 -9.204971e-01 … -1.050159e+00 7.112091e-01

1 2.131628e-14 -2.415845e-13 … -4.440892e-14 1.776357e-14

[2 rows x 10 columns]

Screening parameters have been converged

Final KI calculation

====================

Running final/ki_final… done

Workflow complete

###  S4.2 Silicon

A typical input file for silicon is similar to that of ozone, but with some
tweaked workflow settings, a cell_parameters block that denotes that this
system is periodic, and some additional Wannierization settings:

⬇

{

”workflow”: {

”functional”: ”ki”,

”base_functional”: ”pbe”,

”method”: ”dscf”,

”init_orbitals”: ”mlwfs”,

”alpha_guess”: 0.1,

”eps_inf”: 13.02,

”orbital_groups_self_hartree_tol”: 0.1,

”pseudo_library”: ”sg15”,

”mp_correction”: false

},

”atoms”: {

”atomic_positions”: {

”units”: ”crystal”,

”positions”: [[”Si”, 0.00, 0.00, 0.00],

[”Si”, 0.25, 0.25, 0.25]]

},

”cell_parameters”: {

”periodic”: true,

”ibrav”: 2,

”celldms”: {”1”: 10.2622}

}

},

”kpoints”: {

”grid”: [4, 4, 4]

},

”calculator_parameters”: {

”ecutwfc”: 60.0,

”pw”: {

”nbnd”: 20

},

”w90”: {

”projections”: [

[{”fsite”: [0.25, 0.25, 0.25], ”ang_mtm”: ”sp3”}],

[{”fsite”: [0.25, 0.25, 0.25], ”ang_mtm”: ”sp3”}]],

”dis_froz_max”: 11.0,

”dis_win_max”: 16.5

},

”ui”: {

”smooth_int_factor”: 2

}

}

}

For a full explanation of the meaning of the Wannier90 keywords we refer the
reader to the Wannier90 documentation.

Running koopmans si.json generates the following output. This command prompts
a sequence of Quantum ESPRESSO calculations necessary to initialize the
density and variational orbitals using Wannier functions (lines 20-29), fold
these $k$-resolved functions to the equivalent $\Gamma$-only supercell (lines
31-36) calculate the screening parameters (lines 38-73), run a final KI
calculation (lines 75-77), and finally perform a second DFT Wannierization on
a finer $k$-grid to produce a smoothly interpolated band structure (lines
80-97).

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in si.bib in work involving this calculation

UserWarning: Makov-Payne corrections are not being used; do this with caution
for periodic systems

Initialization of density and variational orbitals

==================================================

Wannierization

==============

Running wannier/scf… done

Running wannier/nscf… done

Running wannier/block_1/wann_preproc… done

Running wannier/block_1/pw2wan… done

Running wannier/block_1/wann… done

Running wannier/block_2/wann_preproc… done

Running wannier/block_2/pw2wan… done

Running wannier/block_2/wann… done

Folding to supercell

——————–

Running block_1/w2kcp… done

Running block_2/w2kcp… done

Running init/dft_dummy… done

Running init/dft_init… done

Calculating screening parameters

================================

Running calc_alpha/ki… done

Orbitals 1-255

————–

Skipping; will use the screening parameter of an equivalent orbital

Orbital 256

———–

Running calc_alpha/orbital_256/dft_n-1… done

Orbital 257

———–

Running calc_alpha/orbital_257/pz_print… done

Running calc_alpha/orbital_257/dft_n+1_dummy… done

Running calc_alpha/orbital_257/dft_n+1… done

Orbitals 258-512

—————-

Skipping; will use the screening parameter of an equivalent orbital

alpha

1 2 3 … 510 511 512

0 0.100000 0.100000 0.100000 … 0.100000 0.100000 0.100000

1 0.104648 0.104648 0.104648 … 0.047709 0.047709 0.047709

[2 rows x 512 columns]

Delta E_i - epsilon_i (eV)

1 2 3 … 510 511 512

0 -0.023309 -0.023309 -0.023309 … -0.139624 -0.139624 -0.139624

[1 rows x 512 columns]

Screening parameters have been determined but are not necessarily converged

Final KI calculation

====================

Running final/ki_final… done

Postprocessing

===============

Wannierization

==============

Running wannier/scf… done

Running wannier/nscf… done

Running wannier/block_1/wann_preproc… done

Running wannier/block_1/pw2wan… done

Running wannier/block_1/wann… done

Running wannier/block_2/wann_preproc… done

Running wannier/block_2/pw2wan… done

Running wannier/block_2/wann… done

Running wannier/bands… done

UserWarning: Some of the pseudopotentials do not have PP_PSWFC blocks, which
means a projected DOS calculation is not possible. Skipping…

Running occ/ki… done

Running emp/ki… done

UserWarning: The DOS will not be plotted, because the Brillouin zone is too
poorly sampled for the specified value of smearing. In order to generate a
DOS, increase the k-point density (”kpath_density” in the ”setup” ”k_points”
subblock) and/or the smearing (”degauss” in the ”plot” block)

Workflow complete

...

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in si.bib in work involving this calculation

UserWarning: Makov-Payne corrections are not being used; do this with caution
for periodic systems

Initialization of density and variational orbitals

==================================================

Wannierization

==============

Running wannier/scf… done

Running wannier/nscf… done

Running wannier/block_1/wann_preproc… done

Running wannier/block_1/pw2wan… done

Running wannier/block_1/wann… done

Running wannier/block_2/wann_preproc… done

Running wannier/block_2/pw2wan… done

Running wannier/block_2/wann… done

Folding to supercell

——————–

Running block_1/w2kcp… done

Running block_2/w2kcp… done

Running init/dft_dummy… done

Running init/dft_init… done

Calculating screening parameters

================================

Running calc_alpha/ki… done

Orbitals 1-255

————–

Skipping; will use the screening parameter of an equivalent orbital

Orbital 256

———–

Running calc_alpha/orbital_256/dft_n-1… done

Orbital 257

———–

Running calc_alpha/orbital_257/pz_print… done

Running calc_alpha/orbital_257/dft_n+1_dummy… done

Running calc_alpha/orbital_257/dft_n+1… done

Orbitals 258-512

—————-

Skipping; will use the screening parameter of an equivalent orbital

alpha

1 2 3 … 510 511 512

0 0.100000 0.100000 0.100000 … 0.100000 0.100000 0.100000

1 0.104648 0.104648 0.104648 … 0.047709 0.047709 0.047709

[2 rows x 512 columns]

Delta E_i - epsilon_i (eV)

1 2 3 … 510 511 512

0 -0.023309 -0.023309 -0.023309 … -0.139624 -0.139624 -0.139624

[1 rows x 512 columns]

Screening parameters have been determined but are not necessarily converged

Final KI calculation

====================

Running final/ki_final… done

Postprocessing

===============

Wannierization

==============

Running wannier/scf… done

Running wannier/nscf… done

Running wannier/block_1/wann_preproc… done

Running wannier/block_1/pw2wan… done

Running wannier/block_1/wann… done

Running wannier/block_2/wann_preproc… done

Running wannier/block_2/pw2wan… done

Running wannier/block_2/wann… done

Running wannier/bands… done

UserWarning: Some of the pseudopotentials do not have PP_PSWFC blocks, which
means a projected DOS calculation is not possible. Skipping…

Running occ/ki… done

Running emp/ki… done

UserWarning: The DOS will not be plotted, because the Brillouin zone is too
poorly sampled for the specified value of smearing. In order to generate a
DOS, increase the k-point density (”kpath_density” in the ”setup” ”k_points”
subblock) and/or the smearing (”degauss” in the ”plot” block)

Workflow complete

###  S4.3 Zinc oxide

The input file is as follows:

⬇

{

”workflow”: {

”task”: ”singlepoint”,

”functional”: ”ki”,

”base_functional”: ”lda”,

”method”: ”dfpt”,

”init_orbitals”: ”mlwfs”,

”calculate_alpha” : true,

”pseudo_library”: ”pseudo_dojo_standard”,

”gb_correction” : true,

”eps_inf”: 5.3,

”from_scratch”: true,

”npool”: 4,

”dfpt_coarse_grid”: [2, 2, 2],

”orbital_groups_spread_tol”: 0.0005

},

”atoms”: {

”cell_parameters”: {

”periodic”: true,

”ibrav”: 4,

”celldms”: {”1”: 6.14057, ”3”: 1.60204}

},

”atomic_positions”: {

”units”: ”crystal”,

”positions”: [

[”Zn”, 0.33330, 0.66670, 0.50000 ],

[”Zn”, 0.66670, 0.33330, 0.00000 ],

[”O” , 0.33330, 0.66670, 0.11725 ],

[”O” , 0.66670, 0.33330, 0.61725 ]

]

}

},

”kpoints”: {

”grid”: [

4,

4,

4

],

”path”: ”ALMGAHK”

},

”calculator_parameters”: {

”ecutwfc”: 50.0,

”pw”: {

”system”: {

”nbnd”: 52

}

},

”w90”: {

”conv_tol”: 1e-8,

”projections”: [

[{”site”: ”Zn”, ”ang_mtm”: ”l=0”}],

[{”site”: ”Zn”, ”ang_mtm”: ”l=1”}],

[{”site”: ”O”, ”ang_mtm”: ”l=0”}],

[{”site”: ”Zn”, ”ang_mtm”: ”l=2”},

{”site”: ”O”, ”ang_mtm”: ”l=1”}],

[{”site”: ”Zn”, ”ang_mtm”: ”l=0”}]

],

”dis_froz_max”: 14.5,

”dis_win_max”: 17.0

}

},

”plotting”: {

”Emin”: -10,

”Emax”: 10

}

}

Here we can see the selection of DFPT for calculating screening parameters
(line 6), the choice of MLWFs as the variational orbitals (line 7), and a
criterion for grouping variational orbitals together based on their spreads
(line 15). We also have specified a coarse $2\times 2\times 2$ $k$-point grid
on which to calculate the screening parameters (line 14) relative to the
$4\times 4\times 4$ grid upon which we construct the Hamiltonian (specified
elsewhere in the output file). Again, it is worth stressing that these
calculations are not fully converged.

The output of koopmans zno.json is below. This command prompts a sequence of
Quantum ESPRESSO calculations which, using a coarse grid, calculates the
density and variational orbitals using Wannier functions (lines 20-40) and
then calculates the screening parameters using DFPT (lines 46-57). Then, on a
regular grid, it repeats a Wannierization (62-82) and then constructs and
diagonalizes the Koopmans Hamiltonian (lines 88-90), using the screening
parameters calculated on the coarse grid.

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in zno.bib in work involving this calculation

UserWarning: This system is not cubic and will therefore not have a uniform
dielectric tensor. However, the image-correction schemes that are currently
implemented assume a uniform dielectric. Proceed with caution

Coarse grid calculations

========================

Wannierization

==============

Running wannier/scf… done

Running wannier/nscf… done

Running wannier/block_1/wann_preproc… done

Running wannier/block_1/pw2wan… done

Running wannier/block_1/wann… done

Running wannier/block_2/wann_preproc… done

Running wannier/block_2/pw2wan… done

Running wannier/block_2/wann… done

Running wannier/block_3/wann_preproc… done

Running wannier/block_3/pw2wan… done

Running wannier/block_3/wann… done

Running wannier/block_4/wann_preproc… done

Running wannier/block_4/pw2wan… done

Running wannier/block_4/wann… done

Running wannier/block_5/wann_preproc… done

Running wannier/block_5/pw2wan… done

Running wannier/block_5/wann… done

Running wannier/bands… done

Running pdos/projwfc… done

Conversion to Koopmans format

—————————–

Running wannier/kc… done

Calculation of screening parameters

===================================

Running screening/band_2/kc… done

Running screening/band_6/kc… done

Running screening/band_8/kc… done

Running screening/band_10/kc… done

Running screening/band_16/kc… done

Running screening/band_18/kc… done

Running screening/band_20/kc… done

Running screening/band_24/kc… done

Running screening/band_26/kc… done

Running screening/band_27/kc… done

Regular grid calculations

=========================

Wannierization

==============

Running wannier/scf… done

Running wannier/nscf… done

Running wannier/block_1/wann_preproc… done

Running wannier/block_1/pw2wan… done

Running wannier/block_1/wann… done

Running wannier/block_2/wann_preproc… done

Running wannier/block_2/pw2wan… done

Running wannier/block_2/wann… done

Running wannier/block_3/wann_preproc… done

Running wannier/block_3/pw2wan… done

Running wannier/block_3/wann… done

Running wannier/block_4/wann_preproc… done

Running wannier/block_4/pw2wan… done

Running wannier/block_4/wann… done

Running wannier/block_5/wann_preproc… done

Running wannier/block_5/pw2wan… done

Running wannier/block_5/wann… done

Running wannier/bands… done

Running pdos/projwfc… done

Conversion to Koopmans format

—————————–

Running wannier/kc… done

Construction of the Hamiltonian

===============================

Running hamiltonian/kc… done

Workflow complete

...

⬇

_

| | _____ ___ _ __ _ __ ___ __ _ _ __ ___

| |/ / _ \ / _ \| ’_ \| ’_ ‘ _ \ / _‘ | ’_ \/ __|

| < (_) | (_) | |_) | | | | | | (_| | | | \__ \

|_|\_\___/ \___/| .__/|_| |_| |_|\__,_|_| |_|___/

|_|

Koopmans spectral functional calculations with Quantum ESPRESSO

version 1.0.0b6

Written by Edward Linscott, Riccardo De Gennaro, and Nicola Colonna

Please cite the papers listed in zno.bib in work involving this calculation

UserWarning: This system is not cubic and will therefore not have a uniform
dielectric tensor. However, the image-correction schemes that are currently
implemented assume a uniform dielectric. Proceed with caution

Coarse grid calculations

========================

Wannierization

==============

Running wannier/scf… done

Running wannier/nscf… done

Running wannier/block_1/wann_preproc… done

Running wannier/block_1/pw2wan… done

Running wannier/block_1/wann… done

Running wannier/block_2/wann_preproc… done

Running wannier/block_2/pw2wan… done

Running wannier/block_2/wann… done

Running wannier/block_3/wann_preproc… done

Running wannier/block_3/pw2wan… done

Running wannier/block_3/wann… done

Running wannier/block_4/wann_preproc… done

Running wannier/block_4/pw2wan… done

Running wannier/block_4/wann… done

Running wannier/block_5/wann_preproc… done

Running wannier/block_5/pw2wan… done

Running wannier/block_5/wann… done

Running wannier/bands… done

Running pdos/projwfc… done

Conversion to Koopmans format

—————————–

Running wannier/kc… done

Calculation of screening parameters

===================================

Running screening/band_2/kc… done

Running screening/band_6/kc… done

Running screening/band_8/kc… done

Running screening/band_10/kc… done

Running screening/band_16/kc… done

Running screening/band_18/kc… done

Running screening/band_20/kc… done

Running screening/band_24/kc… done

Running screening/band_26/kc… done

Running screening/band_27/kc… done

Regular grid calculations

=========================

Wannierization

==============

Running wannier/scf… done

Running wannier/nscf… done

Running wannier/block_1/wann_preproc… done

Running wannier/block_1/pw2wan… done

Running wannier/block_1/wann… done

Running wannier/block_2/wann_preproc… done

Running wannier/block_2/pw2wan… done

Running wannier/block_2/wann… done

Running wannier/block_3/wann_preproc… done

Running wannier/block_3/pw2wan… done

Running wannier/block_3/wann… done

Running wannier/block_4/wann_preproc… done

Running wannier/block_4/pw2wan… done

Running wannier/block_4/wann… done

Running wannier/block_5/wann_preproc… done

Running wannier/block_5/pw2wan… done

Running wannier/block_5/wann… done

Running wannier/bands… done

Running pdos/projwfc… done

Conversion to Koopmans format

—————————–

Running wannier/kc… done

Construction of the Hamiltonian

===============================

Running hamiltonian/kc… done

Workflow complete

##  S5 Details of the koopmans package

###  S5.1 Code structure

koopmans is built on top of the ASE python package (the Atomic Simulation
Environment) 100. Under the hood, it defines various Workflow classes, which
look like

⬇

class Workflow:

parameters: Dict[str, Any]

calculations: List[Calculator]

…

where the parameters attribute is a dictionary that stores the workflow
parameters as specified in the input file, and calculations is a list of the
calculations in the workflow. The individual entries in the calculations list
correspond to Calculator objects:

⬇

class Calculator(ASE_Calculator):

atoms: Atoms

parameters: Dict[str, Any]

results: Dict[str, Any]

…

which are subclasses of corresponding classes defined by ASE. ASE provides the
calculator with the functionality to read and write input and output files
(among many other things).

A Calculator object has — among others — an atoms attribute that stores the
details of the atoms and the simulation cell. The atoms attribute is itself an
instance of the Atoms class from ASE. We note that this hierarchy (namely,
that the atoms object is an attribute of a Calculator, and not the other way
around) is the reverse of the philosophy of ASE, where Atoms objects are the
principal object, and they may or may not have an associated calc attribute.

In addition to an atoms attribute, Calculator objects also have a parameters
attribute where calculator-specific settings are stored, as well as a results
attribute, where the results of the calculation are stored — just like in ASE.

###  S5.2 Scriptability

Because koopmans is written in python, integrating it within a script is
straightforward. For example, here is a script that runs the ozone calculation
from Section 4.1:

⬇

from ase import build

from koopmans.workflows import SinglepointWorkflow

# Create an Atoms object

atoms = build.molecule(’O3’, vacuum=5.0)

# Create a koopmans Workflow object

workflow = SinglepointWorkflow(atoms=atoms, ecutwfc = 65.0, nbnd = 10)

# Run the workflow

workflow.run()

# Fetch the IP and EA

results = workflow.calculations[-1].results

ip = -results[’homo_energy’]

ea = -results[’lumo_energy’]

# Print the IP and EA to screen

print(f’ IP = {ip:.2f} eV’)

print(f’ EA = {ea:.2f} eV’)

Of course, printing the IP and EA to screen is of limited value — in reality
at this stage the user would then generate plots, or feed these results to
another code.

Often a user will want to run workflows and analyse data separately — for
example, they might run their workflow on remote high performance computing
resources, and then, days later, analyse the results on their laptop. To
permit this, koopmans generates a .kwf file when a workflow is run. This file
can be loaded into python in order to recover the Workflow python object. For
example, we could perform exactly the same analysis on our previous ozone
calculation by replacing lines 1-11 with