Exact exchange-correlation potentials from ground-state electron densities

The quest for accurate exchange-correlation functionals has long remained a grand challenge in density functional theory (DFT), as it describes the many-electron quantum mechanical behavior through a computationally tractable quantity—the electron density—without resorting to multi-electron wave functions. The inverse DFT problem of mapping the ground-state density to its exchange-correlation potential is instrumental in aiding functional development in DFT. However, the lack of an accurate and systematically convergent approach has left the problem unresolved, heretofore. This work presents a numerically robust and accurate scheme to evaluate the exact exchange-correlation potentials from correlated ab-initio densities. We cast the inverse DFT problem as a constrained optimization problem and employ a finite-element basis—a systematically convergent and complete basis—to discretize the problem. We demonstrate the accuracy and efficacy of our approach for both weakly and strongly correlated molecular systems, including up to 58 electrons, showing relevance to realistic polyatomic molecules.

from electron densities for small systems does not pose much difficulty. If the input density is generated using Slater-type basis functions or a very dense numerical grid, then the methods of Refs. 12, 20, and 32 would produce at least as good results as the proposed method in Fig. 2 (which, by the way, does not show in detail what happens at very small and very large r).
2. For the same accurate potential, the L2 norm density error used in this work is orders of magnitude smaller than the L1 norm used in other methods such as that of Ref. 32. The convergence thresholds of 10**(-5)e**2 and 10**(-6)e**2 used in the paper are probably equivalent to the L1 norm error of ~0.01e, which is much greater than the values used in Refs. 32 and elsewhere. It is very misleading to claim "unprecedented" L2-norm errors when most readers familiar with the problem are thinking in terms of L1 norms. 3. The authors attempt to use Gaussian basis-set densities to produce xc potentials that do not oscillate and decay as -1/r. Such attempts are misguided because Kohn-Sham potentials corresponding to Gaussian densities in the basis-set limit must oscillate and increase as r**2 at large r, as shown in Refs. 31 and 32. Any potential from a Gaussian density that does not oscillate (at least near the nucleus) or does not blow up is not the exact answer to the inversion problem. 4. The proposed method gets rid of the oscillations by "adding a small correction, Delta rho(r), to rho_data(r), so as to correct for the missing cusp at the nuclei" (eq. 6). In other words, the "solution" is to fudge the input density. This is hardly "an advance that is likely to influence thinking in the field". Moreover, the principle of this correction is clearly inspired by the trick of Ref. 33, a fact that is not acknowledged. 5. The authors apply the following adjectives to describe their method: "rigorous", "unambiguous", "unprecedented", etc. None of these characterizations is justified. Modifying an input density is neither "rigorous" nor "unambiguous". Convergence with respect to the L2 norm of 10**(-6)e**2 is anything but "unprecedented", especially after rho_data is "corrected". Wobbly potentials seen in Figs. 5 and 6c are poorer than could be obtained by other methods. 6. The general idea of using constrained optimization for mapping densities to xc potentials is not new. Refs. 12, 13, 18, 20 employ similar principles.
In summary, the major claims of the paper are not novel, the conclusions are not original, and the work is not convincing overall. Therefore, it does not merit publication.
Reviewer #3 (Remarks to the Author): In "Exact exchange-correlation potentials from ground-state electron densities" the authors present an algorithm to solve the 'inversion' problem in DFT. Where given an electron density, one can obtain the Kohn-Sham potential that will result in the density. This has been an open problem in the field for many years, and an efficient solution to this problem is an important step for a variety of different fields. The paper is well written and clear. Therefore, I believe this work should be considered for publication in Nature Communications; however, I believe there should be several points to address before the manuscript is published.
1.As the authors note, partial differential equation constrained optimization has been explored by other authors before but only in model/test systems in 1-dimension. In many ways, the difficulty in the inverse problem is the numerical solution. The Wu and Yang approach with regularization should provide an 'optimal' solution and in relatively simple systems one can numerically find that solution. However, when one tries to apply it to large systems, then it becomes numerical intractable. Therefore, my largest concern about this work, is that this approach works great for small systems, but is still numerically intractable to any sizeable system. The largest systems presented here are 10 electrons. I think to showcase that this is a significant step forward, they authors should provide numerical results for a large system (50-100 electrons). If one can still solve the inverse problem using this partial differential equation constrained optimization methodology, it will really show that this work is a breakthrough in this field. 2.As stated above, it is the numerical difficulty which really makes this problem hard. Therefore, how is the convergence of these equations? Specifically, on convergence of ||\rho_data -\rho||_{L^2}: (a) How many cycles does the density take to converge? (b) How is the convergence affected by the initial Vxc guess? (c) How is the convergence affected by the number of electrons in the system? (d) How is convergence affected by the weak, eg H2(eq), and strong correlation, eg H2(d), systems? 3.In the H2 stretched/dissociated case, what is the orbital energy of the LUMO going to? Are the significant difficulties in solving the strong correlation system with multiple near-degenerate orbitals? Naively, I would have expected that during the optimization cycles you could have potentials, which make the original LUMO become the HOMO as you try to force the strongcorrelation solution. Does the optimization process continue to be well behaved?
Again, I want to emphasize that I believe this is important work. However, there have been several inverse methods that have been proposed over the last 15 years that have been applied to small systems that work relatively well. The difficulty in the problem is in the numerics and I think the authors should discuss and show that this optimization scheme is numerically tractable and simple to be able to deal with large and complex systems.

