A Kohn–Sham (KS) inversion determines a KS potential and orbitals corresponding to a given electron density, a procedure that has applications in developing and evaluating functionals used in density functional theory. Despite the utility of KS inversions, application of these methods among the research community is disproportionately small. We implement the KS inversion methods of Zhao–Morrison–Parr and Wu–Yang in a framework that simplifies analysis and conversion of the resulting potential in real-space. Fully documented Python scripts integrate with PySCF, a popular electronic structure prediction software, and Fortran alternatives are provided for computational hot spots.

## I. INTRODUCTION

Future progress in improving Kohn–Sham (KS) density functional theory (DFT), one of the most popular computational techniques for materials and molecules, depends upon improved exchange–correlation (XC) functionals.^{1} KS DFT assumes that non-interacting electrons and a local, multiplicative potential (i.e., KS potential) approximate the electron density of the real system’s interacting electrons. No systematic method for improving XC functionals has been identified within the KS system. Instead, functional development draws upon a broad array of techniques, methods, and data sources for guidance on how to produce more accurate approximations. One resourceful method is the KS inversion, which produces a KS potential from a provided electron density. Insight from KS inversions has been used in functional development since 1996^{2} and has seen a recent revival in its use today.^{3–5} It has also become a useful tool for studying DFT methods, such as time-dependent DFT,^{6,7} density-corrected DFT,^{8} inter-molecular interactions in partition-DFT,^{9,10} and embedded-DFT,^{11–13} and even in symmetry-adapted perturbation theory.^{14}

A KS inversion can construct an exact potential from an exact electron density as defined by the one-to-one density-to-potential mapping stated in the Hohenberg–Kohn theorem.^{15} However, in practice, the finite number of localized basis sets needed to expand KS orbitals destroys the advantageous one-to-one mapping,^{16–18} and the problem becomes ill-posed.^{18} Although an exact KS inversion is no longer possible, several approximate methods^{13,19–28} have been proposed.

Despite the developed theory and applications of KS inversions, publicly available software for routinely performing these calculations is not common.^{29} In this work, we introduce KS-pies, an open source implementation of the most frequently cited KS inversion methods, Zhao–Morrison–Parr^{23} (ZMP) and Wu–Yang^{24} (WY). Our WY implementation supports user-defined potential basis sets and Hamiltonians. This software includes a utility module that helps simplify inversion calculations by providing the input file conversion function and real-space evaluation of inversion potentials. The software is distributed under an open source Apache 2.0 Licence and developer community is hosted though GitHub.

Our Python implementation utilizes NumPy,^{30} SciPy libraries,^{31} and features from PySCF^{32} that are familiar to the DFT community. A section of Python code can alternatively call kspies_fort, a compiled Fortran version that decreases computational cost.

In this paper, we include a theoretical summary, details on implementation, validation, and an Appendix with a discussion based user-guide connecting technical features with examples that highlight the simplicity of incorporating KS inversions into modern functional development workflows. Emphasis on the theoretical background seeks to provide a simplified and practical entry point for theorists unfamiliar with the methodology.

## II. BACKGROUND

### A. Kohn–Sham density functional theory and its inverse

The conventional (forward) KS procedure solves a single particle equation

where *ψ*_{i} and *ε*_{i} are *i*th KS orbitals and orbital energies, *v*_{S}[*n*](**r**) is the KS potential, and *n* is the electron density. Atomic units are used throughout unless specified. Typically, *v*_{S}[*n*](**r**) is written as

where *v*_{ext}(**r**) is the external potential, *v*_{H}[*n*](**r**) is the Hartree potential,

and *v*_{XC}[*n*](**r**) is an XC potential, which is approximate in practice. The electron density is determined by the sum of the occupied orbital densities as

As a functional of the density, *v*_{S}[*n*](**r**) is used in Eq. (1) in a self-consistent field (SCF) procedure, until a converged electron density is determined.

