Kohn-Sham Inversion with Mathematical Guarantees

Basic Information

• Paper ID: 2409.04372
• Title: Kohn-Sham inversion with mathematical guarantees
• Authors: Michael F. Herbst (EPFL), Vebjørn H. Bakkestuen (Oslo Metropolitan University), Andre Laestadius (Oslo Metropolitan University & University of Oslo)
• Classification: physics.chem-ph, math-ph, math.MP, quant-ph
• Publication Date: September 2024 (arXiv v3: May 5, 2025)
• Paper Link: https://arxiv.org/abs/2409.04372

Abstract

This paper employs the exact Moreau-Yosida regularization method to obtain the exchange-correlation potential (xc potential) for periodic systems. The research reveals profound connections between rigorous mathematical principles and efficient numerical implementation, representing the first application of Moreau-Yosida-based inversion methods to actual physical systems. The authors developed a mathematically rigorous inversion algorithm and validated it on representative bulk materials (bulk silicon, gallium arsenide, and potassium chloride). The algorithm constructs rigorous error bounds and provides numerical verification, opening new avenues for analyzing Kohn-Sham inversion methods and potentially advancing mathematical approaches for developing approximate functionals.

Research Background and Motivation

Core Problems

Density functional theory (DFT) is an indispensable tool in chemistry, materials science, and solid-state physics. In the Kohn-Sham (KS) formulation, all unknowns in DFT are concentrated in the exchange-correlation (xc) functional, which requires approximation. Although DFT is exact in principle, KS-DFT faces challenges in certain physical scenarios, particularly:

1. Fractional Charge Problem: Difficulty in accurately describing processes involving fractional electronic charge (e.g., dissociation or charge-transfer excitations)
2. Band Gap Problem: Semiconductor band gaps are systematically underestimated

Research Significance

1. Theoretical Deficiency: Lack of mathematical understanding between the exact universal density functional and commonly used approximations, making rigorous construction of new, better functionals difficult
2. Inversion Problem: KS inversion (determining the exact xc potential from a given ground-state density) is far less studied than the forward KS-DFT problem
3. Mathematical Connections: Earlier work established profound links between fractional charge and band gap problems and the non-differentiability of KS-DFT

Limitations of Existing Methods

• Existing KS inversion methods (e.g., van Leeuwen-Baerends, Zhao-Morrison-Parr, Wu-Yang) lack rigorous mathematical guarantees
• Absence of robust and efficient numerical schemes
• No rigorous error bounds and convergence analysis

Research Motivation

Recent theoretical results 49 demonstrated that the xc potential can be obtained through the mathematical limit of Moreau-Yosida (MY) regularization of the exact universal functional. MY regularization addresses the non-differentiability of the exact universal functional, which is closely related to the aforementioned physical problems. This paper represents the first application of this theoretical framework to actual physical systems.

Core Contributions

1. First Implementation: First realization of MY framework-based KS inversion on actual physical systems
2. Mathematical Rigor: Development of an inversion algorithm with rigorous mathematical guarantees, establishing exact inversion formulas (Equation 10)
3. Error Bounds: First establishment of rigorous error bounds for the inverse KS problem (Equations 14, 16, 17) with numerical verification
4. Non-expansivity Proof: Proof of the (firmly) non-expansive property of the proximal mapping (Equation 12)
5. Practical Application: Successful validation of the algorithm on three representative bulk materials (Si, GaAs, KCl)
6. Theory-Practice Bridge: Establishment of connections between rigorous mathematical theory and practical numerical implementation

Methodology Details

Task Definition

Input: Exact ground-state density $\rho_{gs}$ (possibly from experimental data, full configuration interaction, coupled cluster, or quantum Monte Carlo calculations)

Output: Corresponding exchange-correlation potential $v_{xc}$ that reproduces this density in an auxiliary non-interacting system

Constraint: Assumption that the density is non-interacting v-representable (i.e., there exists a potential such that $\rho_{gs}$ is also a non-interacting ground-state density)

Mathematical Framework

1. Function Space Selection

• Density Space: $D = H^{-1}_{per}(\Omega, \mathbb{C})$ (periodic Sobolev space)
• Potential Space: $V = H^1_{per}(\Omega, \mathbb{C})$ (dual space of D)
• Norm Definition: $\|u\|^2_{H^s_{per}} = \sum_G (1 + |G|^2)^s |\hat{u}_G|^2$ where G denotes reciprocal lattice vectors