Response to Reviewers
Bikash Kanungo, Paul M. Zimmerman, Vikram Gavini We thank the reviewers for the careful review of the manuscript, and their comments and suggestions to improve the manuscript. Below we are providing a detailed response to the reviewers' comments and the changes made to the manuscript. We hope we have satisfactorily addressed all the questions and concerns of the reviewers.

Response to Reviewer #1
1. Solving the inverse DFT problem, i.e. obtaining (numerically) exact exchange-correlation potentials from electron densities, is one of the most fundamental and at the same time numerically challenging problems in density functional theory. As a consequence, it has attracted great research interest over the past decades. While several important steps forward have been made, essentially all approaches in use so far suffer from numerical issues, which prevent a robust and systematically convergent solution of this problem.
The great value and novelty of the present work lies in the fact that it properly and simultaneously addresses all the numerical challenges preventing an accurate inverse solution so far. This includes, in particular, (i) input electron densities of very high quality (incremental Full Configuration Interaction, iFCI), (ii) a formulation of the potential evaluation as a PDE-constrained optimization problem, which (iii) is solved in a finite-element basis that can be systematically converged to completeness, (iv) an initial guess for the XC potential with correct asymptotic behavior, which is preserved 1 in the calculation through boundary conditions, and (v) a cusp-correction for the input electron densities, which is necessary since the iFCI densities are represented in a Gaussian basis set.
Response: We thank the reviewer for commenting on the fundamental as well as challenging nature of the inverse DFT problem, and for summarizing the value and novelty of this work in solving this problem. We remark that this limitation, in the context of systems with long-range static correlation, is a consequence of the wrong far-field decay of Gaussian basis-set densities, and not a shortcoming of the inverse procedure itself. We expect better accuracy for systems with such long-range (static) correlation with the availability of more accurate densities. 2 4. The introduction implies that this approach "unlocks the door to ... highly accurate exchange-correlation functionals that provide precise energies". However, the road from potentials for specific densities to general density functionals is not obvious, and the discussion remains rather general and vague about this point. Can the authors briefly elaborate on the question of how their approach can specifically help to obtain accurate energy functionals from the potentials? Or should the statement in the introduction be adapted?
Response: We have added a brief discussion following this sentence in the revised manuscript. The path we envision involves generating {ρ (i) (r), v 6. How were the reference ionization potentials obtained (Table 1)? 3 Response: The reference ionization potentials were obtained using the same method as for the reference electron densities (incremental full CI). Specifically, one electron was removed from the system and the total energy recomputed, the ionization energy being the difference between the original energy and the minus-one-electron energy.
This point has been noted in the methods section.
Response to Reviewer #2 Response: The completeness of the basis-in which the Kohn-Sham orbitals (ψ i (r)) as well as the exchange-correlation potential (v xc (r)) are discretizedis of utmost importance in an inverse DFT calculation, irrespective of the basis-set used in generating ρ data . It has been highlighted, both theoretically and numerically, that the lack of completeness in the basis-set employed in discretizing ψ i 's and v xc can result in an ill-posed problem, thereby producing non-unique solutions and/or spurious oscillations in v xc [4,5,6]. Thus, the use of finite-element basis adopted in this work, which can systematically be converged to completeness, is crucial to resolving the numerical challenges associated with the inverse DFT problem. However, we recognize the reviewer's concern that the range of systems demonstrated in the previous manuscript do not show the full power of our methodology. To address this concern, we have added two larger molecules to the revised manuscript and confirmed that accurate xc potentials can still be straightforwardly constructed using our method. The 1,3,-dimethylbenzene (C 8 H 10 ) and ortho-benzyne (C 6 H 4 ) molecules have 58 and 40 electrons, respectively, which are unprecedented for inverse DFT calculations. The ease in generating xc potentials for these molecules confirms the advances presented herein.
For the atomic systems considered in our verification studies (cf. Fig. 2 Response: While the v xc corresponding to Gaussian basis-set densities exhibit oscillations and r 2 growth at large r, these are precisely Gaussian-basis set artifacts and not the features of the physical v xc . The focus of the present work is not to reproduce the Gaussian basis-set artifacts (which we demonstrate in Figure 4 of the revised manuscript for H 2 ), but rather to use Gaussian basis-set densities (as they are efficient for CI calculations) and still be able to evaluate the physically exact v xc . To this end, our approach of using ∆ρ cusp-correction, as well as the far-field boundary condition ensuring correct asymptotics for v xc , provide a robust means to correct for the Gaussian basis-set artifacts and compute the physically exact v xc .
4. The proposed method gets rid of the oscillations by "adding a small correction, ∆ρ, to ρ data , so as to correct for the missing cusp at the nuclei" (eq. 6). In other words, the "solution" is to fudge the input density. This is hardly "an advance that is likely to influence thinking in the field". Moreover, the principle of this correction is clearly inspired by the trick of Ref. 33, a fact that is not acknowledged.
Response: We note that the ∆ρ-correction 'corrects' for the missing cusp in the Gaussian input density, and is far from being an ad-hoc fudging of the input density.
We hypothesize that the ∆ρ correction encapsulates the Gaussian basis-set error near the nuclei. Evidence to this hypothesis is provided in Fig. 5 of the revised manuscript, wherein we observe only weak sensitivity of the inverted v xc to the xc functional used in constructing ∆ρ. This, in our view, constitutes a robust numerical strategy to mitigate the wild oscillations induced by the missing cusp at the nuclei in Gaussian densities.
We agree that there is similarity of our ∆ρ correction to that of ∆v osc in [9]. We have acknowledged this on page 8 of the revised manuscript.

5.
The authors apply the following adjectives to describe their method: "rigorous", "unambiguous", "unprecedented", etc. None of these characterizations is justified. Modifying an input density is neither "rigorous" nor "unambiguous". Convergence with respect to the L2 norm of 10 −6 e 2 is anything but "unprecedented", especially after ρ data is "corrected". Wobbly potentials seen in Figs. 5 and 6c are poorer than could be obtained by other methods.
Response: We first note that we did not use "rigorous" to describe our work. We refer to the responses in #2,#3 and #4 for more detailed justification for the 'unambiguous' solution to the inverse DFT problem made possible by this work, and 'unprecedented' accuracy obtained. Briefly, we qualify our approach as both "robust" and "unambiguous" owing to the stability of the algorithm, non-sensitivity to the initial guess, and above all the uniqueness of the converged v xc . Furthermore, we qualify our work as "unprecedented" on the basis of attaining L 2 and L 1 errors of O(10 −5 − 10 −4 ) in the density for 3D problems, as well as having an excellent agreement between H and −I p .
We note that the slight wobbly nature of the v xc only appears in the case of the dissociated hydrogen molecule (H 2 (d)), and not in other molecules. We emphasize that H 2 (d) is a particularly difficult system for inverse DFT, and this work is the first to attempt to accurately evaluate the exact v xc for a dissociated molecule. We suspect the slight wobbliness to be a consequence of the Gaussian basis-set density. However, given the difficulty in solving an inverse problem for a dissociated system, we consider the quality of our v xc to be quite good. 6. The general idea of using constrained optimization for mapping densities to xc potentials is not new. Refs. 12, 13, 18, 20 employ similar principles.
Response: We agree that the idea of constrained optimization, in the context of inverse DFT, has been explored in the past, and we do not claim the PDE-constrained optimization in itself as the novelty of this work. As noted by Reviewer 1, the novelty of the work lies in combining the following to solve the outstanding inverse DFT problem: PDE-constrained optimization, use of finite-elements that can be systematically converged to completeness, the cusp-correction in Gaussian densities, and use of far-field boundary conditions that ensures the correct asymptotic behavior for v xc .
In total, the proposed technique provides the ability to invert xc potentials for molecules that are far larger than ever before (40 or more electrons), demonstrating a systematic means to determine v xc for systems well beyond the scope of prior methods.
Response to Reviewer #3 1. In "Exact exchange-correlation potentials from ground-state electron densities" the authors present an algorithm to solve the 'inversion' problem in DFT. Where given an electron density, one can obtain the Kohn-Sham potential that will result in the density.
This has been an open problem in the field for many years, and an efficient solution to this problem is an important step for a variety of different fields. The paper is well written and clear. Therefore, I believe this work should be considered for publication in Nature Communications; however, I believe there should be several points to address before the manuscript is published.
Response: We thank the reviewer for commenting on the importance of this problem, and for being cautiously positive on the manuscript. 8 2. As the authors note, partial differential equation constrained optimization has been explored by other authors before but only in model/test systems in 1-dimension. In many ways, the difficulty in the inverse problem is the numerical solution. The Wu and Yang approach with regularization should provide an 'optimal' solution and in relatively simple systems one can numerically find that solution. However, when one tries to apply it to large systems, then it becomes numerical intractable. Therefore, my largest concern about this work, is that this approach works great for small systems, but is still numerically intractable to any sizeable system. The largest systems presented here are 10 electrons. I think to showcase that this is a significant step forward, they authors should provide numerical results for a large system (50-100 electrons). If one can still solve the inverse problem using this partial differential equation constrained optimization methodology, it will really show that this work is a breakthrough in this field.
Response: We thank the referee for suggesting that we demonstrate our approach on larger systems. In the revised manuscript, we have demonstrated the efficacy of our approach using two additional systems-(a) LDA density based verification for 1,3-dimethylbenzene (C 8 H 10 , 58 electrons), (b) CI density based exact v xc evaluation for ortho-benzyne radical (C 6 H 4 , 40 electrons), a strongly correlated system. In both the cases, we obtain similar accuracies as the smaller systems reported in this work.
Accurately solving the inverse DFT problem on such large systems is unprecedented.
We hope these additional results demonstrates the robustness and efficacy of the proposed approach.
3. As stated above, it is the numerical difficulty which really makes this problem hard.
Therefore, how is the convergence of these equations? Specifically, on convergence of ||ρ data − ρ|| L 2 : (a) How many cycles does the density take to converge? (b) How is the convergence affected by the initial Vxc guess? (c) How is the convergence affected by the number of electrons in the system? (d) How is convergence affected by the weak, eg H2(eq), and strong correlation, eg H2(d), systems? Response: (a) For all the systems, except H 2 (d), it took 300-500 BFGS iterations to obtain an accuracy of ||ρ data − ρ|| L 2 < 10 −4 .
(b) We have experimented with two different initial guess-PW92 (LDA) potential and LB94 (GGA) potential. While the final result is independent of the initial Response: We thank the reviewer for being cautiously positive on this work. We hope the additional results on C 8 H 10 and C 6 H 4 (strongly correlated radical) demonstrates the efficacy of the proposed approach to deal with large complex systems.

Additional Changes
As part of this revision, we conducted inverse DFT simulations with better FE discretization for the atomic systems in the verification study (results in Figure 2  [ρ data ] is also evident from Figure 2 in SI, which shows the data on a log scale to demonstrate the near-field and far-field agreement. The minor revisions and clarifications made by the authors do not change my assessment: the proposed method is neither more robust nor more accurate than many existing approaches to DFT inversion. For non-Gaussian densities, any careful numerical inversion procedure would give as good potentials as those of Fig. 2. For Gaussian densities, the potentials of this work are still imperfect (crude and wobbly) -see Figs. 6 and 7. Moreover, what makes the method work passably for Gaussian densities is not the proposed constrained minimization of Eq. 5 but an expedient alteration of the input density.
The authors are billing their method as "an advance that overcomes all the aforementioned open issues in the inverse DFT problem", "a novel finite-element basis approach that is systematically convergent and complete, and thereby eliminates the ill-conditioning associated with previous approaches". However, they do not provide a single direct comparison with any of the previous approaches to substantiate those claims. In the absence of such evidence, the claims are not ready for publication.
Reviewer #3 (Remarks to the Author): In the revised manuscript the authors have addressed the concerns I brought up in the original manuscript. The application to significantly larger systems shows the robustness of the method and the ability to perform these inverse calculations on complicated systems that will be of interest to the community.
Given the number of years people have worked on this problem and these authors have provided a relatively simple to implement and accurate method, I recommend publishing this manuscript.
My only comment to the authors about improving the manuscript would be to include the information about convergence that was provided in the response to the reviewer. I found this information very helpful. Particularly the discussion about how the HOMO-LUMO gap affects the convergence. I think other readers would enjoy this information in the main text of SI.
time-consuming task. More importantly, the performance of these approaches and the outstanding issues are well documented in the literature.
Thus, in the revised manuscript, we have used the data and results from the literature to compare and contrast existing methods for DFT inversion with our proposed approach, in terms of accuracy, robustness and computational viability. A new section in the SI (Section 5) is devoted to this comparison. Further, we have made appropriate revisions to ensure that all the claims in the revised manuscript are supported by our objective findings.

Response to Reviewer #3
In the revised manuscript the authors have addressed the concerns I brought up in the original manuscript. The application to significantly larger systems shows the robustness of the method and the ability to perform these inverse calculations on complicated systems that will be of interest to the community. Given the number of years people have worked on this problem and these authors have provided a relatively simple to implement and accurate method, I recommend publishing this manuscript.
Response: We thank the reviewer for the careful review of the manuscript. We are glad that we have addressed all the reviewer's concerns satisfactorily.
My only comment to the authors about improving the manuscript would be to include the information about convergence that was provided in the response to the reviewer. I found this information very helpful. Particularly the discussion about how the HOMO-LUMO gap affects the convergence. I think other readers would enjoy this information in the main text of SI.
Response: In the revised manuscript, we have added a section in the SI discussing the convergence aspects, and its dependence on the HOMO-LUMO gap.