KS inversions operate in reverse using a given density [often referred as target density, *n*^{tar}(**r**)] to determine *v*_{S}. Once determined, *v*_{S} and Eq. (1) produce the KS orbitals and associated eigenvalues. In principle, one-to-one density-to-potential mapping stated in the Hohenberg–Kohn theorem^{15} guarantees that KS inversion can construct an exact KS potential (up to a constant) from an exact electron density. However, in practice, the finite number of localized basis sets needed to expand KS orbitals destroys the advantageous one-to-one mapping,^{16–18} and the problem becomes ill-posed.^{18} Although an exact KS inversion is no longer possible, several approximate methods^{13,19–28} have been proposed. KS inversions have been developed for use with input orbitals^{13,33} or wavefunctions,^{34} but we focus exclusively on density-based methods,^{19–28} including ZMP^{23} and WY,^{24} which will be explained in the following.

### B. Zhao–Morrison–Parr

The ZMP KS inversion method^{23} minimizes an objective self-repulsion functional,

by solving a KS-like equation self-consistently under a given Lagrange multiplier *λ*,

A *λ*-dependent KS potential

includes a guiding potential *v*_{g}(**r**) and a correction potential

In principle, as *λ* → ∞, *C*[*n*^{λ}] → 0 and *n*^{λ} → *n*^{tar}. Only $vC\lambda $ depends on *λ* in Eq. (7) and accommodates all necessary potential modification. If provided, additional potential terms *v*_{ext}(**r**), *v*_{H}(**r**), and *v*_{g}(**r**) accelerate the convergence of *v*_{S} with respect to *λ*. In practice, Eq. (6) is solved self-consistently using a given *λ* value. The obtained orbitals are used as an initial guess for following calculations at larger *λ* values, a process repeated until *λ* become large enough.

