Numerical Methods for the Inverse Problem of Density Functional Theory

TL;DR

This paper reviews numerical methods for solving the inverse problem in density functional theory, focusing on algorithms, error analysis, and comparing different inversion techniques on model systems.

Contribution

It introduces and compares two inversion methods based on PDE-constrained optimization and variational ideas for the inverse DFT problem.

Findings

Two inversion methods are introduced and analyzed.

Comparison of methods on finite and periodic systems.

Insights into error behavior and algorithm performance.

Abstract

The inverse problem of Kohn-Sham density functional theory (DFT) is often solved in an effort to benchmark and design approximate exchange-correlation potentials. The forward and inverse problems of DFT rely on the same equations but the numerical methods for solving each problem are substantially different. We examine both problems in this tutorial with a special emphasis on the algorithms and error analysis needed for solving the inverse problem. Two inversion methods based on partial differential equation constrained optimization and constrained variational ideas are introduced. We compare and contrast several different inversion methods applied to one-dimensional finite and periodic model systems.

Tables

Table 1. Table 1: The labels used to identify each DFT inversion method in Sec.  3 .

Inversion Method	Label
PDE-constrained optimization	PDE
van Leeuwen and Baerends	vLB
Wu and Yang	WY
Constrained Variational	CV

Table 2. Table 2: The notation used in the discrete-adjoint derivations.

Notation	Description
$𝐢$	Discrete spatial indices $i_x,i_y,i_z$
${\sum }_{𝐢=\mathrm{𝟏}}^{N_{𝐱}}$	Sum over all spatial indices ${\sum }_{i_x=1}^{N_x}{\sum }_{i_y=1}^{N_y}{\sum }_{i_z=1}^{N_z}$
$F_n$	Cost functional
$ℒ$	Lagrangian
$\tilde{n}$	Target density
$N_{\text{orbs}}$	Number of Kohn-Sham orbitals
${\varphi }^m$ and ${\chi }^m$	mth orbital and corresponding adjoint orbital
$\underset{¯}{\varphi }$	All $N_{\text{orbs}}$ orbitals grouped into a vector
$w_{𝐢}$	Optional cost-functional weight
$d{x}_{𝐢}$ and $d{k}_{𝐣}$	Integration weights

Equations

$${\epsilon }_j\underline{{\varphi }_j\left(\mathbf{r}\right)}$$

$$\underline{n\left(\mathbf{r}\right)}$$

$$v_{\text{KS}}\left(\left[\underline{n}\right],\mathbf{r}\right)$$

$$v_{\text{H}}\left(\left[n\right],\mathbf{r}\right)=\int \text{d}{\mathbf{r}}^{\prime }n\left({\mathbf{r}}^{\prime }\right)\frac{1}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}$$