2. Duality Mapping

The duality mapping $J: D \to V$ is defined as: $J(\rho) = \{v \in V : \|v\|^2_V = \|\rho\|^2_D = \langle v, \rho \rangle\}$

Under the chosen function spaces, the duality mapping has an explicit form: $J[\rho](r) = (\Phi * \rho)(r) = \int_{\mathbb{R}^3} \frac{\rho(r')}{4\pi|r-r'|} e^{-|r-r'|} d^3r'$ where $\Phi(r) = e^{-|r|}/(4\pi|r|)$ is the Yukawa potential, a form that is numerically tractable.

3. Moreau-Yosida Regularization

The guiding density functional is defined as: $\mathcal{F}(\rho) = T(\rho) + E_H(\rho) + \int_\Omega v_{ext}\rho$ where $T(\rho)$ is the kinetic energy functional and $E_H(\rho)$ is the Hartree contribution.

The key optimization problem: $\mathcal{E}(\rho; \rho_{gs}) = \mathcal{F}(\rho) + \frac{1}{2\varepsilon}\|\rho - \rho_{gs}\|^2_D$

Minimizing this functional yields the proximal density $\rho^\varepsilon_{gs} = \arg\min_\rho \mathcal{E}(\rho; \rho_{gs})$

4. Inversion Formula

The exchange-correlation potential is obtained through the following limit: $v_{xc}(r) = \lim_{\varepsilon \to 0^+} \frac{1}{\varepsilon} \int_{\mathbb{R}^3} \frac{\rho^\varepsilon_{gs}(r') - \rho_{gs}(r')}{4\pi|r-r'|} e^{-|r-r'|} d^3r'$

This is the core theoretical result of the paper, providing an explicit computational formula from the proximal density to the xc potential.

Model Architecture

Numerical Implementation Strategy

1. Orbital Parameterization: Parameterize density using orthonormal orbitals $\Phi = (\psi_1, \ldots, \psi_{N_b})$: $\rho_\Phi(r) = 2\sum_{i=1}^{N_b} |\psi_i(r)|^2$
2. Energy Expression (Equation 15): $\mathcal{E}(\Phi, \rho_{gs}) = \sum_{i=1}^{N_b} \int_\Omega |\nabla\psi_i|^2 + E_H(\rho_\Phi) + \int_\Omega v_{ext}\rho_\Phi + \frac{1}{2\varepsilon}\|\rho_\Phi - \rho_{gs}\|^2_D$
3. Optimization Method:
   • Quasi-Newton scheme based on BFGS
   • Adaptation to Stiefel manifold geometry (maintaining orbital orthogonality)
   • Stopping criterion: optimizer reaches machine precision or iterative change in $\rho^\varepsilon_{gs}$ falls below 0.01ε
4. ε Sequence: Exponentially decreasing sequence ranging from 1 to approximately $10^{-7}$

Technical Innovations

1. Non-expansivity Proof

Proof that the proximal mapping $\rho \mapsto \rho^\varepsilon$ is a (firmly) non-expansive operator: $\|\rho^\varepsilon - \tilde{\rho}^\varepsilon\|_D \leq \|\rho - \tilde{\rho}\|_D$

Proof Strategy:

• Utilize $-\frac{1}{\varepsilon}J(\rho^\varepsilon - \rho) \in \partial\mathcal{F}(\rho^\varepsilon)$
• Apply maximal monotonicity of subdifferentials
• Apply Hölder inequality

2. Error Bound Theory

Define the ratio $Q_\varepsilon(\Delta\rho) := \frac{\|\rho^\varepsilon_{gs} - \tilde{\rho}^\varepsilon_{gs}\|_D}{\|\Delta\rho\|_D} \leq 1$

Primary Error Bound (Equation 14): $\|v^\varepsilon_{xc} - \tilde{v}^\varepsilon_{xc}\|_V \leq \frac{1 + Q_\varepsilon(\Delta\rho)}{\varepsilon}\|\Delta\rho\|_D$

Refined Bound (Equation 16): $\left\|v^\varepsilon_{xc} - \tilde{v}^\varepsilon_{xc} - \frac{1}{\varepsilon}J(\Delta\rho)\right\|_V \leq \frac{Q_\varepsilon(\Delta\rho)}{\varepsilon}\|\Delta\rho\|_D$

Introducing ratios $R_\varepsilon$ and $S_\varepsilon$, the paper proves (Equation 17): $0 \leq 1 - Q_\varepsilon(\Delta\rho) \leq R_\varepsilon(\Delta\rho) \leq 1 + Q_\varepsilon(\Delta\rho) \leq 2$

3. Distinction from Traditional Methods

• Traditional Methods: Lack rigorous mathematical guarantees, typically based on heuristic optimization
• This Paper's Method:
  • Based on convex analysis and Banach space theory
  • Provides convergence guarantees ($\rho^\varepsilon_{gs} \to \rho_{gs}$ as $\varepsilon \to 0^+$)
  • Computable error bounds
  • Handles non-differentiability of functionals

Experimental Setup

Material Systems

Three representative bulk materials were studied:

1. Bulk Silicon (Si): Typical semiconductor
2. Gallium Arsenide (GaAs): Compound semiconductor
3. Potassium Chloride (KCl): Ionic crystal

Computational Details

Forward Calculations (Reference)

• xc Functional: PBE functional
• Pseudopotentials: PBE pseudodojo standard pseudopotentials (with nonlinear core corrections)
• k-point Spacing: Maximum 0.12 Å$^{-1}$
• Kinetic Energy Cutoff: Approximately twice the recommended value (ensuring high accuracy)
• Software: Density-Functional ToolKit (DFTK)

Inversion Calculations

• Same pseudopotential approximation used (including Kleiman-Bylander nonlocal potential terms)
• ε Sequence: Exponentially decreasing from 1 to approximately $10^{-7}$
• Optimization stopping criterion: Machine precision or $\Delta\rho^\varepsilon_{gs} < 0.01\varepsilon$

Error Introduction Method

To test error bounds, controlled perturbations $\Delta\rho$ were introduced through Fourier basis truncation:

• Different cutoff energies $E_{cut}$ (15, 25, 35, 45 Ha)
• $E_{cut} = 45$ Ha as unperturbed reference
• Calculation of corresponding $\|\Delta\rho\|_D$

Experimental Setup Remarks

The authors acknowledge the use of an "inverse crime" setup (forward and inversion using the same model and discretization basis), but emphasize this is to:

1. Verify the rigor of mathematical theory
2. Enable direct comparison of inverted densities and potentials with reference values
3. Future work will employ reference densities from other high-accuracy methods

Experimental Results

Main Results: Exact Inversion (Noise-free)

1. Bulk Silicon (Figure 2)

• Potential Recovery: Potential plotted along crystal high-symmetry paths (O → (001) → O' → (110) → O'' → (111) → O)
• Convergence Performance:
  • $\varepsilon \sim 10^{-6}$: Relative error below 10%
  • One order of magnitude decrease in ε: Error decreases by another order of magnitude
• Spatial Features: Near the sharpest features of the potential, pointwise convergence is slower with larger relative errors

2. Gallium Arsenide (Figure 3)

• Potential plotted along similar paths (starting between Ga-Ga bonds)
• Absolute relative errors slightly larger than silicon at the same ε values
• Overall, reference potential is accurately recovered

3. Potassium Chloride (Figure 4)

• Path starting from potassium (K) atom
• Error characteristics similar to GaAs
• Reference potential successfully recovered for all three materials

Key Finding: Without additional noise ($\Delta\rho = 0$), the algorithm accurately recovers the xc potential for all three materials, validating method effectiveness.

Noisy Inversion and Error Estimation

1. Convergence Behavior (Figure 5)

Investigation of the effect of different basis truncations on convergence:

• Key Observation: As long as $\varepsilon > \|\Delta\rho\|_{L^2_{per}}$, potential convergence properties remain unchanged
• For smaller ε, potential begins to deviate from reference (in V-norm)
• Different truncation energies (15, 25, 35 Ha) correspond to different $\|\Delta\rho\|_D$

2. Non-expansivity Verification (Figure 6)

Computation of ratio $Q_\varepsilon(\Delta\rho) = \|\rho^\varepsilon_{gs} - \tilde{\rho}^\varepsilon_{gs}\|_D / \|\Delta\rho\|_D$:

• Theoretical Bound: $0 \leq Q_\varepsilon \leq 1$ (guaranteed by non-expansivity of proximal mapping)
• Numerical Results:
  • Large ε values: $Q_\varepsilon \ll 1$
  • $\varepsilon \to 0^+$: $Q_\varepsilon \to 1^-$
  • Perfect agreement with theoretical predictions

3. Error Bound Verification (Figure 7)

Ratio $S_\varepsilon$ (Figure 7 top):

• Definition: $S_\varepsilon(\Delta\rho) := \varepsilon \|v^\varepsilon_{xc} - \tilde{v}^\varepsilon_{xc} - \frac{1}{\varepsilon}J(\Delta\rho)\|_V / \|\Delta\rho\|_D$
• Theoretical bound: $0 \leq S_\varepsilon \leq Q_\varepsilon \leq 1$
• Numerical performance:
  • Large ε: $S_\varepsilon$ very small (slight deviation from bound due to differences of nearly zero quantities)
  • Small ε: Perfect agreement with bounds set by $Q_\varepsilon$
  • $\varepsilon \to 0^+$: $S_\varepsilon \to 1$

Ratio $R_\varepsilon$ (Figure 7 bottom):

• Definition: $R_\varepsilon(\Delta\rho) := \varepsilon \|v^\varepsilon_{xc} - \tilde{v}^\varepsilon_{xc}\|_V / \|\Delta\rho\|_D$
• Theoretical bound: $1 - Q_\varepsilon \leq R_\varepsilon \leq 1 + Q_\varepsilon$
• Numerical performance:
  • Strictly obeys bounds set by $Q_\varepsilon$
  • Closely follows lower bound $R_\varepsilon \geq 1 - Q_\varepsilon$
  • Large ε values: $R_\varepsilon \approx 1$
  • Small ε: Approaches lower bound

Numerical Challenge: For $\varepsilon \lesssim 5 \times 10^{-6}$, the problem becomes numerically challenging, manifesting as small oscillations in the trends of two quantities.

Summary of Experimental Findings

1. Theory Verification: Numerical calculations perfectly align with theoretical predictions of error bounds and non-expansivity
2. Robustness: Algorithm shows good robustness to density perturbations (within $\varepsilon > \|\Delta\rho\|$ range)
3. Universality Indication: $Q_\varepsilon$ may be estimated by a constant independent of $\Delta\rho$ (though parameters depend on ε and guiding functional)
4. Practical Applicability: Successful application to three different types of material systems
5. Numerical Precision: Maintains stable computation even at $\varepsilon \sim 10^{-7}$

Related Work

Development History of KS Inversion Methods

Early Work (1988-1994)

• Aryasetiawan & Stott (1988): Effective potential method
• Knorr & Godby (1992): Quantum Monte Carlo study of model semiconductors
• Görling (1992): Determination of KS potential and wave functions from electron density
• van Leeuwen & Baerends (1994): xc potential with correct asymptotic behavior

Optimization Methods (2002-2007)

• Wu-Yang (2002, 2003): Direct optimization method
• Zhao-Morrison-Parr (1994): ZMP method
• Bulat et al. (2007): Optimized effective potential in finite basis sets

Recent Progress (2013-2024)

• Software Development:
  • n2v (Shi, Chávez, Wasserman, 2022)
  • KS-pies (Nam et al., 2021)
• Solid-State System Extensions:
  • Aouina et al. (2023): Exact KS auxiliary systems for solid ground-state densities
  • Ravindran et al. (2024): Density inversion of local xc potential in solids
• Theoretical Analysis:
  • Burke group: Density-driven error analysis
  • Gould (2023): "Lieb-response" method

Mathematical Theory Foundation

Application of Moreau-Yosida Regularization in DFT

• Kvaal et al. (2014): Differentiable yet exact DFT form
• Laestadius et al. (2018, 2019): Generalized KS iteration on Banach spaces
• Penz et al. (2019): Guaranteed convergence of regularized KS iteration in finite dimensions
• Penz, Csirik, Laestadius (2023): Density-potential inversion from MY regularization (direct theoretical basis of this paper)

Relative Advantages of This Paper

1. First Implementation: MY framework first applied to actual physical systems
2. Rigorous Guarantees: Provides mathematically rigorous error bounds (unprecedented)
3. Theory-Practice Integration: Transforms abstract mathematical theory into computable numerical schemes
4. General Framework: Applicable to periodic systems, extensible to more complex systems
5. Error Analysis: Goes beyond heuristic error estimates of existing methods

Conclusions and Discussion

Main Conclusions

1. Method Effectiveness: Successfully developed and validated a rigorous KS inversion algorithm based on MY regularization
2. Theoretical Contributions:
   • Establishment of explicit inversion formula (Equation 10)
   • Proof of non-expansivity of proximal mapping (Equation 12)
   • Derivation of first rigorous error bounds (Equations 14, 16, 17)
3. Numerical Verification: Theory predictions verified on three representative bulk materials
4. Bridge Function: Establishment of connections between mathematical analysis, numerical schemes, and physical approximations

Limitations

1. Theoretical Aspects

• Nonlocal Potentials: Current theoretical framework does not yet include nonlocal effects of pseudopotentials (though numerical implementation uses them)
• Function Space Selection: While the choice of $H^{-1}_{per}$ and $H^1_{per}$ is reasonable, other choices might be superior

2. Numerical Aspects

• Inverse Crime: Forward and inversion use the same model; future work needs independent reference densities
• ε Sequence Optimization: Current simple exponential decay may not be optimal
• Stopping Criterion: Heuristic criterion of 0.01ε could be further optimized
• Computational Cost: Requires solving optimization problem for each ε value

3. Application Scope

• Currently limited to periodic insulator systems
• Validation on only three relatively simple materials
• System size limited to hundreds of electrons

Future Directions

Short-term Goals

1. Reference Density Sources: Application to densities from mean-field theory (beyond semi-local DFT)
2. Nonlocal Potential Theory: Refinement of theoretical framework to include nonlocal effects
3. Function Space Optimization: Exploration of other function space choices
4. Approximate Error Bounds: Development of more practical error estimates based on observation that $Q_\varepsilon$ may be constant

Long-term Goals

1. Functional Development: Utilization of rigorous inversion scheme to aid development of new approximate functionals
2. Hohenberg-Kohn Mapping: Deeper understanding of density-potential mapping
3. Quantum Embedding: Application to quantum embedding techniques
4. Optimized Effective Potential: Improvement of optimized effective potential methods
5. Complex Systems: Extension to larger, more complex material systems

In-Depth Evaluation

Strengths

1. Method Innovation (★★★★★)

• Theoretical Breakthrough: First successful application of MY regularization theory to actual physical systems, bridging theory and practice
• Mathematical Rigor: Provides unprecedented mathematical guarantees in KS inversion field
• Error Bounds: First establishment of computable, verifiable rigorous error bounds
• Non-expansivity Utilization: Clever use of non-expansivity from convex analysis to establish error theory

2. Experimental Sufficiency (★★★★☆)

• Multi-material Validation: Verification on three different material types (semiconductors, ionic crystals)
• Systematic Testing:
  • Exact inversion (noise-free)
  • Noisy inversion
  • Error bound verification
  • Convergence analysis
• Clear Visualization: Intuitive presentation through real-space potential plots and error diagrams
• Quantitative Analysis: Detailed numerical data and ratio analysis

Room for Improvement:

• Could include more material types (metals, strongly correlated systems)
• Quantitative comparison with other inversion methods
• Computational efficiency analysis

3. Result Convincingness (★★★★★)

• Theory-Experiment Consistency: Numerical results perfectly align with theoretical predictions
• Rigorous Error Bounds: All ratios remain within theoretical bounds
• Clear Convergence: Clear demonstration of convergence behavior as $\varepsilon \to 0$
• Robustness Verification: Proof of method stability against density perturbations

4. Writing Clarity (★★★★★)

• Logical Structure: Clear progression from theory → numerical implementation → experimental verification
• Mathematical Expression: Rigorous yet readable, with appropriate use of physical intuition
• Figure Quality: High-quality potential plots and error analysis diagrams
• Reproducibility: Complete open-source code (GitHub) and data (Zenodo) provided

Weaknesses

1. Method Limitations

• Computational Cost: Requires solving optimization problems for a series of ε values, potentially more expensive than traditional methods
• ε Selection: Lacks theoretical guidance for adaptive ε sequence selection
• Function Space Dependence: Results depend on specific function space choice; optimality not fully explored

2. Experimental Setup

• Inverse Crime: Acknowledged limitation requiring resolution in future work
• Material Diversity: Only three relatively simple materials tested
• Benchmark Comparison: Lacks direct quantitative comparison with other inversion methods (Wu-Yang, ZMP)

3. Theory Completeness

• Nonlocal Potentials: Theoretical framework not yet covering nonlocal pseudopotentials actually used
• Approximate Bounds: Current error bounds require computing $Q_\varepsilon$, potentially infeasible in practical applications
• Optimality: No proof that proposed method is optimal in any sense

Impact

1. Field Contribution (★★★★★)

• Paradigm Shift: Introduction of rigorous mathematical framework to KS inversion, potentially changing field research approaches
• Theoretical Foundation: Establishes solid foundation for developing more reliable inversion schemes
• New Error Analysis Path: Opens new direction for rigorous error estimation in density-potential inversion
• Interdisciplinary Bridge: Connects functional analysis, convex optimization, and quantum chemistry

2. Practical Value (★★★☆☆)

• Current Stage: Primarily proof-of-concept; direct practical utility limited
• Future Potential:
  • May improve functional development process
  • Provides tools for quantum embedding
  • Aids understanding of DFT fundamentals (band gaps, fractional charges)
• Computational Cost: Requires further optimization for routine computational use

3. Reproducibility (★★★★★)

• Open Source Code: Complete Julia implementation (based on DFTK)
• Public Data: Raw data publicly available on Zenodo (DOI: 10.5281/zenodo.14894064)
• Detailed Documentation: Comprehensive method and parameter descriptions
• Software Ecosystem: Built on mature DFTK platform, easily extensible

Applicable Scenarios

Ideal Application Scenarios

1. Functional Development: Approximate functional construction requiring rigorous error control
2. Benchmark Testing: Providing rigorous reference standards for other inversion methods
3. Theoretical Research: Exploration of DFT fundamentals (non-differentiability, v-representability)
4. Methodological Research: Development of new numerical inversion techniques

Current Limitation Scenarios

1. Large-scale Systems: Computational cost may limit applicability
2. Metallic Systems: Current implementation limited to insulators
3. Strongly Correlated Systems: Not tested on such systems
4. Real-time Applications: Unsuitable for scenarios requiring rapid inversion

Future Extension Scenarios

1. Quantum Embedding: As core component of embedding methods
2. Machine Learning: Providing high-quality training data for ML functionals
3. Uncertainty Quantification: Utilizing error bounds for uncertainty analysis
4. Multiscale Simulation: Information transfer between different accuracy levels

Key References

Theoretical Foundation

1. Penz, Csirik, Laestadius (2023): "Density-potential inversion from Moreau–Yosida regularization", Electron. Struct. 5, 014009 - Direct theoretical basis of this paper
2. Penz et al. (2019): "Guaranteed convergence of a regularized Kohn-Sham iteration in finite dimensions", Phys. Rev. Lett. 123, 037401
3. Laestadius et al. (2018): "Generalized Kohn–Sham iteration on Banach spaces", J. Chem. Phys. 149, 164103

DFT Fundamentals

4. Hohenberg & Kohn (1964): "Inhomogeneous electron gas", Phys. Rev. 136, B864
5. Kohn & Sham (1965): "Self-consistent equations including exchange and correlation effects", Phys. Rev. 140, A1133
6. Levy (1979): "Universal variational functionals of electron densities", Proc. Natl. Acad. Sci. USA 76, 6062

Inversion Methods

7. van Leeuwen & Baerends (1994): "Exchange-correlation potential with correct asymptotic behavior", Phys. Rev. A 49, 2421
8. Wu & Yang (2003): "A direct optimization method for calculating density functionals", J. Chem. Phys. 118, 2498
9. Shi & Wasserman (2021): "Inverse Kohn–Sham Density Functional Theory: Progress and challenges", J. Phys. Chem. Lett. 12, 5308

Numerical Tools

10. Herbst, Levitt, Cancès (2021): "DFTK: A Julian approach for simulating electrons in solids", Proceedings of the JuliaCon Conference 3, 69

Overall Evaluation

Recommended Audience: DFT theory researchers, quantum chemistry methodologists, computational materials scientists, researchers interested in numerical analysis and convex optimization.