The guiding potential *v*_{g}(**r**) mimics the XC potential, for which a variety of potentials can be used. Typically, *v*_{g}(**r**) is initially formulated to mimic the asymptotic decay of XC potential, −(1/*N*)*v*_{H}(**r**), where *N* is the number of electrons. We refer to −(1/*N*)*v*_{H}(**r**) as FAXC, the non-Hartree portion of the Fermi–Amaldi potential.^{35} In principle, any potential can be used for *v*_{g}(**r**) when the asymptotic decay of XC potential is not important.^{8}

For open-shell systems, ZMP is used as a spin-unrestricted formalism.^{36} The correction potential is spin-dependent, and Eq. (8) is rewritten for *α* or *β* spin as

where *σ* denotes the spin index. The factor of 2 is required for consistent results in closed-shell systems for restricted and unrestricted schemes. The guiding potential is spin-dependent for standard DFT XC potentials, but not for FAXC.

### C. Wu–Yang

The WY approach^{24} maximizes an objective functional,

where $vSb(r)$ is a KS potential similar to Eq. (7) in ZMP, except with *v*_{C}(**r**) represented as a linear combination of potential basis functions *g*_{t}(**r**), given as

making $vSb(r)$ exclusively a functional of the target density. KS orbitals $\psi ib(r)$ are determined by solving a KS-like equation,

which does not require the SCF procedure. The objective functional *W*_{S} is maximized by adjusting {*b*_{t}}. The gradient and Hessian of *W*_{S} with respect to {*b*_{t}} are given in an analytical form

where *N*_{occ} (*N*_{vir}) denotes the number of occupied (virtual) orbitals, simplifying maximization of Eq. (10).

In a spin-unrestricted formalism, Eq. (10) becomes

Spin dependence of XC potentials is reflected in the correction (*v*_{C}) potential due to its spin dependence on {*b*_{t}} and standard DFT XC guiding (*v*_{g}) potentials.

Highly oscillatory XC potentials are a well-known drawback of the WY method that often arises when using a large potential basis.^{37} One solution to this problem is regularization with the objective functional

where *η* is a regularization strength hyperparameter, and the smoothness of the correction potential is measured with

An optimal value of *η* must be selected. If *η* is too small, the potential may contain severe oscillations, but if *η* is too large, the optimized potential becomes too smooth and misses physically important features. Heaton-Burgess and Yang^{38} describe how this optimal selection can be made.

## III. IMPLEMENTATION

### A. Architecture and general workflow of KS-pies

KS-pies consists of three sub-modules; zmp and wy, which are the modules for KS inversion, and util module that supports additional utility functions to help prepare inputs to zmp and wy or analyze inversion results. To start inversion calculation (ZMP or WY) with KS-pies, two inputs are required: Mole object and density matrix of the target density. Mole object defines standard details needed for quantum chemical calculations, such as atomic coordinates, number of electrons, and basis sets. This is exactly the Mole object defined in PySCF, which can be defined with only a few lines of code or loaded from a common geometry file such as xyz. The second input, density matrix, can be generated with PySCF calculations or loaded from an external file format such as molden.^{39} To extend accessibility to KS-pies and PySCF, the util sub-module in KS-pies also supports loading the wfn file format, which is supported for variety of quantum chemistry programs, such as Gaussian,^{40} ORCA,^{41} GAMESS,^{42} and Molpro.^{43} Target densities generated with quantum chemistry programs that support molden or wfn formats can be used as inputs for KS-pies.

Performing an inversion calculation with KS-pies is very similar to running KS-DFT calculation with PySCF. Users can specify algorithm-dependent options for running ZMP and WY calculations. All inversion calculations are defined as a class object, and KS-pies manages data using Python instance variables, storing the majority of results in memory. Data from one instance can be used as an initial guess for subsequent calculations or analysis within KS-pies or PySCF. util module supports additional utility functions to help prepare inputs to zmp and wy and analyze inversion results.

KS-pies makes use of analytical functions within the PySCF integral library^{44} to convert Eqs. (6) and (12) to solvable matrix equations. The Hartree and guiding potentials are constructed from the target density using a matrix representation by default. KS-pies then uses a method-specific procedure to optimize the correction potential. For ZMP, a self-consistent calculation with Eq. (6) and a user provided *λ* is performed. For WY, the obtained gradient [Eq. (13)] and Hessian [Eq. (14)], {*b*_{t}} are adjusted to maximize Eq. (10) using the SciPy optimizer.

WY can accommodate non-Gaussian potential basis sets. Although a Gaussian function is typically used to expand potentials due to its integration efficiency, Eq. (11) is not limited to Gaussian type basis set. KS-pies accommodates this uncommon feature and can handle the user-defined potential basis set. When arbitrary user-defined potential basis sets are encountered, a numerical integration of the three-center overlap integral used in Eq. (10) is calculated with

where *ϕ*_{i}(**r**) is *i*th orbital basis function. Numerical integration of Eq. (18) adds a substantial computational overhead at the beginning of the calculation.

Several features reduce the computational cost of determining the correction potential. For ZMP, recalculating the Hartree potential at each self-consistent iteration can be accelerated using a density fitting procedure. Density fitting results in minor differences in the inversion while greatly reducing the computational cost. For WY, Eqs. (13) and (14) are calculated using a Fortran module kspies_fort to provide a substantial time savings over a Python implementation. Compiled binaries and the source code are provided.

### B. Real-space potentials

The KS potential in real-space provides valuable insight beyond a target density calculation. For example, a major motivation for performing a KS inversion is the visualization of the exact KS potential that can show the non-intuitive step structures present in some potentials.^{45,46} However, the use of finite atom-centered Gaussian basis sets prevents the XC potentials produced by ZMP and WY from being directly converted into real-space representation. For a given system, the basis functions {*ϕ*_{i}(**r**)}, real-space function *v*(**r**), and its matrix representation *V*_{ij} have the relation

However, Eq. (19b) is only exact under the basis set limit and, in practice, would result in large errors in the real-space function *v*(**r**). Therefore, a method other than Eq. (19b) is required to obtain real-space values of *v*(**r**). This is necessary when the guiding or correction potential needs to be evaluated in real-space. The method of Franchini *et al.*^{47} for converting Hartree potentials can be applied to the FAXC guiding potential and the ZMP correlation potential. A DFT XC potential on real-space is necessary when ZMP or WY utilizes a DFT XC guiding potential.

KS-pies evaluates real-space Hartree potentials following Franchini *et al.*^{47} In this approach, the density is decomposed into one-center (i.e., atomic) contributions and into different angular contributions as

where *Z* represents real spherical harmonics and *s* is a cubic spline interpolation of radial density at the *i*th atom radial grid. The summation inside Eq. (20) can be used to calculate the Hartree potential for each angular contribution of the *i*th atom with the analytical form

This allows conversion of a density to a real-space Hartree potential directly, without going through matrix representation, and is implemented in kspies.util.eval_vH.

The kspies.util.eval_vxc function evaluates the DFT XC potential on user-defined grid points. Our implementation is based on the numerical differentiation. For local density approximation (LDA) functionals, we use

where $\epsilon XCLDA$ is the XC density of LDA and can be directly obtained from pyscf.dft.libxc.eval_xc. For generalized gradient approximation (GGA) functionals, we use

where $vn=\u2202\epsilon XCGGA/\u2202n$, $v\gamma =\u2202\epsilon XCGGA/\u2202\gamma $, and *γ* = ∇*n*·∇*n*. Although *v*_{n} and *v*_{γ} are obtainable from pyscf.dft.libxc.eval_xc, ∇*v*_{γ} should be evaluated using a numerical derivative of *v*_{γ}. For spin-polarized densities, Eq. (23) for *α* spin becomes

## IV. VALIDATION AND PERFORMANCE

Valid implementation is confirmed with accurate KS inversions of densities obtained from HF or correlated wavefunction methods, in both restricted and unrestricted schemes. Run time benchmarks are reported as wall time and were performed with an Intel® Xeon® Gold 6142 CPU using eight processors at 2.6 GHz. Calls to PySCF and kspies_fort can take advantage of parallelization with OpenMP. In our KS inversion examples below, benzene used the most memory, ∼1.2 GB. The PySCF coupled-cluster singles-and-doubles (CCSD) calculation of molecular O_{2} used substantially more memory; however, the inversion using KS-pies required less than that of benzene. Restricted and unrestricted inversion benchmarks are included for convenience within the software repository.

### A. Restricted inversion

Restricted inversion performance was validated using ZMP and WY on benzene (*R*_{CC} = 1.3936 Å and *R*_{CH} = 1.0852 Å). The target density was generated with HF/cc-pVTZ. Resulting inversion potentials from ZMP and WY can be used to accurately reproduce the target density. Furthermore, ZMP qualitatively indicates that *C* [Eq. (8)] is approaching 0 as *λ* increases.

kspies.zmp.RZMP (restricted ZMP) was used with a FAXC guiding potential and tested with and without density fitting. Using the self-repulsion functional value *C* from Eq. (5) and the integrated density difference *dN* = ∫|*n*^{λ}(**r**) − *n*^{tar}(**r**)|*d***r**, we confirmed an expected decrease in *C* as *λ* increased. This is plotted in Fig. 1(a). As *λ* increases, *dN* also decreases. Our implementation shows negligible difference in accuracy with and without density fitting.

Computationally, the density fitting method decreases the cost of ZMP. Figure 1(b) plots computational performance in terms of SCF iterations and run time per *λ*, highlighting negligible difference in the number of SCF cycles. Per SCF iteration, ZMP takes ∼0.25 s with density fitting and 3 s without, requiring around 700 iterations to reach convergence on small molecules.

WY performance was evaluated using the same target density as above and cc-pVTZ for the potential basis set. Optimization was complete after eight iterations, taking less than 1 s on default settings. The maximum gradient element was 3 × 10^{−8} with *dN* = 170.8 me. The *dN* agrees with the values from ZMP at *λ* = 128, indicating an accurately reproduced target density. Implementation verification is also incidentally discussed in Subsection 5 of the Appendix where WY reproduces a target density from a user-defined harmonic potential Hamiltonian.

### B. Unrestricted inversion

To validate the unrestricted calculations of ZMP and WY, we use a coupled-cluster singles-and-doubles (CCSD) target density of molecular oxygen (*R*_{OO} = 1.208 Å) obtained with UHF-UCCSD/cc-pVQZ. We used the FAXC guiding potential for both cases.

Benchmark values for ZMP with *λ* = 2048 are *C* = 1.10 × 10^{−6} and *dN* = 5.75 me. A small and decreasing *dN* as *λ* is increased confirms the unrestricted ZMP implementation. As a secondary validation benchmark, identical results are produced from spin-restricted and unrestricted ZMP calculation on closed-shell benzene.

Benchmark values for WY with a cc-pVQZ potential basis set converges after five optimization steps, taking 0.07 s, with a maximum gradient element of 3 × 10^{−8}, and dN = 36.3 me, which agrees with ZMP *dN* at *λ* = 128, verifying the WY implementation.

### C. Inversion of user-defined systems

WY inversions can be calculated on user-defined systems using kspies.wy and a properly defined, user supplied, Hamiltonian component. Figure 2 shows a calculated potential from a user-defined Hamiltonian of a harmonic potential [*v*_{ext}(*x*) = 1/8*x*^{2}, blue solid curve in Fig. 2] with four electrons. Finite-difference HF with soft electron–electron repulsion $w(x1,x2)=1/(x2\u2212x1)2+0.52$ was used to generate the target density within the domain *x* ∈ [−10, 10] and a grid spacing of 0.02. The WY calculation was then performed using the finite-difference method by providing a finite-difference Hamiltonian and using a grid to expand the potential.

Subsection 5 of the Appendix and the online KS-pies documentation expand upon the present example with information on how to obtain necessary inputs for user-defined Hamiltonians.

The resulting WY KS potential for the user-defined harmonic potential is displayed in orange in Fig. 2. The KS potential includes small oscillations near *x* = ±7 due to a known issue of WY^{18} arising from the small electron density in these regions. The oscillations are not a fault of our implementation. The resulting potential and its oscillations are not unique.^{48} In the general use of WY, these oscillations can be influenced by conditions, such as the optimization scheme, convergence criteria, or selection of a basis set. Nevertheless, the agreement of the density as represented by the red and black curves in Fig. 2 highlights the accuracy of WY theoretically, as well as our successful implementation.

## V. USAGE

To encourage the use of KS inversions in functional development, kspies requires only a few inputs. In the Appendix, we provide several examples that highlight the relative simplicity of routine KS inversions with and without real-space conversions. The ZMP section details basic ZMP use, and the WY section details basic use as well as an example using regularization and another using a user-defined harmonic potential. Utility module examples demonstrate importing target densities from other software and analysis tools for evaluating and visualizing resulting KS potentials.

The package can be easily installed from PyPI using pip install kspies, and full documentation, examples from the Appendix, and test scripts are available online. Frequent users can compile a provide Fortran subroutine that improves the speed of WY calculations.

## VI. CONCLUSION

KS-pies presents an open source, publicly available code for performing KS inversions of electron densities into KS potentials. Since every KS inversion method is approximate, we implemented the two most cited inversion methods, ZMP and WY. Our software integrates with PySCF, an environment familiar to the theoretical development community. Our framework provides a starting point for the implementation of future KS inversion methods. This publication presents the theoretical context and examples that highlight the simplicity of running KS inversions. Incorporating KS inversion methods to determine real-space potentials should be beneficial for XC functional development and testing.

With two implemented methods, users are able to choose and compare results, leveraging advantages of each method. ZMP requires many SCF iterations and is computationally intensive relative to WY, but the result of inversion can be systematically improved by increasing *λ*. Alternatively, WY can perform inversions on user-defined Hamiltonians and is computationally efficient. Potentials produced by both methods can be converted into real-space representations using our software.

## ACKNOWLEDGMENTS

We thank Kieron Burke for thoughtful manuscript feedback. S.N., H.P., and E.S. are thankful for support from the National Research Foundation of Korea (Grant Nos. NRF-2020R1A2C2007468 and NRF-2020R1A4A1017737). R.J.M is thankful for support from the University of California President’s Postdoctoral Fellowship and the National Science Foundation (Grant No. CHE 1856165).

## DATA AVAILABILITY

The KS-pies code is openly available in GitHub at https://github.com/ssnam92/KSPies and can be referenced via https://doi.org/10.25351/V3.KS-PIES.2020.

### APPENDIX: BASIC KS-PIES USAGE AND EXAMPLES

#### 1. Zhao–Morrison–Parr method

ZMP calculations proceed by instantiating a kspies.zmp.RZMP object and then calling zscf(l), which generates a KS potential with the user specified l (*λ*) value. The instance requires a Mole() object that defines the basis set and an atomic orbital representation of the target density. Calculation results are printed to the terminal, and a matrix representation form of the potential is stored as instance attributes. This can be converted into real-space representation with the later discussed kspies.util module.

Figure 3 demonstrates generating a Ne HF density using PySCF, generating a KS potential with kspies.zmp, and the run outputs printed to the terminal. Creating a RZMP instance requires a Mole() object (mol_1) and density matrix of the target density (P_tar_1). In the output, niter is the number of SCF iterations needed to reach convergence, gap is the HOMO–LUMO gap in atomic unit, dN is the integrated density difference in millielectrons, and C is the minimized value from Eq. (5). In unrestricted ZMP, resulting C values account for contributions from both spins, 2(*C*_{α} + *C*_{β}).

Additional options in ZMP include the following:

diis_space (int): DIIS space size. Default is 40.

max_cycle (int): Maximum SCF iterations. Default is 400.

guide (character): Guiding potential. None sets

*v*_{g}(**r**) = −*v*_{H}(**r**) [i.e., only the external and correction potential cover*v*_{S}see Eq. (7)], “faxc” sets*v*_{g}(**r**) = −*v*_{H}(**r**)/*N*, or any DFT XC functional defined in PySCF can be specified at any additive value. For example, setting ‘b3lyp-0.2^{*}hf+0.2^{*}faxc’ will remove the exact exchange portion from the b3lyp potential and add or replace it with the faxc guiding potential.level_shift (float): Amount of level shift used during SCF cycles. Default is 0.2.

with_df (boolean): Use of density fitting. Default is False.

conv_tol_dm (float): Density matrix convergence criteria. Default is 1 × 10

^{−7}.conv_tol_dm (float): Convergence criteria of DIIS error. Default is 1 × 10

^{−4}.

The results from zscf(l) are stored as the following class attributes:

converged (boolean): If convergence criteria was met during SCF iterations;

dm (array): Density matrix;

mo_coeff (array): Molecular orbital coefficient;

mo_energy (array): Molecular orbital energy; and

mo_occ (array): Orbital occupation numbers.

For kspies.zmp.UZMP, alpha and beta values for dm, mo_coeff, mo_energy, and mo_occ are stored as an ordered pair tuple. The instance only stores one set of results; previous values will be overwritten, and only the most recent calculations results are accessible.

On the first zscf() SCF cycle, the target density is used as the initial density matrix guess, as demonstrated Fig. 4 with P_tar_1. The guiding potential is constructed from user specifications during the first call to the function. Subsequent changes to the user specifications for the guiding potential after the first call are ignored.

The direct inversion of the iterative subspace (DIIS) procedure^{49} for convergence stabilization implemented in ks-pies is independent of the PySCF DIIS options. We recommend setting diis_space ≤40 to avoid a matrix singularity during DIIS extrapolation. DIIS convergence stabilization is not used when diis_space ≤1. ZMP often fails to converge without DIIS, and we recommend using this feature, even for small *λ*.

When specified, virtual orbital energies are increased by level_shift, which can aid convergence. In general, setting level_shift = 0.1 × *λ* is likely to be sufficient for most systems. level_shift is inactive when set to 0. When *λ* is larger than 10–20, level_shift values of 1–2 or larger are needed for convergence. The larger values for level_shift reduce the mixing between the occupied and virtual orbits, which slows the orbital rotation with the intent of assisting convergence. In some systems, the initial ZMP iterations will significantly perturb the potential from a trajectory toward convergence. Increasing level_shift to slow the orbital rotation as shown in Fig. 4 can minimize or remove this effect. As level_shift is an artificial value added to aid in convergence, its contribution is ignored when calculating properties such as the HOMO–LUMO gap and orbital energies.

#### 2. Wu–Yang method

Performing a KS inversion with kspies.wy requires an input Mole() object and target density, as demonstrated in Fig. 5. The potential basis (pbas) defaults to the atomic orbital basis, and Sijt is integrated analytically. At the conclusion of Fig. 5, the KS potential produced by maximizing *W*_{S} is now accessible as an instance attribute.

A second WY example in Fig. 6 demonstrates additional options and uses non-equivalent basis functions. The user specifies the potential basis set with pbas, demonstrated in the example with aug-cc-pV5Z. The terminal output (Fig. 6, bottom) shows the HOMO–LUMO gap and the smoothness of the correction potential [Eq. (17)] for each eta.

Figure 7 calls the instance from Fig. 6 for information on the run, the number of basis functions used in the potential, and confirmation on convergence. The instance overwrites itself after each calculation, and the reported values are for eta = 1 × 10^{−6}.

Options available in kspies.wy include the following:

method (string): Optimization algorithm used in scipy.optimize.minimize. Default is “trust-exact.”

guide (string): Guiding potential. Default is “faxc.” Usage is the same as guide in ZMP.

tol (float): Tolerance of the maximum gradient values used to determine optimization completion [Eq. (13)]. Default is 1 × 10

^{−6}.reg (float): Potential regularization weight [

*η*in Eq. (16)]. Default is 0.

The molecular orbital coefficients, energies, occupation numbers, density matrices, and convergence are all stored as instance attributes. Optimized {*b*_{t}} values are stored as “b” and can be useful as an initial guess for subsequent calculations using different regularization strength *η*, as demonstrated in the for loop of Fig. 6. Other notable methods include make_rdm1() for generating a density matrix in the same format and method as the PySCF function of the same name, info() for printing optimization results, and Dvb() for calculating the smoothness of the correction potential [Eq. (17)].

A WY inversion may fail if the Scipy optimizer failed to find a maximum *W*_{S}. Common reasons for this failure with straightforward solutions include the following: 1. The given electron density is not *v*_{S}-representable. As the target density cannot be reconstructed with a density from single-determinant. Schipper *et al.*^{50} provide further discussion about *v*_{S}-representability. 2. A bad initial guess was used. By default, {*b*_{t}} is initialized as zero. Users can specify an initial {*b*_{t}} guess with the b attribute, i.e., wy_a.b = b_init. For stretched molecules, an initial guess {*b*_{t}} from less stretched calculations may solve this failure. 3. A nearly singular Hessian was encountered. This can occur when the potential basis is very large. We recommend a gradient-based optimization algorithm for these problems, such as conjugate gradient (CG) or Broyden–Fletcher–Goldfarb–Shanno (BFGS). 4. An unsuitable guiding potential was used. For example, using a semi-local DFT XC guiding potential for the stretched ionic bond system would lead to failure in a similar nature to how a forward KS calculation using the same functional may have trouble identifying a converged density.^{8}

#### 3. Regularized WY

Regularized WY can compensate for issues resulting from unbalanced potential basis sets.^{37} Figure 8 provides example inputs on the a N_{2} molecule, and Fig. 9 displays example plots used to select the regularization parameter. kspies.wy.RWY requires a Mole() object, a target density, and a large even-tempered Gaussian basis set for the potential calculation, all of which can be produced using PySCF.

The optimal regularization parameter, *η*, can be identified using the L-curve method of Ref. 51, where the maximum value of the reciprocal derivative should be selected. As *η* decreases, $W\xafS\eta $ approaches *W*_{S} and the curvature of the KS potential increases (i.e., increasing $\Vert \u2207vCb\Vert 2$). The need for regularization arises because an unbalanced potential basis set may cause $\Vert \u2207vCb\Vert 2$ to increases sharply but is not guaranteed to correspond with a significant decreases in $WS\u2212W\xafS\eta $ [see Fig. 9(a)]. This implies that the increase in potential curvature results from unphysical oscillation. According to Ref. 37, an optimal *η* gives a maximum reciprocal slope of the L-curve, which has the analytical form

#### 4. WY for user-defined potential basis

Although Gaussian functions are typically used to expand correction potential in WY due to their integration efficiency, this basis choice is not mandatory. KS-pies can perform WY calculation on any user-defined potential basis. Figure 10 shows the example of using Slater-type basis functions to expand the WY correction potential. A user can designate pbas as a list of functions that take Bohr-unit xyz coordinates as inputs and return a value for each coordinates. Using this pbas, kspies.wy will then numerically calculate the three-center overlap integral and perform a WY calculation with it. Note that numerical integration of the three-center overlap integral adds an initial computing cost, which can be substantial for large systems.

#### 5. WY for user-defined systems

KS-pies can calculate a WY KS potential for user-defined Hamiltonians in the restricted and unrestricted formalism, as demonstrated in Fig. 11. A user must provide a Mole() object (mol), target density density matrix (P_tar), and a Sijt array to instantiate a user-defined system, as well as the following instance attributes: kinetic energy T, external potential V, overlap matrix S, kinetic energy matrix of the potential basis Tp, and a None override to prevent default use of a guiding potential.

#### 6. Utility for cross-platform inputs

The kspies.util.wfnreader function loads wfn file formats as generated by many quantum chemistry packages, such as Gaussian,^{40} ORCA,^{41} GAMESS,^{42} and Molpro,^{43} and converts it into the PySCF format used in kspies. Figure 12 shows an example usage of kspies.util.wfnreader to read wavefunction information from a ccsd.wfn file. Note that system information such as nuclear coordinates, number of electrons, and basis sets stored in wfn format should be the same as in mol. The dm_tar output in the second line of Fig. 12 is ready to serve as a target density for ZMP or WY.

Of course, users also can use PySCF default conversion tools to load densities that are generated with other platforms, which can extend the usability of KS-pies.

#### 7. Evaluation utility

The usage of the kspies.util.eval_vh function for evaluating XC potentials [Fig. 15(a)] is presented in Fig. 13. Referring to the example, kspies.util.eval_vh requires a Mole() object (mol), density matrix to calculate Hartree potential (dmxc), and specified xyz-coordinates (coords) that describe the positions of grid points that are to be calculated. At a given *λ*, the ZMP XC potential can be written as

indicating that the XC potential obtained from ZMP at the specific *λ* is the Hartree potential of the density *λn*^{λ} − ((1/*N*) + *λ*)*n*^{tar}.

An example of kspies.util.eval_vxc is used in Fig. 14 to create the visualization of the WY XC potential in Fig. 15(b). The finite difference required for numerical differentiation of *v*_{γ} in Eqs. (23) and (24) is set with delta (atomic units) in eval_vxc and defaults to 1 × 10^{−7}.

Either ZMP or WY methods can be used with kspies.utils.eval_vh and kspies.utils.eval_vxc. Some useful possibilities beyond our current examples include using ZMP with PBE XC guiding potential to draw PBE XC potential with kspies.utils.eval_vxc and ZMP correction potential with kspies.utils.eval_vh. Visualization with kspies.utils is not limited to XC potential obtained from KS inversion but can be used independently with KS inversion.