KS-pies: Kohn–Sham inversion toolkit

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.¹ 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² 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,⁸ inter-molecular interactions in partition-DFT,^9,10 and embedded-DFT,^11–13 and even in symmetry-adapted perturbation theory.¹⁴ 

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.¹⁵ 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.¹⁸ 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.²⁹ In this work, we introduce KS-pies, an open source implementation of the most frequently cited KS inversion methods, Zhao–Morrison–Parr²³ (ZMP) and Wu–Yang²⁴ (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,³⁰ SciPy libraries,³¹ and features from PySCF³² 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.

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

where ψ_i and ε_i are ith 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¹⁵ 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.¹⁸ 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,³⁴ but we focus exclusively on density-based methods,^19–28 including ZMP²³ and WY,²⁴ which will be explained in the following.

The ZMP KS inversion method²³ 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λ$ 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.³⁵ In principle, any potential can be used for v_g(r) when the asymptotic decay of XC potential is not important.⁸ 

For open-shell systems, ZMP is used as a spin-unrestricted formalism.³⁶ 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.

The WY approach²⁴ 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 $ψ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.³⁷ 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³⁸ describe how this optimal selection can be made.

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.³⁹ 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,⁴⁰ ORCA,⁴¹ GAMESS,⁴² and Molpro.⁴³ 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⁴⁴ 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 ith 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.

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.⁴⁷ 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.⁴⁷ 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 ith atom radial grid. The summation inside Eq. (20) can be used to calculate the Hartree potential for each angular contribution of the ith 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 $ε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=∂εXCGGA/∂n$⁠, $vγ=∂εXCGGA/∂γ$⁠, 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

and the formulation for β spin requires a trivial swapping of respective spins. Equations (23) and (24) are implemented in kspies.util.eval_vxc. Subsection 7 of the  Appendix and associated figures, Figs. 13-15, provide an example of obtaining real-space representation using these approaches.

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

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)|dr, 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⁻⁸ 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.

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⁻⁶ 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⁻⁸, and dN = 36.3 me, which agrees with ZMP dN at λ = 128, verifying the WY implementation.

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/8x², blue solid curve in Fig. 2] with four electrons. Finite-difference HF with soft electron–electron repulsion $w(x1,x2)=1/(x2−x1)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¹⁸ 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.⁴⁸ 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.

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.

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.

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

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.

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⁻⁷.

• conv_tol_dm (float): Convergence criteria of DIIS error. Default is 1 × 10⁻⁴.

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⁴⁹ 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.

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⁻⁶.

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⁻⁶.

• 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.⁵⁰ 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.⁸ 

Regularized WY can compensate for issues resulting from unbalanced potential basis sets.³⁷ Figure 8 provides example inputs on the a N₂ 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¯Sη$ approaches W_S and the curvature of the KS potential increases (i.e., increasing $‖∇vCb‖2$⁠). The need for regularization arises because an unbalanced potential basis set may cause $‖∇vCb‖2$ to increases sharply but is not guaranteed to correspond with a significant decreases in $WS−W¯Sη$ [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

as identified by Bulat et al.⁵¹ Using an N₂ molecule with the cc-pVDZ basis and a HF target density, the reciprocal derivative of the L-curve plotted for each η in Fig. 9(b), highlights the usefulness of this approach in selecting an optimal value; η ≈ 10^−4.2 in the given example.

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.

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.

The kspies.util.wfnreader function loads wfn file formats as generated by many quantum chemistry packages, such as Gaussian,⁴⁰ ORCA,⁴¹ GAMESS,⁴² and Molpro,⁴³ 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.

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⁻⁷.

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.