koopmans: An Open-Source Package for Accurately and Efficiently Predicting Spectral Properties with Koopmans Functionals

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 semilocal 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.


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 If we take the derivative with respect to the occupancy of the j th variational orbital then we have where we assumed that the cross-term derivatives dΠ i /df j vanish, and because E Koopmans ought to be linear in f j , we replaced its derivative with some yet-to-be deterimined constant η j .For the second equality we invoked Janak's theorem, and f is some number between 0 and 1.
Assuming that the energy correction Π j is zero at integer occupancies, is independent of f i for i = j, and neglecting for the moment any orbital relaxation as the orbital occupancies change, it follows that 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 {α i } such that Π j = α j Π u j .Having done this, we arrive at eq. 5 of the main text, the final result.

S2 KIPZ details
In previous works, KIPZ has been presented in slightly different ways.In eq.27 of Ref. S1, KIPZ was introduced as )   where ĤPZ i (s) = ĤDFT (s) − vDFT  S5)   and in that same paper it was also stated that One can prove that these three definitions are equivalent via the identity from which it follows that 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, S3-S5 but it would be interesting to explore alternative prescriptions for the scaling of the PZ correction that are decoupled from the generalized PWL condition. S-4

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.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 γ → 0 limit of the "KIγ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 Γ-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.S6 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.xcodes (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.S7 The one final task is to map the Wannier functions in the k-sampled primitive cell to the equivalent Γ-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 eqs. 10 and 11 from the main text, with which we calculate these parameters: for occupied orbitals and 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 α 0 (node g).This gives us access to the energy of the N -electron system S-7 E(N ) (node m) as well as all of the requisite expectation values Koopmans Hamiltonian on the variational orbitals λ ii (α, f ) (nodes k, l, n, and o).The second and third calculations (nodes h and i) are calculations on the N ± 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 ± 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.S8 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, ∆E i is independent of α, as are occupied variational orbitals, and consequently λ KI ii (α, 1) is linear in α.This is not the case for empty orbitals for the KI functional (for which λ KI ii (α, 0) is not strictly linear) or for the KIPZ functional (where additionally ∆E i is dependent on α).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 Γ-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. S6.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 The Koopmans potential is very smooth and slowly-varying in k-space, applying an almostconstant 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 k-indexed Hamitonian as where now {R } corresponds to a much larger supercell or, equivalently, a much denser kpoint 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. S6.

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.

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 eq.21 (from the main text) are then performed by a single kcw.x run (node x).These calculations yield the screening parameters (node y).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)

S-21
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: 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.

S-22
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 .kwffile 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 from koopmans import io workflow = io.read('ozone.kwf') where ozone.kwfhas been generated by some previously completed koopmans calculation.
In the above, we used the SinglepointWorkflow for running a Koopmans workflow from start to finish.koopmans implements several other workflows that automate tasks that are useful when performing Koopmans calculations, such as convergence testing, standalone Wannierization, and DFT calculations.

S5.3 Code quality and testing
koopmans contains an extensive test suite implemented with pytest.S10 It also has typing annotations which allow it to be statically typechecked using mypy.
Hxc [s|ϕ i (r)| 2 ]. (Ref.S1 included an erroneous sum over i in the definition of this Hamiltonian.)Meanwhile, eq.6 of Ref S2 defined KIPZ as 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.

Figure S2 :
Figure S2: The DFPT workflow.Individual nodes are explained in the text.
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).

##
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) Print the IP and EA to screen print(f' IP = {ip:.2f}eV') print(f' EA = {ea:.2f}eV') 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: where