$$⟨{\varphi }_j⟩{\varphi }_k={\delta }_{j,k}=\{\begin{array}{ll}1& \text{if }j=k\\ 0& \text{otherwise}\end{array}.$$

$${\epsilon }_j\underline{{\varphi }_j\left(\mathbf{r}\right)}$$

$$n\left(\mathbf{r}\right)$$

$$v_{\text{KS}}\left(\left[n\right],\mathbf{r}\right)$$

$$\frac{1}{h^2}\left[\begin{array}{ccccccc}\frac{15}{4}& -\frac{77}{6}& \frac{107}{6}& -13& \frac{61}{12}& -\frac{5}{6}& 0\\ \frac{5}{6}& -\frac{5}{4}& -\frac{1}{3}& \frac{7}{6}& -\frac{1}{2}& \frac{1}{12}& 0\\ -\frac{1}{12}& \frac{4}{3}& -\frac{5}{2}& \frac{4}{3}& -\frac{1}{12}& 0& 0\\ 0& -\frac{1}{12}& \frac{4}{3}& -\frac{5}{2}& \frac{4}{3}& -\frac{1}{12}& 0\\ 0& 0& -\frac{1}{12}& \frac{4}{3}& -\frac{5}{2}& \frac{4}{3}& -\frac{1}{12}\\ 0& \frac{1}{12}& -\frac{1}{2}& \frac{7}{6}& -\frac{1}{3}& -\frac{5}{4}& \frac{5}{6}\\ 0& -\frac{5}{6}& \frac{61}{12}& -13& \frac{107}{6}& -\frac{77}{6}& \frac{15}{4}\end{array}\right],$$

$$\Delta v\left(\mathbf{r}\right)=-4\pi n\left(\mathbf{r}\right)$$

$$v_{\text{KS}}\left(\left[n\right],\mathbf{r}\right)$$

$$v\left(x\right)=$$

$$+\left[{\delta }^{\prime }\left(x-4\right)-{\delta }^{\prime }\left(x+4\right)\right]/2,$$

$$n_0^{\text{exact}}\left(x\right)=\frac{2}{L}{\sin }^2\left(\frac{\pi x}{L}\right)$$

$$v\left(x\right)=\{\begin{array}{ll}0& \text{if }0<x<L\\ \infty & \text{otherwise}\end{array}.$$

$$n_0^{\text{noise}}=n_0^{\text{exact}}\left[1+N\left(0,1\times 10^{-10}\right)\right],$$

$$v\left(\mathbf{r}\right)=\frac{2}{n\left(\mathbf{r}\right)}\sum _{j=1}^N{\varphi }_j^{\ast }\left(\mathbf{r}\right)\frac{\Delta }{2}{\varphi }_j\left(\mathbf{r}\right)+{\epsilon }_j{\left|{\varphi }_j\left(\mathbf{r}\right)\right|}^2.$$

$$v^{k+1}\left(\mathbf{r}\right)$$

$$=\frac{1}{\tilde{n}\left(\mathbf{r}\right)}\left[n^k\left(\mathbf{r}\right)v^k\left(\mathbf{r}\right)\right]=\frac{n^k\left(\mathbf{r}\right)}{\tilde{n}\left(\mathbf{r}\right)}{v}^k\left(\mathbf{r}\right),$$

$$\underset{\mathbf{r}}{\max }\left|1-\frac{n^k\left(\mathbf{r}\right)}{\tilde{n}\left(\mathbf{r}\right)}\right|<ϵ$$

$$1-\delta <\gamma \frac{n^0\left(\mathbf{r}\right)}{\tilde{n}\left(\mathbf{r}\right)}<1+\delta ,$$

$$v^{k+1}\left(\mathbf{r}\right)=v^k\left(\mathbf{r}\right)+\gamma \frac{n^k\left(\mathbf{r}\right)-\tilde{n}\left(\mathbf{r}\right)}{\tilde{n}\left(\mathbf{r}\right)}$$

$$W_s\left[v\left(\mathbf{r}\right)\right]=2\sum _{j=1}^N⟨{\varphi }_j⟩\hat{T}{\varphi }_j+\int \text{d}\mathbf{r}v\left(\mathbf{r}\right)\left[n\left(\mathbf{r}\right)-\tilde{n}\left(\mathbf{r}\right)\right]$$

$$\frac{\delta W_s\left[v\left(\mathbf{r}\right)\right]}{\delta v\left(\mathbf{r}\right)}=n\left(\mathbf{r}\right)-\tilde{n}\left(\mathbf{r}\right).$$

$$F\left[v\right]=\frac{1}{2}{\Vert \sqrt{w\left(\mathbf{r}\right)}\left[n\left(\left[v\right],\mathbf{r}\right)-\tilde{n}\left(\mathbf{r}\right)\right]\Vert }_2^2$$

$${\epsilon }_j\sqrt{\tilde{n}\left(\mathbf{r}\right)}{g}_j\left(\mathbf{r}\right)$$

$$n\left(\mathbf{r}\right)$$

$${\epsilon }_j{g}_j\left(\mathbf{r}\right)=-\frac{1}{2}\left[\frac{\Delta \sqrt{\tilde{n}\left(\mathbf{r}\right)}}{\sqrt{\tilde{n}\left(\mathbf{r}\right)}}+2\frac{\nabla \sqrt{\tilde{n}\left(\mathbf{r}\right)}}{\sqrt{\tilde{n}\left(\mathbf{r}\right)}}\cdot \nabla +\Delta +v\left(\mathbf{r}\right)\right]g_j\left(\mathbf{r}\right).$$

$$\sqrt{\tilde{n}\left(\mathbf{r}\right)}=\exp \left\{\frac{1}{2}\log \tilde{n}\left(\mathbf{r}\right)\right\}$$

$$\frac{\nabla \sqrt{\tilde{n}\left(\mathbf{r}\right)}}{\sqrt{\tilde{n}\left(\mathbf{r}\right)}}$$

$$\frac{\Delta \sqrt{\tilde{n}\left(\mathbf{r}\right)}}{\sqrt{\tilde{n}\left(\mathbf{r}\right)}}$$

$$\begin{gathered} J\left(\left\{{\varphi }_j\right\}\right)=\int \text{d}\mathbf{r}\{\sum _{j=1}^{N_{\text{orbs}}}{\left|\nabla {\varphi }_j\right|}^2+v_{\text{KS}}\left(\mathbf{r}\right)\left[\sum _{j=1}^{N_{\text{orbs}}}{\varphi }_j^2\left(\mathbf{r}\right)-\tilde{n}\left(\mathbf{r}\right)\right] \\ +\sum _{j=1}^{N_{\text{orbs}}}\sum _{k=j}^{N_{\text{orbs}}}{\epsilon }_{j,k}\left[{\varphi }_j\left(\mathbf{r}\right){\varphi }_k\left(\mathbf{r}\right)-{\delta }_{j,k}\right]\}, \end{gathered}$$

$$R\left(\left\{g_j\right\},\alpha \right)=\alpha \sum _{j=1}^{N_{\text{orbs}}}\int \text{d}\mathbf{r}{\left|\nabla g_j\left(\mathbf{r}\right)\right|}^2,$$

$$\int \text{d}\mathbf{r}\left[\sum _{j=1}^{N_{\text{orbs}}}{g}_j^2\left(\mathbf{r}\right)\right]-1.$$

Full text

1 Introduction

Density functional theory (DFT) is a very popular and formally exact quantum many-body theory that has successfully been used to model both finite and periodic systems.1, 2 Much of this success stems from the availability of approximate exchange-correlation functionals used as input in the direct problem of DFT to produce the ground-state electronic density of a given system. Although far less common, this process can be reversed and the true exchange-correlation potential can be recovered from a given electronic density by solving the inverse problem of Kohn-Sham (KS) DFT.111Other flavors of DFT exist but we focus on the standard ground-state version of KS DFT and drop the KS prefix for brevity throughout this tutorial. Although both the forward and inverse problems use the same equations, the difference in input and output variables leads to major differences in solution methods for the two problems. We refer to the inverse problem of DFT as a density-to-potential inversion and the numerical methods required for solving this problem are the main focus of this tutorial. In this section we give a brief introduction to inverse problems and explain how density-to-potential inversions can be used within DFT.

Inverse problems are common in science and have been central in quantum mechanics since its inception. Much of what we know about the structure of matter has come from scattering experiments,3 which can be described mathematically as inverse problems. A variety of mathematical and numerical techniques exist for solving scattering problems as well as other inverse problems.4 These methods are often very different from the methods used for solving the direct problems due to the differences in mathematical structure and input data.

According to Hadamard,5 a problem is well-posed if a solution exists, it is unique, and it depends continuously on the data. If any of the three properties listed above is violated then the problem is ill-posed. The forward or direct problem of DFT is widely considered to be well-posed thanks to several proofs including the Hohenberg-Kohn Theorems and the constrained-search formulation of Levy and Lieb.1 The inverse problem of DFT can also be well-posed for discretized systems as proven in Ref. 6 but often errors and missing information in the input density lead to violations of Hadamard’s well-posed conditions. In such cases, special precautions must be taken in designing inversion methods that converge to the true solution of the inverse problem.

We have found that understanding the numerical errors involved in the inverse problem of DFT is extremely helpful in developing successful DFT inversion methods. Appendix A and Sec. 2.2.1 contain a number of examples showing how the numerical error in both the inversion algorithm and the target density being inverted set bounds on the inverted potential’s accuracy. A central theme in this tutorial is the creation of robust inversion methods that use information about the numerical errors involved in order to produce accurate potentials. This information can be used to properly smooth target densities or choose a form of regularization. Regularization is the general term used to describe a method for solving ill-posed problems. Some regularization methods are as simple as stopping an optimization according to an a posteriori parameter choice rule while others, such as total variation, can be sophisticated and very challenging mathematical and programming exercises.7 Often the form of regularization is chosen to penalize features known to be incorrect based on a priori information about the solution. For a thorough review of regularization and inverse problem theory see Refs. 7, 8, 9. The density-to-potential inversion methods presented in this tutorial are designed explicitly to include regularization.

We make a special effort to verify the integrity and robustness of the inversion algorithms presented here. Our test cases are all noninteracting systems because we can focus on recovering only the known external potential as explained in Sec. 3. These systems may seem trivial but their corresponding inverse problems are just as difficult as the inverse problems corresponding to interacting densities and provide great insight into the design of inversion algorithms. We also strive to use different algorithms and grids for the direct problems in the examples with approximate target densities to avoid any reliance on inverse crimes as described in App. A.1. The numerical errors in both the inversion methods and the target densities used as input are also tracked carefully to avoid problems with overfitting. In general, a great deal of numerical analysis is used throughout this work to make sure that our inversion methods are both correct and efficient.

One possible motivation for studying the inverse problem of DFT is to benchmark the many approximate exchange-correlation functionals currently available. Benchmarking is usually done on model systems where the density is either known analytically or found numerically through a direct solution of the electronic Schrödinger equation. As seen in many works on inverse problems in DFT,10, 11, 12, 13 major features of the exact exchange-correlation potential are often found to be missing in popular approximations and thus the inverse problem can help guide the development of new approximations.

Perhaps an even more practical reason to study the inverse problem of DFT comes from its connection with the optimized effective potential (OEP) method. In the OEP method, the exchange-correlation functional is an explicit functional of the KS orbitals and therefore an implicit functional of the density.14 As pointed out in Ref. 13, the optimization procedures used in the OEP method can be adapted to perform density-to-potential inversions and vice versa. This connection with the OEP method shows that density-to-potential inversions are also important predictive tools in addition to their use in benchmarking or guiding the construction of potentials.

2 Direct and Inverse Methods

We begin our tutorial by highlighting the differences between the forward and inverse problems of DFT. Although the direct and inverse problems of DFT share the same KS equations, there are several important differences between them that lead to very different algorithms and convergence criteria. We briefly describe common numerical methods for solving the direct problem of DFT and then describe in detail several inverse-problem methods. All of the examples presented are one-dimensional systems for simplicity but the formulas are written for multidimensional systems to aid future development. We also limit our studies to spin-compensated systems and drop the spin index on most formulas.

The KS equations for a closed-shell system with $2N$ electrons can be written as

[TABLE]

where the underlined quantities are unknowns, $v_{\text{xc}}=\frac{\delta E_{\text{xc}}\left[n\right]}{\delta n}$ is the exchange-correlation potential,

[TABLE]

is the Hartree potential, the external potential $v_{\text{ext}}\!\left({\bf r}\right)$ is usually the electron-nuclei potential, and the orbitals are constrained to be orthonormal:

[TABLE]

The orbitals are usually chosen such that their corresponding eigenvalues $\left\{\varepsilon_{j}\right\}$ form the set of lowest possible eigenvalues. The KS equations for the inverse problem of DFT

[TABLE]

only differ from the direct problem in the unknowns as indicated by underlining.

Equation (1a) is a nonlinear eigenvalue problem due to the KS potential’s density dependence whereas Eq. (4a) is a linear eigenvalue problem. The inverse problem is still nonlinear but the nonlinearity comes from the definition of the density shown in Eq. (4b). This subtle difference in nonlinearity means that the mixing schemes prevalent in numerical methods for the direct problem are not needed in solving the eigenvalue problem in the inverse problem. Additionally, the definition of the density in Eq. (4b) puts constraints on the values of the KS orbitals when solving the inverse problem that aren’t present in the direct problem. These additional constraints are exploited in the inversion method described in Sec. 2.2.4. Finally, the nonlinearity in both problems makes the choice of initial guess very important as convergence is not guaranteed in many numerical methods when poor initial guesses are used.15

The KS equations apply to both finite and extended systems provided that the appropriate boundary conditions are applied. In the direct problem of DFT the orbitals are usually constrained to satisfy periodic or box-type (Dirichlet) boundary conditions as dictated by the system being studied.16 In the inverse problem the boundary conditions may not be known exactly but the density constraint uniquely determines the orbitals’ behavior at the boundaries. If boundary conditions are imposed then it is important that they agree with the density constraint as explained in App. A.1. We show in this tutorial the alternative of letting the inversion method enforce the orbital boundary conditions to agree with the density constraint. This method has the advantage that no detailed knowledge of the boundary conditions is needed and the disadvantage that the values of the potential at the boundaries may need to be discarded after the inversion if the density at the boundaries contains some inherent error. We call this method density-constrained boundary conditions (DCBC) and discuss it further in Sec. 2.2.1.

As mentioned in Sec. 1, uniqueness is a key ingredient in solving both direct and inverse problems. Although the density in DFT is a unique quantity in both the forward and inverse problems, the KS orbitals that it is built from are not unique. As seen in Eq. (1b), the density is formed from the absolute value squared of the orbitals and therefore the orbitals can only be unique up to a phase factor. Furthermore, any unitary transformation of the KS orbitals will also produce the same density. The KS potential in Eq. (1c) is also a source of nonuniqueness as it is only unique up to a constant.6 As seen in Sec. 2.2, most density-to-potential inversion methods force the unknown potential to be unique by imposing restrictions of some form.

Numerical convergence in the direct problem of DFT involves both a convergence in the iterations to self-consistency and a convergence for the eigenvalue problem solved at each iteration. This convergence is described in great detail in Ref. 15 and we only point out here that the potentials are usually given by analytic formulas. This means that the KS equations can be solved self-consistently and the error in the resulting density will depend mainly on the level of discretization employed and the choice of approximate exchange-correlation functional. Numerical convergence in the inverse problem of DFT is also limited by the chosen discretization method but it is the quality of the target density that sets a somewhat surprising limit on convergence as explained in this tutorial.

Although most of the test cases presented in this work have analytic formulas for the target densities, most target densities of real interest are not known analytically. If the numerical error present in a given target density is not accounted for properly in a density-to-potential inversion then overfitting becomes a problem and unphysical features develop in the recovered potential. Examples of overfitting and convergence limits related to the relative error of the target density are given in App. A.

2.1 Direct-Problem Methods

A variety of numerical methods for solving the KS equations exist and can be placed in essentially three different categories: plane waves and grid methods, localized atomic-(like) orbitals, and atomic sphere methods. 1 We use the finite difference grid method for our calculations similar to the Octopus time-dependent density functional theory (TDDFT) code.17 The finite difference method is capable of producing sharp features that would be difficult to find using other common electronic-structure basis sets. In this section we review several aspects of the finite difference formulation of DFT that are necessary for solving the direct problem and also relevant to some of the inverse problem methods. A more thorough review of the finite difference method applied to DFT can be found in Ref. 18.

The orbitals, potentials, and density in the KS equations are represented as discrete points on a grid when using the finite difference method. The potentials are diagonal in the coordinate representation and the Laplacian is a very sparse matrix operator depending on the approximation used. An example fourth-order Laplacian for a regularly spaced grid with seven points is

[TABLE]

where $h$ is the grid spacing and sided finite differences are used at the boundaries. Different orders of approximate finite difference operators can be derived using the algorithm described in Ref. 19. There are many possible boundary conditions and methods to implement them depending on the problem and accuracy desired.20, 21, 22 Custom meshes can also be used that conform to the geometry of a given system as explained in Ref. 16. One disadvantage of the finite difference method applied to DFT is the lack of a variational principle for most approximate finite difference operators. This means that the approximate energy may actually be smaller than the true energy found by solving the KS equations exactly. The approximate energy will approach the true energy as the grid is refined and/or the accuracy of the finite difference operator is improved but it doesn’t have to approach it from above unless the finite difference operator is modified as in Ref. 23.

A variety of numerical algorithms are available for computing the potentials in Eq. (1c). The Hartree potential $v_{\text{H}}$ shown in Eq. (2) satisfies the Poisson equation

[TABLE]

and has been studied extensively. The fast Fourier transform, fast multipole method, and conjugate gradient methods are just a few of the available algorithms for computing $v_{\text{H}}$; a thorough review of these methods can be found in Ref. 24. As mentioned in Sec. 1, there are many approximate exchange-correlation functionals and corresponding potentials. The libxc library contains a wide selection of these approximations and is easily incorporated in DFT programs.25 Exchange-correlation approximations can also be derived for systems with custom electron-electron interactions as shown in Refs. 26 and 12.

The KS equations are usually solved using sparse matrix eigenvalue solvers in contrast to the dense solvers typically used in basis-set methods.18 The examples in this work use the eigs and eig_banded eigenvalue solvers in the SciPy library, which are interfaces to routines in the lapack and arpack libraries respectively.27 These solvers are capable of computing the eigenvalues and eigenvectors of a sparse matrix to machine precision and are very efficient when only a few eigenvalues and eigenvectors are needed. If the potentials used in a given calculation are exact and the solvers are allowed to converge to machine precision then the majority of the numerical error will come from the discretization of the Laplacian. An example of this error is shown in Fig. 1 for the ground state density of the harmonic oscillator using a second-order approximate Laplacian. Although the majority of the error is in the central region, Fig. 2 shows that the relative error is actually largest in the asymptotic regions. The approximate finite difference operators are very accurate in regions where the orbitals are well represented by polynomials and are less accurate in the asymptotic region where the orbitals are dominated by exponential decay. The large relative error is largely inconsequential for most direct DFT problems with almost no effect on the total energies but it does play a strong role in the inverse problem as shown in Sec. 2.2. This example also shows that the errors in most direct DFT problems are systematic rather than random noise typical of experimental data and must also be accounted for in the inverse problem as explained in App. A.

The previous example of a noninteracting system only required solving one eigenvalue problem but a typical DFT calculation for an interacting system requires multiple eigenvalue solves in what is known as the self-consistent field (SCF) method. The initial density guess for a given system is inserted into the right-hand side of Eq. (1a), the eigenvalue problem is solved for a new set of orbitals, the orbitals are inserted into Eq. (1b) to produce a new density, and the cycle repeats until a convergence criterion is reached. This method is usually modified to avoid oscillating densities by mixing the new density with previous densities via the direct inversion in the iterative subspace (DIIS) method28) or some other form of mixing.16 Often the SCF loop is stopped when the change in density or the change in each eigenvalue becomes smaller than a prescribed value but this does not guarantee global convergence. An alternative gradient-based method that avoids the repeated eigenvalue computation is also possible as illustrated in Ref. 29. A more detailed analysis of numerical convergence in DFT calculations can be found in Ref. 15.

2.2 Inverse-Problem Methods

As is the case with the direct problem, there exist an assortment of numerical methods for solving the inverse problem of DFT, each with its own set of advantages and disadvantages. Here we review several density-to-potential algorithms and introduce two new inversion methods. All of the methods are presented using the finite-difference method and may include small adaptations made to convert algorithms that were originally designed for basis-set methods.

2.2.1 One and Two Electrons: One-Orbital Formula

The one-orbital inversion formula

[TABLE]

can be derived simply by setting the energy to zero and solving for the KS potential in Eq. (1) in terms of the ground state density $n_{0}$. (We are allowed to set the energy to zero because the potential is only unique up to a constant.) This formula is exact for one electron or two electrons with opposite spins in the same orbital and is often referred to as the bosonic or one-electron potential.30 Although at first glance the one-orbital formula appears somewhat trivial, it is helpful in unraveling many of the numerical problems common to all density-to-potential algorithms and can even serve as a useful approximation for systems with many electrons.

The one-orbital formula makes it very clear that the KS orbital is uniquely determined up to an arbitrary phase by the target density via the equation $\phi_{0}=\sqrt{n_{0}}$. This means that the values of the orbital in the entire domain of the problem are known including the boundaries. In the forward problem of DFT we usually force the orbitals to satisfy certain boundary conditions but in the inverse problem this is not necessary because the target density uniquely determines their behavior at the boundaries up to an arbitrary unitary transformation of the orbitals. Both Fig. 14 and the $g_{j}^{\prime}\!\left(\pm L\right)\approx 0$ line in Fig. 9 show examples of molecular system inversions with approximate boundary conditions enforced. In both examples the resulting potentials have several incorrect values near the boundaries where the potential actually enforces the approximate boundary conditions. In Fig. 14 the potential should actually be

[TABLE]

because the density is artificially forced to be zero at the boundaries and the “incorrect” potential values near the boundaries are approximations to this discontinuous behavior, where $\delta\left(x\right)$ is the Dirac delta function and $\Theta\left(x\right)$ is the Heaviside function. When approximate boundary conditions are not enforced in these examples then the potential is correct everywhere including at and near the boundaries.

Although the one-orbital formula is exact for one-electron and most closed-shell two-electron systems, it must be carefully implemented numerically to avoid spurious results. As seen in App. A.2, noisy densities can easily produce wildly inaccurate potentials when using Eq. (7) and finite difference operators. Such noisy behavior is unphysical because densities arising from the Schrödinger equation are theoretically continuous even when generated by singular potentials.3 The following examples show how to appropriately smooth target densities based on the noise level before inserting them into the one-orbital formula.

We first illustrate smoothing using the particle-in-a-box3 ground state density

[TABLE]

corresponding to the potential

[TABLE]

If we take the exact ground-state density $n_{0}^{\text{exact}}$ and add a small amount of weighted noise generated from the standard normal distribution $\operatorname{N}\!\left(\mu,\sigma^{2}\right)$31 according to the formula

[TABLE]

then the potential resulting from the one-orbital formula will also be noisy with a level proportional to the number of grid points. This noise can be removed by fitting a cubic spline to the target density weighted by the inverse of the approximate standard deviation $w=1/(\sigma*n_{0}^{\text{approx}})$. We implement this smoothing procedure using the UnivariateSpline routine in the SciPy library and the results are shown in Fig. 3 for a box of length 5 with varying numbers of grid points. This example shows that smoothing is essential for dense grids with errors but not as important on coarse grids since the derivative of the error doesn’t dominate in the one-orbital formula as explained in App. A.3.

The target density in the previous example was well approximated by a cubic spline but other target densities may need additional modifications before smoothing is applied. In particular, exponentially decaying densities are not adequately described by cubic splines unless they are first logarithmically transformed. We show how this can be done using the ground-state of the harmonic oscillator with the same weighted noise given by Eq. (11). In this case we take the logarithm of the noisy density, fit this new quantity to a cubic spline using the weights $w=1/\log(1+\sigma)$, and then exponentiate the result to get the smoothed target density. The potential produced using this logarithmic scaling before fitting is much more accurate in the asympototic region than the potential produced with no scaling before smoothing as seen in Fig. 4. We refer to this smoothing procedure later in this tutorial as the logarithmic smoothing method. Even highly singular potentials can be recovered using this method provided that the locations of zero density are first averaged with their nearest neighbors before applying the scaling and smoothing. (This or a similar procedure must be followed to avoid division-by-zero errors.) Figure 4 shows the case described in Ref. 32 of treating the density of the first excited state of the simple harmonic oscillator as a ground state density. The singular feature of the potential ($\delta\!\left(x\right)/\left|x\right|$) becomes more pronounced with this procedure as the grid is refined.

In the previous examples we were able to correctly smooth the target densities by having estimates of the standard deviation at each grid point. Such information would be typical of an experimentally determined density but most of the error patterns encountered in our density-to-potential inversions come from computationally determined densities using interacting wave function methods. In these cases we can use error estimates to compute the correct weights for the smoothing procedure. We show how this can be done by computing the ground state density of the harmonic oscillator using the infinite-to-finite spatial grid mapping of Ref. 21 with a scaling parameter of $\alpha=2$. (We purposely use a different grid to compute the target density in order to avoid comitting an inverse crime.) We compute the target density on this grid once with a second-order discretization of the Laplacian and 101 points, (labeled $n^{\mathcal{O}\!\left(2\right)}$), and then again with a fourth-order approximate Laplacian and 201 points, (labeled $n^{\mathcal{O}\!\left(4\right)}$). These two densities are linearly interpolated to the same equally-spaced grid shown in Fig. 4 and subtracted from each other to form an error estimate at each grid point. We then apply the logarithmically scaled spline fitting procedure described above with the approximate weights $w=1/\left|\log n^{\mathcal{O}\!\left(2\right)}-\log n^{\mathcal{O}\!\left(4\right)}\right|$ and the exact weights $w=1/\left|\log n^{\mathcal{O}\!\left(2\right)}-\log n^{\mathrm{exact}}\right|$. The results of this calculation are displayed in Fig. 5 and show that the error pattern of the approximate target density is far more systematic than the random noise error patterns in the previous examples. The smoothing using approximate weights underestimates the potential at the edges of the box because the error estimate is too small in that region. If the error estimate is computed using 501 points instead of 201 points in the fourth-order density computation then the potential is recovered correctly in the entire box. This example illustrates the need not only for error estimates but also for the need to put bounds on the error estimates either through computing additional error estimates or through statistical means.

Even if the target density is correctly smoothed before applying the one-orbital formula, the finite difference operator used to approximate the Laplacian can severely limit the accuracy of a density-to-potential inversion. As was the case with the spline fitting procedure described above for exponentially decaying densities, the finite difference method performs poorly in the asymptotic region unless extra precautions are taken. In Sec. 2.2.3 we show how the finite difference operators can be modified to treat much of the exponential decay analytically through a scaling procedure. Figure 6 shows a comparison of the one-orbital formula applied to the exact ground state of the harmonic oscillator with and without scaling the finite difference operator.

Although the one-orbital formula is only exact for one- or two-electron systems, it is a remarkably good approximation for systems with large regions where the target density is dominated by one orbital. We illustrate the power of this approximation by applying it to a target density formed from the sum of the ground state and first two excited states of the harmonic oscillator. The harmonic potential is correctly recovered in the asymptotic region where most of the density comes from the second excited state as seen in Fig. 7. The potential in the central region is not correct but it is a smooth approximation and is the same order of magnitude as the correct potential. Furthermore, the cost of this approximation is simply one sparse matrix-vector multiplication followed by an array division, which is negligible in comparison to the many sparse matrix-vector multiplications typically needed to solve the eigenvalue problem in Eq. (4a).

The examples in this section demonstrate how the convergence criteria for DFT inverse problems are different from those of the direct problem. These examples and App. A clearly show that the size and location of the relative error in the target density is the main constraint on the accuracy of the potential resulting from a density-to-potential inversion. In many cases the target density must first be smoothed before performing an inversion and its standard deviation is the key ingredient in performing this smoothing operation correctly. We chose to use a spline-based smoothing filter because it is easy to incorporate the standard deviation but other methods can be used with similar success such as the Savitzky-Golay filter.33 The accuracy of the finite difference operators also puts a strict limit on the quality of the recovered potential as some operators simply aren’t good approximations in certain regions of rapidly changing density. Furthermore, the computational cost of the inverse problem using the one-orbital formula is much lower in these examples than in the forward problem and can serve as a useful approximation for many-electron systems.

2.2.2 Standard Methods

A number of density-to-potential inversion methods have been developed besides the one-orbital inversion formula. In this section we look at the iterative inversion method of Ref. 11 and the direct optimization method of Ref. 13. A more complete review of the available density-to-potential inversion methods is found in Ref. 13. The two inversion schemes described in this section are fairly representative of the available methods and serve as useful benchmarks later in this chapter. The target density in the remainder of this work is written as $\tilde{n}$ in order to distiguish it from the intermediate densities created in the iterative procedures that we describe here. When it does not lead to confusion, we also drop the subscript KS on the KS potential $v_{\text{KS}}\!\left(\mathbf{r}\right)$ and the subscript orbs on the number of orbitals $N_{\text{orbs}}$ to simplify the presentation.

Van Leeuwen and Baerends: The iterative inversion method of van Leeuwen and Baerends described in Ref. 11 is derived by multiplying both sides of Eq. (4a) by $\phi_{j}^{*}$, summing over $j$, and then dividing by the density given in Eq. (4b) to produce the formula

[TABLE]

This formula is then turned into an iterative scheme by placing the target density in the denominator and iterating until convergence via the formula

[TABLE]

where $k$ is the iteration index. This iterative scheme simply increases $v^{k+1}\!\left(\mathbf{r}\right)$ in regions where $n^{k}\!\left({\bf r}\right)>\tilde{n}\!\left({\bf r}\right)$ and decreases it in regions where $n^{k}\!\left({\bf r}\right)<\tilde{n}\!\left({\bf r}\right)$ as seen in the last line of Eq. (13). The iterations are stopped when

[TABLE]

for some desired threshold $\epsilon$. This scheme is sensitive to the initial potential guess $v^{0}\!\left(\mathbf{r}\right)$ and usually requires a prefactor to avoid wild potential oscillations. Reference 11 suggests choosing this prefactor by enforcing the condition

[TABLE]

where $\delta\approx 0.05$ and $\gamma$ is the prefactor.

This iterative inversion scheme of van Leeuwen and Baerends is very easy to implement as it only requires a partial differential equation (PDE) solver to produce the new density $n^{k}\!\left({\bf r}\right)$ at each iteration. We do, however, make some modifications to the original scheme by not forcing the potential to approach any given limit at the boundaries and by using the formula

[TABLE]

instead of Eq. (13). We don’t enforce boundary conditions on the potential because many of the test cases we use recover the external potential with unknown boundary conditions as opposed to just the exchange-correlation potential with its known asymptotic behavior.11 (Note that the orbitals may still have boundary conditions and the above discussion applies only to fixing values of the potential at the boundaries.) We use Eq. (16) instead of Eq. (13) because it still uses the same principle of increasing (decreasing) the potential in regions where the density is too large (small) but does not need boundary conditions to function properly. Other modifications are also possible such as the use of an approximate density-density response matrix to help guide the iterations.34

This method produces smooth potentials by construction provided that the underlying eigenvalue solver produces smooth orbitals. The reliance on a prefactor, however, requires some manual intervention or heuristics to ensure convergence. Large prefactors may result in wildly oscillating potentials that never converge while small prefactors may produce very small potential changes that take an inordinate amount of time to converge. Choosing a prefactor is usually not a problem for isolated inversions but can become very tedious when repeatedly applying the procedure to cases with many similar inversions such as in the creation of dissociation curves.35

Wu-Yang Algorithm The Wu-Yang inversion algorithm first introduced in Ref. 13 is one of the more sophisticated and versatile DFT inversion methods. In this inversion algorithm the Levy constrained-search formulation of DFT is modified to produce an unconstrained optimization method. The functional

[TABLE]

is maximized to produce the optimal potential corresponding to the target density $\tilde{n}\!\left({\bf r}\right)$. (The kinetic energy is doubled because we are still assuming the orbitals are doubly occupied.) This maximization is usually performed using classical optimization routines and is very efficient because the required first derivatives (Jacobian values) are given by the simple formula

[TABLE]

There is very little extra computational effort required to compute these derivatives because the density is already computed at each iteration while computing Eq. (17). The second derivatives (Hessian values) can also be computed either through first-order perturbation theory13 or via the discrete adjoint-method described in App. B.

The Wu-Yang inversion algorithm was developed originally for the OEP method of DFT and later modified for the DFT inverse problem.36, 13 The method is efficient and robust when implemented properly but can suffer from rounding errors as shown at the end of App. A. It can also produce unphysical highly oscillatory potentials when an unbalanced basis set is used without regularization.37

2.2.3 PDE-Constrained Optimization

The inversion schemes described in the previous section rely on modifications to the KS equations or underlying energy functionals. Such modifications are not essential in creating DFT inversion methods as we show in this section through the use of PDE-constrained optimization. Although conceptually simpler than many other inversion procedures, PDE-constrained optimization can be difficult to implement due to the large number of unknowns and programming challenges involved in the optimization procedure.38 In this section we give a brief introduction to PDE-constrained optimization in the context of DFT inverse problems and address many of these challenges. We also introduce an important scaling concept that plays a critical role in increasing the accuracy of the potential in the asymptotic regions of DFT inverse problems.

PDE-constrained optimization applied to the DFT inverse problem amounts to a fitting procedure in which the unknown potential is optimized until the corresponding density matches the target density. We use the weighted least-squares cost functional

[TABLE]

as our measure of fitness, where the subscript 2 indicates the $L^{2}$ norm and $w\!\left({\bf r}\right)$ is a positive definite weighting function.39 The density $n\!\left(\left[v\right],{\bf r}\right)$ comes from solving Eqs. (4a) and (4b) for a given potential chosen during the optimization procedure. The numerical minimization of the discretized Eq. (19) usually involves a very large number of unknown potential values and generally requires the use of gradient/Jacobian-based optimization algorithms.40 The derivation, programming, and computation of these gradients, [functional derivatives of the cost functional given by Eq. (19)], is the main difficulty in applying PDE-constrained optimization to the DFT inverse problem.

We use the discrete adjoint-method to compute the cost functional derivatives of Eq. (19). This method is also called the discretize-then-differentiate method because the cost functional is first discretized and then differentiated with respect to the potential.41 Appendix B contains a complete derivation of the discrete adjoint equations for the DFT inverse problem in multiple dimensions. These derivatives are then employed in the truncated-Newton (TN) algorithm of Ref. 42 as implemented in the SciPy library27 to optimize the unknown potential until a desired density fit is found. The cost of computing the functional derivatives via the discrete adjoint-method is roughly the same as computing $F\!\left[v\right]$ itself depending on the choice of linear solver.

The PDE-constrained inversion method outlined above can also suffer from rounding errors similar to the Wu-Yang algorithm described in the previous section. These errors can be avoided simply by weighting the densities in the asymptotic region more heavily than in the high-density regions. One possibility is to set $w\!\left({\bf r}\right)=1/\left[\tilde{n}\!\left({\bf r}\right)\right]^{2}$ so that the relative error is minimized instead of the absolute error. This weighting scheme in combination with the rescaling described below allows us to accurately reproduce the unknown potential in all regions of a given density-to-potential problem.

As mentioned in Sec. 2.2.1, the finite difference operators used to approximate the Laplacian are not very accurate in regions where the exponential decay dominates the behavior of the KS orbitals. The finite difference operators would be accurate in all regions if we could apply them to the logarithm of the orbitals similar to the logarithmic spline fitting procedure described earlier. The KS orbitals, however, are not positive semidefinite like the particle density so we instead write them as $\phi_{j}\!\left(\mathbf{r}\right)=\sqrt{\tilde{n}\!\left({\bf r}\right)}g_{j}\!\left(\mathbf{r}\right)$ and solve for the scaled orbitals $\left\{g_{j}\!\left(\mathbf{r}\right)\right\}$ in the KS equations

[TABLE]

After distributing the Laplacian and canceling common factors, Eq. (20a) can be rewritten as

[TABLE]

The scaling factor can also be written as

[TABLE]

and differentiated to produce the terms

[TABLE]

The derivatives of $\log\tilde{n}\!\left({\bf r}\right)$ and $\left\{g_{j}\!\left(\mathbf{r}\right)\right\}$ are well approximated by the finite difference operators described in Sec. 2.1 because they do not decay exponentially. An example of the scaled orbitals for a system consisting of the first three orbitals of the harmonic oscillator can be seen in Fig. 8.

Although it is common to enforce boundary conditions on the orbitals when solving the KS equations, we simply let the numerical optimization routine choose the ideal boundary conditions via the DCBC method described in Sec. 2.2.1. In Fig. 8, the scaled orbitals appear to have derivatives of nearly zero near the boundaries so it may seem reasonable to apply zero derivative boundary conditions when using scaled orbitals. This assumption leads to errors near the boundaries as shown in Fig. 9 in which the PDE-constrained optimization procedure outlined above is applied to a target density formed from the sum of the ground state and first two excited states of the harmonic oscillator. This example clearly shows that enforcing approximate zero derivative boundary conditions on the scaled orbitals results in the incorrect potential at the boundary whereas the potential is recovered correctly when no boundary conditions are applied. If the target density happens to be produced using artificial boundary conditions then the potential recovered when not applying boundary conditions on the orbitals may have several incorrect values at the edge of the simulation box.

2.2.4 Constrained Variational Method

The previous inversion methods described in this chapter all require the solution of an eigenvalue problem at each iteration. In the case of small systems this is not a severe constraint but as the system size increases the solution of the eigenvalue problem begins to dominate the computation.15 Although dismissed in Ref. 13 as “difficult to perform”, the Levy constrained-search formulation can be used directly to perform a density-to-potential inversion without repeatedly solving eigenvalue problems at each iteration. In this section we show how to overcome many of the difficulties involved in this constrained-search inversion method through the use of scaling and regularization.

In the Levy constrained-search formulation of DFT, the noninteracting kinetic energy is minimized subject to the constraints that the orbitals are orthonormal and produce the molecular density of the real system. The functional to be minimized in this search can be written using Lagrange multipliers as

[TABLE]

where the orbitals are assumed to be real for simplicity of presentation. The first term in Eq. (25) is one possible form of the kinetic energy and lacks a factor of $1/2$ because the orbitals are doubly occupied. This expression can be used directly in a constrained numerical optimization program to find the orbitals and KS potential for a given target density without the need to solve an eigenvalue problem at each iteration. The orbitals themselves are being varied by the optimizer to minimize the kinetic energy while satisfying the density and orthonormality constraints. (This inverse problem method is similar to the gradient-based forward method described in Ref. 29.) It does, however, suffer from the same rounding errors involved in the Wu-Yang algorithm as it still involves a minimization of the noninteracting kinetic energy. Furthermore, the orthonormality and density constraints can easily lead to unphysical orbitals as they compete with one another in the optimization.

In order to resolve the issues mentioned above, we use the scaled orbitals introduced in Sec. 2.2.3 and regularize them by adding the penalty functional

[TABLE]

to the cost functional $J\left(\left\{\phi_{j}\right\}\right)$, where $\alpha>0$ is the regularization parameter chosen via the discrepancy principle.7 The density constraint with the scaled orbitals is

[TABLE]

No boundary conditions are applied to the orbitals for reasons explained in Sec. 2.2.3.

Optimizing the noninteracting kinetic energy using the scaled and regularized orbitals produces a set of orthonormal orbitals that are a unitary transformation of the KS orbitals.43 We use the ipopt library44 to perform this constrained minimization. The KS potential can be extracted from these orbitals by solving the set of linear equations

[TABLE]

where $\alpha,\beta=1,\ldots,N_{\text{orbs}}$. This linear system can be solved efficiently using a sparse least-squares solver in which the full matrix is never explicitly created. An example of this density-to-potential inversion procedure with two different regularization parameters is shown in Fig. 10 for a target density formed from the sum of the ground state and first two excited state densities of the harmonic oscillator.

3 Inversion Examples

We have chosen to use well-known noninteracting quantum systems with which most readers are familiar to help illustrate the inversion methods described in Sec. 2.2. The KS potentials recovered for these systems are equal to the known external potential because the systems are noninteracting and this allows us to focus on the numerical properties of each inversion algorithm. More realistic DFT inversions use interacting densities and extract the exchange-correlation potential from the recovered KS potential using Eq. (4c). Unless otherwise noted, the examples we present all use the one-orbital approximation as the initial guess for the unknown potential. The gradients used for the PDE-constrained optimizations are given in App. B. Table 1 displays the abbreviations used in this section to distinguish the different inversion methods.

In these examples we implement scaling, different finite difference approximations, and placement of the boundary conditions in the potential for both the PDE and CV routines. We do this simply because we our more familiar with these methods and already have this capability implemented in our code. It should also be possible to modify both the vLB and WY methods in similar ways to implement these features so that all of the schemes have similar error patterns in the asymptotic regions. All of the error patterns in these examples are shifted relative to one another so as to see the differences and make better comparisons.

3.1 Harmonic Potential Inversion

Our first inversion example is a noninteracting system of six electrons in the harmonic potential $v\!\left(x\right)=\frac{1}{2}x^{2}$. The target density is constructed from the lowest 3 orbitals

[TABLE]

with each orbital doubly occupied as indicated in Eq. (1b). The potentials recovered from this density using several different inversion methods are shown in Fig. 11 using a grid of 101 equally-spaced points running from -8 to 8. The potential is only unique up to a constant so we also plot the differences between the recovered and exact potentials with the more constant differences indicating better recovered potentials.

In this example, the rounding errors in the WY algorithm’s functional dominate in the region $\left|x\right|>4$ as explained in App. A. The unscaled finite-difference scheme used in the vLB method is not accurate in the region $\left|x\right|>4$ and leads to an incorrect higher curvature in that region. Both the CV and PDE methods have similar problems when they are not scaled or use zero boundary conditions so this example is not meant to show the superiority of any one scheme but rather to show the importance of these modifications.

3.2 Morse Potential Inversion

Our second inversion example is a noninteracting system of six electrons in the Morse potential

[TABLE]

where $a=1/2$, $D_{e}=10$, and $r_{e}=3$. Analytic solutions for this system’s wave functions can be found in Ref. 45 and are used as input in Eq. (1b) for forming the target density. The potentials recovered from the target density described above are shown in Fig. 12 using a grid of 51 equally-spaced points running from 1 to 9.

This inversion example illustrates some of the numerical issues described in Sec. 2. One numerical issue is the need to adequately represent the kinetic energy operator on the chosen grid. The PDE and CV methods in this example use a fourth order approximation of the Laplacian while the other methods use a second-order approximation. This difference in accuracy on such a coarse grid leads to several wiggles in the WY and vLB potentials because the inversion routines fit to the errors in the approximate Laplacian. Additionally, the one-orbital approximation is not a sufficiently accurate initial guess for the scaled PDE-constrained optimization to locate the correct potential, hence the optimizer gets stuck in a local minimum. (The potential is so shallow that the scaled density is much too large near the right boundary.) This can be remedied by either raising the boundary values of the initial guess to force a more rapid decay or by including a small amount of regularization via Eq. (26) with $\alpha=1$ to guide the initial optimization and then turning it off with $\alpha=0$ before running the final optimization. The PDE result shown in Fig. 12 uses the regularization option but the result is nearly identical to that found when raising the boundary values of the initial guess. Both of these issues arise even when we have an analytic target density with no added error.

3.3 Kronig-Penney Inversion

Our last inversion example is a noninteracting system of six electrons per unit cell in the Kronig-Penney potential

[TABLE]

where $L=2$ is the length of the unit cell, $a=5/4$ is the width of the well, $V_{0}=-10$ is the depth of the well, and $m$ is an integer. The target density is found by extracting the density from the center cell of a finite system comprised of 15 unit cells, [$m=-7,\ldots,7$ in Eq. (33)], with 400 points per unit cell and a second-order discretization of the Laplacian. By performing a finite calculation of the target density we avoid the risk of committing an inverse crime.

The inverse problem for this example is solved on a grid with 50 points using PDE-constrained optimization without scaling and a fourth-order approximate Laplacian. (We don’t use scaling because the exponential decay doesn’t dominate anywhere in the system’s density.) The Hamiltonian is modified as explained in Ref. 16 to work with the Bloch functions $\left\{\phi_{i,{\bf k}}\!\left({\bf r}\right)\right\}$. In particular, the kinetic energy operator has a $k$-space dependence

[TABLE]

and the total density is found via an integration over the first Brillouin zone

[TABLE]

where $\Omega_{{\bf k}}$ is the volume of the first Brillouin zone. We perform this $k$-space integration using the trapezoid rule and 11 equally spaced points running from $-\pi/L$ to $\pi/L$. (For simplicity we integrate over the whole first Brillouin zone but usually symmetry is used to reduce this to an integral over the irreducible Brillouin zone.) The results of this inversion are shown in Fig. 13 and show that even sharp features such as Heaviside functions can be recovered successfully in these inversions.

Conclusions

We have presented a series of algorithms for performing density-to-potential inversions in DFT. Although the examples in this work are all one-dimensional model systems, the algorithms are written in multiple dimensions and are applicable to more realistic systems. In particular, these methods can all be used to find the exact exchange-correlation potential corresponding to a given interacting density as explained in Sec. 2. We have also presented a detailed numerical analysis that points to the relative error as the main limitation in density-to-potential inversions and the need for error estimates in performing these inversions properly.

The inversion methods presented in this tutorial use the finite difference method for solving the governing PDEs but they can also be applied to the common basis function methods described in Ref. 1. The insight we have gained through our many test cases and numerical analysis with the finite difference method can also be used to improve other inversion methods as hinted at in Sec. 3.1. In particular, the scaled-orbital regularization presented in Sec. 2.2.4 is perhaps a more natural regularization than the potential regularization presented in Ref. 37. This regularization should also make it possible to examine more realistic examples including possibly ill-posed densities from experimental data.

All of the inversion methods used in Sec. 3 have strengths and weaknesses that should be considered when choosing which method to use for a given inverse DFT problem. The PDE-constrained optimization algorithm is very much a black-box method requiring almost no physical intuition but it is only competetive with other methods when the gradient is computed efficiently and used in a good optimization routine. This is not a trivial endeavor as both the derivation and programming of these gradients can be quite involved as shown in App. B. Despite this weakness, it is very straightforward to adapt this method to different inverse problems resulting from extensions of DFT such as TDDFT.38 The method of van Leeuwen and Baerends does have a sound physical motivation and is simple to program because it only needs to solve the Kohn-Sham equations. The biggest disadvantage to this method is its reliance on a prefactor, which requires user intervention and makes automating the routine difficult. The method of Wu and Yang is similar to the van Leeuwen and Baerends method in that it is also based on some physical motivation and only needs to solve the Kohn-Sham equations provided the Hessian isn’t needed. This method is very efficient when paired with a good optimizer but can suffer from rounding and overfitting errors if extra precautions are not taken. The constrained variational method is the only method that doesn’t solve an eigenvalue problem at each iteration and this may be an advantage for large systems with more demanding eigenvalue problems. In the systems that we have studied this method appears to be insensitive to the initial guess with even random numbers for the initial orbitals resulting in rapid convergence. Its largest drawback is its current reliance on a regularization parameter but this choice can be automated using standard inverse problem methods. Since we have only considered model one-dimensional systems in this study we do not know how these methods will perform on more realistic three-dimensional systems and plan to do future studies to see how they scale with larger systems.

Our original intent in developing inversion algorithms was to benchmark approximate exchange-correlation approximations to guide new approximations within DFT. During this development, however, several new and interesting applications have come to light that warrant further study. For example, by allowing the inversion methods to enforce the boundary conditions by modifying the potential near the boundaries in the DCBC method, it should be possible to extract new approximate boundary conditions and compare them to other state-of-the-art methods.22 Likewise, the orbital scaling introduced in Sec. 2.2.3 could be modified for the direct DFT problem as an alternative to using adaptive grid methods.16 Thanks to the very similar structure between the OEP method and the DFT inverse problem, algorithm development in one area can lead directly to improvements in the other as seen in the development of the Wu-Yang inversion algorithm from a previous OEP optimization method.13 We plan to continue exploring these other areas of DFT and refine the density-to-potential inversion methods described in this tutorial.

Acknowledgments

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This work was supported by the National Science Foundation CAREER program under Grant No. CHE- 1149968. A.W. also acknowledges support from the Camille Dreyfus Teacher-Scholar Awards Program. We thank Jonathan Nafziger for many discussions related to DFT inversions. We also thank Andrew Baczewski, Attila Cangi, and the referees for their many helpful comments.

Appendix A Numerical Issues

A number of numerical problems arise in performing DFT inversions that are often very difficult to isolate. In this appendix we use very basic one-dimensional models to show how boundary conditions, relative error, and rounding errors all set convergence limits on density-to-potential inversions. We employ the one-orbital formula shown in Eq. (7) and a grid ranging from $-4$ to $4$ with 51 points in the examples of this section unless otherwise noted.

A.1 Boundary Conditions and Inverse Crimes

Many of the inversion schemes mentioned in Sec. 2.2 impose boundary conditions on the orbitals that don’t have to agree with the target density’s boundary conditions. A mismatch in boundary conditions can place severe constraints on the accuracy of density-to-potential inversions in the boundary regions as we illustrate here with an example. For simplicity, we take the exact density of the simple harmonic oscillator’s ground state and set the values at the boundaries to zero as would be typical for box-type (Dirichlet) boundary conditions. The one-orbital formula should give us the exact potential for this problem everywhere except the boundaries, (where the density was artificially set to zero), but Fig. 14 shows that the potential is also incorrect at several points next to the boundaries. The number of incorrect points near the boundaries increases as higher-order finite difference schemes are used to represent the Laplacian due to the increased nonlocality of the finite difference operators.

One common method for avoiding the mismatch in boundary conditions described above is to compute the target density with the same boundary conditions and numerical method. In Fig. 15 we compute the ground state density of the harmonic oscillator numerically with box-type boundary conditions and a second-order finite difference approximation for the Laplacian. If the same finite difference approximation is used in the one-orbital formula then the resulting potential is almost exact everywhere. This type of inversion is known as an inverse crime. It is easily seen in Fig. 15 that changing to a fourth-order approximate Laplacian with the same box-type boundary conditions does not have the same error cancellation and produces worse values than the second-order method.

The term ‘inverse crime’ is commonly used in inverse problem theory to describe an inversion method relying on a cancellation of errors between the data simulation and reconstruction methods.9 Inverse crimes are only committed in DFT inversions when the same numerical methods are used to solve both the direct and inverse problems. This is somewhat rare in performing DFT inversions because the densities being inverted are usually found using wave-function methods that employ different numerical procedures than those used in DFT inversions.26 Inversion methods that rely on inverse crimes are very limiting in the number of systems that can be studied and will likely fail when applied to experimental data.

This example of an inverse crime is not meant to discourage practitioners from using knowledge about the target density in performing inverse DFT problems. Much of Sec. 2.2.1 is devoted to using this knowledge to better approximate the true density and invert it to find the KS potential. An alternative is to use similar boundary conditions and numerical methods for both the interacting system and inverse problem with the understanding that there will be some cancellation of errors if the methods are similar enough. This is not really an inverse crime since the solution method for solving the interacting system is not identical to the numerical methods for solving the KS equations but it is not a requirement when solving these problems. In fact, the methods presented in Sec. 2.2 can be applied to target densities resulting from a variety of different numerical methods with no need to rely on error cancellations.

A.2 Target Densities with Numerical Noise

The examples in the preceding section assume a perfect knowledge of the target density but this is not a realistic assumption for most inverse problems. As mentioned in the previous example, the wave function numerical methods used in computing interacting densities are often very different from DFT numerical methods and the correct boundary conditions are not always known beforehand or may even be incorrect. It is far more realistic to assume that the target density in a density-to-potential inversion has some numerical error that must be accounted for to avoid overfitting. In order to illustrate the effects of numerical noise in the target density, we take the ground state orbital of the simple harmonic oscillator, add in a small amount of noise generated from the standard normal distribution $\operatorname{N}\!\left(\mu,\sigma^{2}\right)$, 31 and put the density corresponding to this modified orbital in the one-orbital formula. The results of this inversion in Fig. 16 are very poor especially in the asymptotic region. If we add a relative error to the ground state orbital instead of the uniform error described above then the resulting potential is still noisy but the noise level is now the same everywhere as seen in Fig. 17. These two examples show that the relative error of the target density sets a limit on the convergence in density-to-potential inversions. Figure 18 shows that simply increasing the number of grid points actually makes the inverted potential worse when the same amount of error is present in the target density.

A.3 Finite Difference Operators

As mentioned in the previous examples, higher-order finite difference operators don’t always improve density-to-potential inversions and boundary conditions often have a dramatic impact on the quality of an inversion. Another way to study the delicate balance between finite difference operators and boundary conditions is to use the ground state density of a particle in a box in the one-orbital formula. If we compute the derivative in Eq. (7) with one finite difference operator on the left side of the simulation box and a different operator on the right then the resulting second derivative has different error patterns on both sides of the box. The resulting potential has a jump between the left and right sides of the box as seen in Fig. 19. This jump is actually very small, (on the order of $10^{-4}$ as seen in the scale), but makes a significant impact on the error pattern. Although one usually doesn’t mix finite difference operators inside the simulation box, different operators are often used at the boundaries and produce similar jumps even if the operators are the same order of accuracy. An example of this jump at the boundary can be seen in Fig. 19 at the right boundary where a fourth-order backwards finite-difference operator produces a different error pattern than the fourth-order centered operator used elsewhere. (See Eq. (5) for the matrix representation of this fourth-order finite-difference operator.)

Many of the somewhat counterintuitive behaviors shown in the preceding examples can be understood mathematically by studying the effects of errors in the one-orbital formula. If we have the exact density of a one orbital system then the potential is given by the formula

[TABLE]

where $\phi_{i}=\sqrt{n_{i}}$ and $\varepsilon_{i}^{\text{fd}}$ is the error of the approximate second derivative at the point $x_{i}$. In the example shown in Fig. 19 there are different error patterns $\varepsilon_{i}^{\text{fd}}$ for the left- and right-hand sides of the box and this difference is magnified during the inversion because of the division by the square root of the density. A more realistic example includes some numerical error $\varepsilon_{i}^{n}$ in the density being inverted and the inversion formula becomes

[TABLE]

If the noise $\varepsilon_{i}^{n}$ is not smooth then its second derivative will be very noisy and dominate the other terms in Eq. (37). Figure 20 shows an example where a very small normally distributed error severely limits the accuracy of the potential.

A.4 Rounding Errors

Rounding errors play an important role in many of the inversions methods dependent upon numerical optimization methods. The Wu-Yang algorithm discussed in Sec. 2.2 needs to compute the kinetic energy at each iteration of a given numerical optimization routine. If the density of a system decays exponentially then the kinetic energy from the asymptotic regions is very small in comparison to values near the nuclei. Even if the values of the kinetic energy density can be computed accurately in all regions of a simulation, adding all of the values together to get the total kinetic energy will often result in rounding errors due to the large differences in magnitude involved. In Fig. 21 we apply the Wu-Yang inversion algorithm to the ground-state density of the harmonic oscillator on a grid with 51 points ranging from $-8$ to 8. Both the gradient-based TN algorithm42 and the non-gradient-based modified Powell’s method46 fail to resolve the density in the asymptotic region, ($\left|x\right|>5.5$), when used in the Wu-Yang inversion method. The abrupt cutoffs near $-5.5$ and 5.5 are a result of rounding errors when computing the kinetic energy with double-precision floating point numbers. If we use the midpoint integration rule to compute the kinetic energy in the region $\left|x\right|<5.5$ the result in double precision arithmetic is 0.2500000000000231 and the same method used in the region $\left|x\right|>5.5$ produces the value $-2.32460041955568\times 10^{-14}$. When the two numbers above are added together, most of the smaller number’s significant digits are lost and this loss of information limits the accuracy of the recovered potential in the asymptotic region. The same example performed using single-precision floating point numbers gives an even stricter cutoff near $-3.5$ and 3.5 as seen in Fig. 22. In this second example the kinetic energy in the region $\left|x\right|<3.5$ is 0.25000963 when using single-point precision numbers compared to a value of $-9.6271506\times 10^{-6}$ in the region $\left|x\right|>3.5$, which again results in a loss of significant figures when added together.

The preceding examples of time-independent density-to-potential inversions illustrate many of the pitfalls in performing numerical density-to-potential inversions. Simply refining the grid or increasing the accuracy of the kinetic energy operator does not guarantee an improved density-to-potential inversion and in many cases will worsen the solution. In cases where the error pattern is not known explicitly it is actually not correct to fit the density exactly as even small density errors can lead to large artifacts in the potential. Finally, inversions reliant on numerical optimization routines often suffer from rounding errors, which set yet another hard limit on the convergence of density-to-potential inversions.

Appendix B Discrete-Adjoint Method

In this appendix we derive the discrete adjoint equations for the PDE-constrained optimization method presented in Sec. 2.2.3. The discrete-adjoint method is a very efficient method for computing derivatives (discrete functional derivatives) in PDE-constrained problems.47 Computationally, these derivatives amount to a sparse linear solve for each occupied orbital and are about the same computational cost as computing the cost functional. The derivations give here are only valid for systems without degeneracies because this method is based on the algebraic method described in Ref. 48. A similar derivation can be performed for systems with degeneracies using the method of Ref. 49.

The notation used in the following derivations is the same notation used in Ref. 38 for TDDFT inverse problems. It is very similar to the continuous adjoint derivation found in quantum optimal control derivations applied to TDDFT.50 The discrete adjoint method is discussed in depth in Ref. 51. The notation used here is simply a modification of other discrete adjoint derivations applied to the DFT inverse problem.

The discretized cost functional

[TABLE]

can be differentiated with respect to an arbitrary orbital $\phi_{\bm{\alpha}}^{\gamma}$ to produce the expression

[TABLE]

where $\gamma$ is the orbital index and $\bm{\alpha}=\alpha_{x},\alpha_{y},\alpha_{z}$ is the spatial index of the orbital being differentiated. We introduce the Lagrangian

[TABLE]

where the Lagrange multipliers $\left\{\chi^{m}\right\}$ and $\left\{\lambda^{m}\right\}$ constrain the orbitals to satisfy Eqs. (4a) and (3). The Hamiltonian $H\left(v\right)=-\frac{\Delta}{2}+v$ in Eq. (40) is is a discretized version of Eq. (4a). We subtract the Lagrangian from the cost functional to form the total functional

[TABLE]

Differentiating the total functional $J$ with respect to the adjoint orbitals $\left\{\chi^{m}\right\}$ and normalization multipliers $\left\{\lambda^{m}\right\}$ gives the standard KS equations

[TABLE]

The adjoint equations are found by differentiating with respect to the KS orbitals and eigenvalues to produce a set of sparse linear equations

[TABLE]

where $\frac{\partial F}{\partial\varepsilon^{\gamma}}=0$. After solving the direct and adjoint equations in Eqs. (42)-(44), the total derivative is given by the expression

[TABLE]

If the potential is placed diagonally along $H$ then this last expression reduces to

[TABLE]

A very similar derivation is found when using the scaled orbitals described in Sec. 2.2.3. In this case, the discretized cost functional is

[TABLE]

where $\phi_{i}^{m}=\sqrt{\tilde{n}_{{\bf i}}}g_{{\bf i}}^{m}$ and

[TABLE]

The forward equations are

[TABLE]

The corresponding adjoint equation is

[TABLE]

and the total derivative is

[TABLE]

where $H$ is a discretization of the scaled Hamiltonian given in Eq. (21).

The discrete adjoint equations for a periodic system follow the same pattern as the above derivations but with the added complication of a numerical integration in $k$-space. The discretized cost functional is

[TABLE]

where $u_{{\bf j}}=\frac{1}{\Omega_{{\bf k}}}dk_{{\bf j}}$ is the $j^{\text{th}}$ integration weight in $k$-space divided by the volume of the first Brillouin zone and $\phi_{{\bf i},{\bf j}}^{m}$ corresponds to the point $\phi_{{\bf k}_{{\bf j}}}^{m}\left({\bf x}_{{\bf i}}\right)$. The forward equations are

[TABLE]

where $H_{{\bf j}}\!\left(v\right)$ is a discretization of the Hamiltonian $-\frac{1}{2}\left(\Delta+2\imath{\bf k}_{{\bf j}}\cdot\bm{\nabla}-k_{{\bf j}}^{2}\right)+v\!\left({\bf r}\right)$. The corresponding adjoint equations are

[TABLE]

and the total derivative is

[TABLE]