Exchange Correlation Potentials from Full Configuration Interaction in a Slater Orbital Basis

TL;DR

This paper develops and tests the Ryabinkin-Kohut-Staroverov (RKS) method with Slater atomic orbitals, enabling accurate exchange-correlation potentials from full configuration interaction calculations, crucial for density functional theory.

Contribution

The work introduces the first implementation of the RKS method with Slater orbitals, demonstrating its efficiency and importance of cusp conditions for accurate potentials.

Findings

RKS method produces artifact-free exchange-correlation potentials.

Enforcing nuclear cusp conditions is essential for success.

Applicable to both weakly and strongly correlated systems.

Abstract

Ryabinkin-Kohut-Staroverov (RKS) theory builds a bridge between wave function theory and density functional theory by using quantities from the former to produce accurate exchange-correlation potentials needed by the latter. In this work, the RKS method is developed and tested alongside Slater atomic orbital basis functions for the first time. To evaluate this approach, Full Configuration Interaction computations in the Slater orbitals are employed to give quality input to RKS method, allowing full correlation to be present along with correct nuclei cusps and asymptotic decay of the wavefunction. The RKS method will be shown to be an efficient algorithm to arrive at exchange correlation potentials without unphysical artifacts in moderately-sized basis sets. Furthermore, enforcement of the nuclear cusp conditions will be shown to be vital for the success of the Slater-basis RKS method. Examples of weakly and strongly correlated molecular systems will demonstrate the main features of Slater RKS.

Tables

Table 1. Table 1: Normalized L 1 subscript 𝐿 1 L_{1} errors and L 2 subscript 𝐿 2 L_{2} errors of density from RKS with respect to CI densities in the QZ4P basis set

System	$\Delta {\rho }_{L_1}/N_e$	$\Delta {\rho }_{L_2}$
$H_2$(eq.)	$4.06\times {10}^{-3}$	$1.58\times {10}^{-3}$
$H_2$(2eq.)	$1.37\times {10}^{-2}$	$3.70\times {10}^{-3}$
$H_2$(d)	$8.86\times {10}^{-3}$	$1.49\times {10}^{-3}$
LiH	$4.66\times {10}^{-3}$	$4.68\times {10}^{-3}$
${\mathrm{CH}}_2$	$4.91\times {10}^{-3}$	$8.33\times {10}^{-3}$
$H_{2}O$	$3.28\times {10}^{-3}$	$1.88\times {10}^{-2}$

Equations

$$\bigg{[}-\frac{1}{2}\nabla^{2}+v(\mathbf{r})+v_{H}(\mathbf{r})+v_{XC}(\mathbf{r})\bigg{]}\phi_{i}(\mathbf{r})=\epsilon_{i}\phi_{i}(\mathbf{r}),$$

$$v_{XC}(\mathbf{r})=\frac{\delta E_{XC}[\rho (\mathbf{r})]}{\delta \rho (\mathbf{r})}.$$

$$\frac{\partial\phi_{i}}{\partial\mathbf{r}}\bigg{|}_{\mathbf{r}=R_{B}}=-Z_{B}\phi_{i}(R_{B}),$$

$$v_{XC}(\mathbf{r})=v_{XC,Slater}^{WF}(\mathbf{r})+\frac{{\tau }^{WF}(\mathbf{r})}{{\rho }^{WF}(\mathbf{r})}-\frac{{\tau }^{KS}(\mathbf{r})}{{\rho }^{KS}(\mathbf{r})}+{ϵ}^{KS}(\mathbf{r})-{ϵ}^{WF}(\mathbf{r}),$$

$$v_{XC,Slater}^{WF}(\mathbf{r})=\int \frac{{\rho }_{XC}^{WF}(\mathbf{r},{\mathbf{r}}_2)}{|\mathbf{r}-{\mathbf{r}}_2|}d{\mathbf{r}}_2,$$

$$\Gamma (\mathbf{r},{\mathbf{r}}_2;{\mathbf{r}}^{\prime },{\mathbf{r}}_2^{\prime }){|}_{\mathbf{r}={\mathbf{r}}^{\prime };{\mathbf{r}}_2={\mathbf{r}}_2^{\prime }}=\sum _{pqrs}{\Gamma }_{pqrs}{\psi }_p(\mathbf{r}){\psi }_q({\mathbf{r}}_2){\psi }_r^{\ast }({\mathbf{r}}^{\prime }){\psi }_s^{\ast }({\mathbf{r}}_2^{\prime }),$$

$${\tau }^{WF}(\mathbf{r})=\frac{1}{2}\sum _i{n}_i|\nabla {\psi }_i(\mathbf{r}){|}^2,$$

$${\tau }^{KS}(\mathbf{r})=\frac{1}{2}\sum _i{n}_i|\nabla {\varphi }_i(\mathbf{r}){|}^2,$$

$${ϵ}^{WF}(\mathbf{r})=\frac{1}{{\rho }^{WF}(\mathbf{r})}\sum _j{\lambda }_j|{\psi }_j(\mathbf{r}){|}^2,$$

$${ϵ}^{KS}(\mathbf{r})=\frac{2}{{\rho }^{KS}(\mathbf{r})}\sum _{i=1}^{N/2}{ϵ}_i|{\varphi }_i(\mathbf{r}){|}^2,$$

$$v_{XC}(\mathbf{r}{)}^{smooth}=F(\mathbf{r})v_{XC}(\mathbf{r})+(1-F(\mathbf{r}))v_{XC,Slater}^{WF}(\mathbf{r})$$

$$F(\mathbf{r})=\frac{\rho (\mathbf{r})}{\rho (\mathbf{r})+\theta }$$

$$\sum _{\alpha A}{p}_{\alpha A,B}{c}_{\alpha A,i}=0,$$

$$p_{\alpha A,B}=\frac{\partial\chi_{\alpha B}}{\partial\mathbf{r}_{B}}\bigg{|}_{\mathbf{r}_{B}=0}\delta_{AB}+Z_{B}\chi_{\alpha A}(\mathbf{r}_{B})$$

$$[(\mathbf{I}-\mathbf{A})(\mathbf{F}-ϵ\mathbf{S})(\mathbf{I}-\mathbf{A})]\mathbf{c}=0,$$

$$\Delta X_{L_{1}}=\int\big{|}X^{A}(\mathbf{r})-X^{B}(\mathbf{r})\big{|}d\mathbf{r}$$

$$\Delta X_{L_2}=\sqrt{\int (X^A(\mathbf{r})-X^B(\mathbf{r}){)}^{2}d\mathbf{r}}$$

$$(1a_1{)}^2(2a_1{)}^2(1b_2{)}^2(3a_1{)}^2$$

$$(1a_1{)}^2(2a_1{)}^2(1b_2{)}^2(1b_1{)}^2$$

Full text

I Introduction

Density Functional Theory (DFT)Mardirossian and Head-Gordon (2017); Becke (2014); Burke and Wagner (2013), particularly within the Kohn Sham (KS) formalism,Hohenberg and Kohn (1964) provides a cost-effective, scalable means for approximating the quantum behavior of electronic states. The KS equations are

[TABLE]

and $\phi_{i}$ are the KS orbitals, which together comprise a system of non-interacting electrons having the density $\rho^{KS}(\mathbf{r})=\sum_{i}n_{i}|\phi_{i}(\mathbf{r})|^{2}$.Kohn and Sham (1965) The orbital energies, $\epsilon_{i}$, arise from the kinetic energy and three potentials: $v(\mathbf{r})$ is the external potential (usually the nuclear potential), $v_{H}(\mathbf{r})$ is the Hartree potential which accounts for the classical Coulomb interaction between electrons, and $v_{XC}(\mathbf{r})$ is the exchange-correlation potential which contains the non-classical contribution from the kinetic energy and electron-electron repulsion. The unknown component of DFT, $v_{XC}(\mathbf{r})$, is defined as the functional derivative of the XC energy functional,

[TABLE]

The functional of the density, $E_{XC}[\rho(\mathbf{r})]$, and its derivative in Eq. 2 are approximated in practice. Together, Eqs. 1 and 2 are used to solve for KS orbitals, which provide a description of the electron density and its corresponding energy.

The choice of the approximate form for $E_{XC}[\rho(\mathbf{r})]$ is made from the various rungs of the ’Jacob’s Ladder’ of functionals that have been developed over the past half century.Perdew et al. (2005) The capabilities of DFT for molecules Cohen and Handy (2000); Kümmel and Kronik (2008); Sun et al. (2016) and solidsDiStasio et al. (2014); Marsman et al. (2008); Isaacs and Wolverton (2018) have undoubtedly increased over this time, but questions still remain regarding well-documented problems with self-interaction, delocalization errors, and strong correlation. Cohen, Mori-Sánchez, and Yang (2012); Verma and Truhlar (2020); Crisostomo et al. (2022); Bryenton et al. (2022) Since these errors are tied closely to the electron density, it is clear that improvements to $v_{XC}(\mathbf{r})$ are needed to improve upon current DFT approximations.

The focus of this article is therefore on the exchange-correlation potential, since knowledge about this potential could lead to improved functional approximations. Green, Tozer, and Handy (1998); Tozer and Handy (1998); Wilson, Bradley, and Tozer (2001); Menconi, Wilson, and Tozer (2001); Gaiduk, Mizzi, and Staroverov (2012) This point has been emphasized by recent work showing that widely used functionals give rise to potentials that are far from the exact (even for SCAN0 which has $v_{XC}$ closest to the exact, the errors are in the range of $O(10^{-1}-10^{0})$).Kanungo, Zimmerman, and Gavini (2021) These errors may be attributed to the standard practice of functional training to reproduce energies, while neglecting errors in the exchange-correlation potential. Providing this information, however, requires means to accurately compute $v_{XC}(\mathbf{r})$, which is a highly nontrivial task.

For a given reference density, the corresponding XC potential can in principle be obtained numerically by the inverse DFT approach.Shi and Wasserman (2021) The inverse problem maps a given density to its corresponding potential through a unique one-to-one mapping.Hohenberg and Kohn (1964); Runge and Gross (1984) The inverse relation, however, is well posed only in a complete basis.Kohn (1983); Hadamard (1902) In an incomplete basis the unique mapping of potential to density does not hold true, i.e., different XC potentials can map to the same density.Harriman (1986, 1990); Staroverov, Scuseria, and Davidson (2006) In practice, inverse calculations are often performed using incomplete basis sets (such as finite Gaussian basis sets) and persistent numerical problems result. Recently, Kanungo et al. implemented a complete, finite-element basis and used it to obtain highly accurate potentials,Kanungo, Zimmerman, and Gavini (2019) and Stückrath et al. alleviated the problem using a multiresolution wavelet basis to similar effect.Stückrath and Bischoff (2021) These recent results show the possibilities of complete basis sets in tackling the inverse DFT problem, but have not improved upon the situation for finite basis sets.

Elegantly skipping past the inverse DFT problem, Staroverov and coworkers developed means to obtain $v_{XC}(\mathbf{r})$ from the reduced density matrices (RDMs) of any wavefunction method.Ryabinkin, Kohut, and Staroverov (2015); Cuevas-Saavedra, Ayers, and Staroverov (2015) The RKS method evaluates $v_{XC}(\mathbf{r})$ from the comparison of two local energy balance equations, one originating from the KS equations and the other from wavefunction theory. Since the RKS method uses the wavefunction instead of the density to obtain the potential and employs systematic approximations to the kinetic energy density, it is less prone to finite-basis-set errors as compared to related techniques. Kumar, Singh, and Harbola (2020); Shi, Chávez, and Wasserman (2022) Using RKS, potentials can be straightforwardly derived for atoms and molecules in moderate or larger Gaussian basis sets. Even though densities expanded in Gaussian type orbitals (GTOs) have inherent shortcomings, the stability of the RKS (and mRKS) methods in obtaining accurate exchange-correlation potentials is a promising step forward for finite-basis DFT. To go further with the RKS method, we envisioned improving the wave function, electron densities, and $v_{XC}(\mathbf{r})$ by the use of Slater type orbitals (STO). Schipper, Gritsenko, and Baerends (1997)

STOs have good physical properties that are expected for electronic states, for instance the ability to describe cusps at the nuclei and at long range, exponentially decaying tails. Kato (1957); Helgaker, Jørgensen, and Olsen (2000); Reinhardt and Hoggan (2009) Particularly in the construction of densities and the corresponding potentials, these properties make STOs a better choice than GTOs, which suffer from unphysical oscillatory behaviour under the Laplacian.Schipper, Gritsenko, and Baerends (1997) While STOs have seen less use than GTOs in quantum chemistry, their limitations can be traced to the difficult two-electron integral evaluation, which must be done numerically in 6 dimensions. Recently Zimmerman and coworkers developed an efficient GPU-accelerated algorithm to evaluate STO integrals in the resolution-of-the-identity (RI) approximation.Dang, Wilson, and Zimmerman (2022) This opened up the opportunity for us to carry out the RKS procedure with STO basis sets.

By using STO basis sets, one important property of the electron density can be naturally incorporated into the RKS procedure. The Kato cusp condition specifies how the density must behave near the nucleus.Kato (1957) This condition is met when every occupied molecular orbital satisfies

[TABLE]

where $R_{B}$ is the position of a nucleus $B$ with atomic number $Z_{B}$. With STOs, the cusp condition can be enforced by using a modified SCF procedureHandy (2004) which will be described in the Methods section. As a result, the singularity due to electron-nuclear attraction at the nucleus will be compensated by the kinetic energy. This property appears unreachable in finite GTO basis sets, and its effect on $v_{XC}(\mathbf{r})$ near the nucleus will be an interesting question to examine (Section IV).

This article presents exchange correlation potentials for atoms and molecules, derived from highly accurate full configuration interaction (FCI) densities using Slater type orbitals. FCI wave functions in STO basis sets capture all dynamic and static correlation available to the basis, with the additional benefit of having the correct nuclear cusp and long range asymptotic behavior. The ability to calculate XC potentials at this high level of theory will therefore be examined for the first time in the results that follow.

II Method

The RKS method utilizes the one particle and two particle reduced density matrices (1-RDM and 2-RDM respectively) of the wavefunction to evaluate $v_{XC}$.Ryabinkin, Kohut, and Staroverov (2015) The notations $WF$ or $KS$ are used throughout this manuscript to denote the terms derived from wavefunction or KS-DFT methods, respectively. The RKS working equation,

[TABLE]

is derived from local energy balance equations under the condition $\rho^{WF}(\mathbf{r})=\rho^{KS}(\mathbf{r})$. Here, $\rho^{WF}(\mathbf{r})$ is fixed to that of the wavefunction reference, and $\rho^{KS}(\mathbf{r})$ is the density from the self-consistent solutions of the KS equations (Eq. 1). The Slater exchange-correlation charge potential, $v_{XC,Slater}^{WF}(\mathbf{r})$,Slater (1953) is expressed as,

[TABLE]

where, $\rho^{WF}_{XC}(\mathbf{r})$ is the exchange-correlation hole density derived from the relation, \Gamma(\mathbf{r},\mathbf{r}_{2};\mathbf{r},\mathbf{r}_{2})=P(\mathbf{r},\mathbf{r}_{2})=\frac{1}{2}\rho^{WF}(\mathbf{r})\big{[}\rho^{WF}(\mathbf{r}_{2})+\rho^{WF}_{XC}(\mathbf{r},\mathbf{r}_{2})\big{]}. Here $\Gamma(\mathbf{r},\mathbf{r}_{2};\mathbf{r},\mathbf{r}_{2})$ is the special case of the coordinate representation of the 2-RDM,

[TABLE]

where $\Gamma_{pqrs}$ are the matrix elements of the orbital representation of the 2-RDM. The positive-definite kinetic energy densities, $\tau(\mathbf{r})$, and the average local electron energies, $\epsilon(\mathbf{r})$, are defined as,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the orbitals $\psi_{j}$ are associated with generalized Fock eigenvalues $\lambda_{j}$, and $\phi_{i}$ are the set of KS orbitals with eigenvalues $\epsilon_{i}$. Other choices for the kinetic energy densities are the Laplacian and Pauli forms. In the original RKS method, the positive-definite kinetic energy was used rather than the Laplacian kinetic energy ($\tau_{L}^{KS}=-\frac{1}{2}\sum_{i}n_{i}\phi_{i}^{*}(\mathbf{r})\nabla^{2}\phi_{i}(\mathbf{r})$) which are related by $\tau_{L}=\tau-0.25\nabla^{2}\rho$.Ryabinkin, Kohut, and Staroverov (2015); Cuevas-Saavedra, Ayers, and Staroverov (2015) This choice can be motivated by equating KS and WF densities (i.e. $\rho^{KS}=\rho^{WF}$) or by observing that the positive-definite kinetic energy is more numerically stable near an atom. In the modified RKS method, the Pauli kinetic energy density ($\tau_{P}=\tau-|\nabla\rho|^{2}/8\rho$)Ospadov, Ryabinkin, and Staroverov (2017) is used to proceed further along this path. For Slater-based RKS, two forms of kinetic energy densities were evaluated (Laplacian and positive definite), and it was found that the positive definite form has better properties (vide infra).

The above RKS equations give $v_{XC}(\mathbf{r})$ up to an arbitrary constant. At the asymptotic limit of $\mathbf{r}\to\infty$, $v_{XC}(\mathbf{r})\to v_{XC,Slater}^{WF}(\mathbf{r})\sim-1/\mathbf{r}$. $\epsilon^{KS}$ and $-\tau^{KS}/\rho^{KS}$ approach $\epsilon_{HOMO}$, whereas the analogous wavefunction terms approach $-I_{min}$, the first ionization energy from extended Koopman’s theorem. In order to enforce the asymptotic decay of $v_{XC}(\mathbf{r})$ (and determine the constant), all the terms to the right of Eq. 4 except for the $v_{XC,Slater}^{WF}$ must cancel at $\mathbf{r}\to\infty$. Therefore all the eigenvalues of the KS orbitals are shifted such that $\epsilon_{HOMO}=-I_{min}$ is satisfied. Further, in the far field, where the densities are almost zero, the numerical evaluation of the kinetic energy densities, i.e. the second and third terms in Eq. 4 yields unphysically large values. Since at large distances $v_{XC}$ and $v_{XC,Slater}^{WF}$ are expected to be identical, a function ($F$) is used to smoothly transition $v_{XC}$ into $v_{XC,Slater}^{WF}$ at the threshold of low density ($\theta=10^{-5}$),

[TABLE]

where,

[TABLE]

The algorithm of the RKS method (Fig. 1) entails evaluating the wavefunction terms first, i.e., $v_{XC,Slater}^{WF}(\mathbf{r})$, $\rho^{WF}(\mathbf{r})$, $\tau^{WF}(\mathbf{r})$, and $\epsilon^{WF}(\mathbf{r})$. Next, an initial guess for the KS orbitals is chosen and the corresponding terms, $\rho^{KS}(\mathbf{r})$, $\tau^{KS}(\mathbf{r})$ and $\epsilon^{KS}(\mathbf{r})$ are evaluated and the $\epsilon_{i}$ are shifted with respect to $I_{min}$. Eq. 4 then provides $v_{XC}(\mathbf{r})$. This potential is then used to solve for new KS orbitals from Eq. 1, and the KS terms are updated in Eq. 4 until $v_{XC}(\mathbf{r})$ and the KS orbitals become self-consistent.

In all, the SlaterRKS procedure follows the original RKS procedure that used a Gaussian basis set, except for a few details. In particular, it is notable that due to the correct short- and long-range behavior of Slater orbitals, the Laplacian kinetic energy operator, $\tau^{KS}_{L}(\mathbf{r})$ may be a tractable choice. This will be the case only when the conditions of the next paragraph are under consideration.

The molecular orbitals used for the WF and KS theories are constrained to obey Kato’s cusp condition (Eq. 3) using Handy’s method.Handy (2004) Each molecular orbital is expanded in terms of the basis functions ($\chi$), i.e. $\phi_{i}=\sum_{\alpha A}c_{\alpha A,i}\chi_{\alpha A}$, where the sum is over exponent $\alpha$ and the atom $A$ on which the basis functions are centered, and $c_{i}$ are the MO coefficients corresponding to the $i^{th}$ MO. Kato’s condition gives

[TABLE]

where,

[TABLE]

The condition of Eq. 13 is enforced by using a modified SCF procedure with,

[TABLE]

where $\mathbf{A}=\sum_{BC}\hat{\mathbf{p}}_{B}(\hat{\mathbf{p}}_{B}^{T}\hat{\mathbf{p}}_{C})^{-1}\hat{\mathbf{p}}_{C}^{T}$. This procedure ensures that all MOs in the correlated wavefunction have correct electron-nuclear cusps, and the same will be true when solving the KS equations.

III Computational Details

The SlaterRKS method was implemented in C++ and interfaced with the SlaterGPUDang, Wilson, and Zimmerman (2022) library to evaluate integrals for Slater basis functions using GPU acceleration. All Coulomb integrals will therefore be evaluated using the Resolution of the Identity (RI) approximation.Dunlap, Connolly, and Sabin (1979); Werner, Manby, and Knowles (2003); Distasio Jr. et al. (2007) The orbitals and RDMs for the reference WFs come from the heat-bath configuration interaction (HBCI) procedureHolmes, Tubman, and Umrigar (2016); Sharma et al. (2017); Li et al. (2018); Dang, Kammeraad, and Zimmerman (2023); Chien et al. (2018) using a tight threshold (He, $\mathrm{H_{2}}$, LiH: $10^{-5}$ Ha and $\mathrm{H_{2}O}$, $\mathrm{CH_{2}}$: $10^{-4}$ Ha) for accuracy of the variational wavefunction. For atoms and small molecules, the RDMs from this approach will be essentially FCI quality.

The primary working equation in the SlaterRKS method is Eq. 4, where all the quantities are evaluated according to their expressions given in Eqs 5-10. Evaluation of these equations is accelerated on GPU using OpenACC to generate their contribution to $v_{XC}$ on the grid. For Eqs 7-10, a simple acc parallel directive is prepended to each for-loop over the grid points. Evaluation of $v_{XC,Slater}^{WF}(\mathbf{r})$ via Eq. 5 is the most expensive step due to the 6-dimensional integration over two electrons. To minimize wall-time, OpenMP and OpenACC are jointly used to compute $v_{XC,Slater}^{WF}(\mathbf{r})$. Therefore, contributions from Eq. 5 involve two nested for-loops, each over the entire grid. The outer for-loop is parallelized with OpenMP where one thread is launched for each GPU. The inner for-loop is then parallelized with OpenACC using the acc parallel directive with a reduction clause. The reduction clause indicates a summation in the inner for-loop, which reflects the numerical evaluation of the integral in Eq. 5. Since the same grid and integral weights are used for each inner for-loop, they only need to be generated once before entering the outer for-loop of Eq. 5.

The Slater basis sets were taken from the set developed by Baerends and coworkers.Van Lenthe and Baerends (2003) Four types of basis set are included in this set, namely DZP, TZP, TZ2P and QZ4P. The DZP, TZP and TZ2P basis sets have double zeta-core functions, while the QZ4P has a triple-zeta core. In the valence region, the DZ is double zeta, TZP and TZ2P triple zeta, and QZ4P quadruple zeta, all with valence polarization functions. For example, the QZ4P basis has 3$\times$1s, 4$\times$2s, 4$\times$2p, 2$\times$3d and 2$\times$4f functions for C atom. Larger basis sets were needed to test the cusp condition implementation on the helium atom (vide infra). 5Z6P and 6Z6P basis sets were thus created in an even-tempered manner, following the procedure of Baerends and coworkersChong et al. (2004) (see Supporting Information for full specification of these basis sets).

The geometries and ionization energies ($I_{min}$) are taken from the CCCBDB databaseJohnson (1999) and provided in the Supporting Information. The value of $I_{min}$ can be set to the ionization energy from simulation or to the experimental value. Herein, $I_{min}$ is set to the experimental value (given in the SI), as the FCI result is expected to be similar. The STO integral evaluation as well as the RKS procedure is carried out on numerical atom-centered grid.Becke (1988); Mura and Knowles (1996); Murray, Handy, and Laming (1993) The three dimensional grid is composed of products of radialMura and Knowles (1996) and angularLebedev (1976) points with weights according to the Becke partitioning scheme.Becke (1988) In all the calculations in this manuscript, 50 radial and 5810 angular points are employed for each atom.

The convergence of each SlaterRKS run is verified by the $L_{1}$ and $L_{2}$ errors in the KS density compared to the reference WF. The norms, $\Delta X_{L_{1}}$ and $\Delta X_{L_{2}}$ of a property $X$ are defined as,

[TABLE]

[TABLE]

where $A$ and $B$ are the reference and calculated values. The calculation is deemed converged when the $\Delta\rho_{L_{1}}$ does not change more than $10^{-5}$ from one iteration to the next.

IV Results and Discussion

In this section the SlaterRKS method is applied to a few prototypical test systems along with a few strongly correlated test cases. The simplest example is the two electron case of hydrogen molecule at three separate bond distances–equilibrium, twice the equilibrium and fully dissociated. The next example is the heteronuclear LiH molecule. In addition the water molecule is examined as a small polyatomic, followed by the more challenging multireference test case of singlet methylene ($\mathrm{CH_{2}}$). For all of these examples, FCI wavefunctions in the Slater orbital basis sets will be used (unless otherwise noted) as the reference for the RKS procedure.

Before delineating these examples, it is emphasized that the enforcement of the long-range asymptotic behavior is important to reach meaningful exchange correlation potentials. The correct asymptotic decay of $v_{XC}$ follows $-1/\mathbf{r}$,Schmidt et al. (2014); Wu, Ayers, and Yang (2003), though LDA and GGA functionals fail to satisfy this condition.Staroverov et al. (2004) Properties such as the energy of the HOMO, which is tied closely to the ionization energy, will only be accurate with correct asymptotics.Kraisler (2020) In the RKS method, this condition is enforced by shifting $v_{XC}$ such that $\epsilon_{HOMO}=-I_{min}$. In SlaterRKS, an additional condition (see method section) ensures a smooth transition $v_{XC}(\mathbf{r})\to v_{XC,Slater}^{WF}(\mathbf{r})$ as $r\to\infty$. The Supporting Information shows the $\mathrm{H_{2}}$(eq.) molecule, where $v_{XC}$ and $v_{XC,Slater}^{WF}$ decay as $-1/\mathbf{r}$ at low density regions far from the nuclei.

While the cusp condition (Eq. 3) will be applied for most of the SlaterRKS results of this work, it will first be evaluated for the helium atom in a range of basis sets. Numerical evaluation of Kato’s equation on a few grid points near the nucleus is shown in Tables LABEL:tab:cusptz-LABEL:tab:cusp6z in the SI. These results confirm the cusp enforcement reduces the error at $\mathbf{r}=10^{-6}$ (i.e. at the nuclear position) by two orders of magnitude or more compared to densities without the cusp condition. Obtaining this accuracy near the nucleus, however, is not free because one degree of freedom per atom is lost in the basis when the constraint is applied. Therefore the dependence of the SlaterRKS results on the basis set size needs to be examined before applying the cusp condition more widely.

Fig. 2 shows the $L_{1}$ error of density from the SlaterRKS procedure (blue) and the corresponding FCI correlation energy ($E_{corr}$) (red) with cusp condition (solid line) and without (dashed line) for the helium atom. With increase in basis set size, the correlation energy improves systematically and converges.not In all basis sets, the $E_{corr}$ is lower when the cusp condition is enforced. In the smaller TZ2P basis set this effect is quite significant, resulting in a 5 mHa decrease in correlation energy at the FCI level. The effect diminishes to the sub-mHa level with larger basis sets (QZP and higher). The errors in the KS density (compared to the WF reference) behave similarly. With larger basis sets, the $L_{1}$ errors for the pairs of densities agree within $O(10^{-4})$ a.u. and converge to a common point (see Table LABEL:tab:cuspbasis in the supporting information). Differences remain in $L_{1}$ errors between the cusp-enforced vs. cusp-not-enforced densities at the TZ2P level, reflecting the smallness of this basis. For the larger basis set sizes, the cusp condition can readily be applied, giving an improved description of the density at the nucleus.

Having established how the SlaterRKS procedure behaves with and without a cusp condition, the next item to examine is the effect of basis set size, this time in a diatomic molecule. Results for the $\mathrm{H_{2}}$ molecule (at its equilibrium geometry) with various basis sets are therefore shown in Fig. 3. The DZP basis set is too small to be meaningfully used alongside the cusp condition, but the others include the Kato cusp. Far from the nuclei where the density is low, the three different basis sets result in similar potentials. As more correlation is present near the bonding region, the larger basis sets (TZ, QZ) produce deeper wells in the exchange correlation potential, though the three basis sets give qualitatively similar structures consisting of double-well potentials. For comparison, the SlaterRKS potential from the HF wavefunction is also shown in Fig. 3. Notably, $v_{XC}$ from HF is missing the characteristic maximum in the middle of the bond, due to the lack of electron correlation.

Fig. 4 shows the differences in density between FCI and SlaterRKS for the QZ4P, TZ2P and DZP basis sets. The corresponding $L_{1}$ and $L_{2}$ errors in the density for these calculations are given in the Supporting Information. As reflected in the original RKS procedure in a Gaussian basis, larger basis sets do a better job at reproducing WF electron densities.Cuevas-Saavedra, Ayers, and Staroverov (2015); Ryabinkin, Kohut, and Staroverov (2015) The Laplacian kinetic energy ($\tau_{L}$) can also be used in place of the positive-definite form, which is shown for the QZ4P basis set (see Supporting Information). The sharp features in the $\tau_{L}-v_{XC}$ are due to the divergence of the Laplacian near the nucleus, which arise due to the differences in WF and KS KE density via Eq. 4. In all, because the larger (QZ4P) Slater basis set and positive-definite KE operator produced the best density (as well as potential), the remaining SlaterRKS results in this work were carried out using the QZ4P basis set, FCI RDM, and the positive-definite KE operator.

To show the behavior of SlaterRKS as strong correlation is introduced to the reference WF, three geometries of the $\mathrm{H_{2}}$ molecule were examined. The $L_{1}$ errors of densities normalized by the number of electrons ($\Delta\rho_{L_{1}}/N_{e}$) are given in Table 1 for the QZ4P basis set. Excellent agreement of the RKS and FCI densities in the case of the equilibrium geometry of $\mathrm{H_{2}}$ is demonstrated by the normalized $L_{1}$ error in the order of $10^{-3}$. The increase in the $L_{1}$ errors at the elongated bond lengths is expected as strong correlation increases since the FCI density has significantly noninteger orbital occupancies, which are challenging to reproduce with a single determinant in a finite basis set. The features of the $v_{XC}$s (shown along the H–H bond in Fig. 5) are largely representative of the exact potentials obtained using the finite-element inverse calculation.Kanungo, Zimmerman, and Gavini (2021) For the strongly correlated cases of stretched $\mathrm{H_{2}}$ molecules, where most functionals fail, the SlaterRKS potentials are promising. The depth of the potential at the nucleus and the height of the maxima between the nuclei gradually increase as the bond length increases, signifying gradual depletion of electron density between the bond.

Moving on to a diatomic with more electrons, SlaterRKS analysis for the minimum energy geometry of LiH is shown in Fig. 6 (top). Various commonly used DFT XC potentials have a deeper well at Li, whereas the SlaterRKS $v_{XC}$ is shallower, closer to the accurate finite-element inverse calculation.Kanungo, Zimmerman, and Gavini (2021) A comparison with $v_{XC}$ profiles obtained from FCI, CASCI with a (4e,4o) active space and Hartree Fock (HF) reveals the role of fully correlated wavefunctions in RKS. The FCI has deeper wells at Li and H whereas both CASCI and HF have shallower wells and are almost indistinguishable (Fig. 6). To the left of the Li potential well, an intershell feature is present, typical of an accurate $v_{XC}$.Kanungo, Zimmerman, and Gavini (2021) The intershell feature gradually shifts away from Li in the sequence FCI to CAS to HF. The errors in the RKS densities constructed from various wavefunctions from the FCI density in the QZ4P basis set are shown in Fig. 6 (bottom). The FCI RKS density has largest errors near the Li nucleus. This is likely the case because the QZ4P basis has only three core 1s orbitals, limiting the ability of SlaterRKS to resolve the density to higher accuracy. On the other hand, the CAS and HF RKS densities have larger errors at both Li and H due to lack of correlation. The corresponding $\Delta\rho_{L_{1}}/N_{e}$ of the CAS and HF RKS densities with respect to the FCI density are $1.32\times 10^{-2}$ and $1.39\times 10^{-2}$, respectively. Overall, the $L_{1}$ errors in the RKS density from the FCI reference are still relatively low, as reflected in Table 1. These might be improved further with the availability of larger Slater basis sets.

The polyatomic $\mathrm{H_{2}O}$ molecule was next subjected to SlaterRKS analysis. The difference between QZ4P-FCI and TZ2P-FCI $v_{XC}$ is shown in Fig. 7 in the molecular plane. The largest differences are near the O nuclei as the QZ4P basis set has a deeper well near $\mathbf{r}=0$. This difference gradually decreases as distance from the O nucleus increases until about 0.3-0.4 Bohr, after which a diffuse yellow band depicts the difference in the intershell feature in the two basis sets. The bonding region however is quite similar in both QZ and TZ with slightly higher differences appearing near the H nuclei.

The importance of having correct nuclear cusps was emphasized in the introduction, and an explicit enforcement of these cusps is present in the SlaterRKS algorithm. The $\mathrm{H_{2}O}$ molecule provides a good example to demonstrate the nature of SlaterRKS $v_{XC}$ with and without the cusp condition being enforced (Eq. 3). Fig. 8 shows $v_{XC}$ in a small region around the oxygen nucleus. The peak that forms without the cusp condition is unexpected, and likely an artifact of using Slater orbitals that have steep slopes near the nucleus (see Eq. 8). Since Gaussian orbitals have zero slope at the nucleus, such behavior is less likely with finite-sized Gaussian basis sets. Fortunately, enforcement of the nuclear cusp condition dramatically remedies this situation, producing a more physical, single well potential near the oxygen nucleus. The middle and right side of Fig. 8 show that the better behavior in the cusp-enforced case stems from smoother contributions to $v_{XC}$ from $\epsilon$ and $\tau$, since the Slater potentials are monotonic (see SI). The contributions from $\epsilon$ and $\tau$ without the cusp condition have more features—which do not cancel out—and overall result in a significant, likely incorrect effect on $v_{XC}$.

Finally, SlaterRKS is used to examine a strongly correlated polyatomic. The lowest energy singlet state (${}^{1}A_{1}$) of the $\mathrm{CH_{2}}$ (methylene) molecule has significant contributions from two electron configurations,

[TABLE]

and,

[TABLE]

giving this state a multiconfigurational nature.Sherrill et al. (1998); Zimmerman et al. (2009); Chien et al. (2018) Fig. 9 shows the SlaterRKS XC potential along one of the C–H bonds. Using FCI the depth of the potential at the carbon nucelus is -4.7 a.u. which is in good agreement to previous reports.Morrison and Zhao (1995); Schipper, Gritsenko, and Baerends (1998) Notably, the potential at both the nuclei is shallower when CASCI with a (2e,2o) active space is used. The intershell structure distinguishing the core and valence regions around the C atom is also present, as expected. This feature, however, is quantitatively different when the two methods are compared, with deviations on the order of 0.1 a.u. The potential at the H nucleus resembles the trailing valence region of the C atom. The nature of the $v_{XC}$ in the bonding region of the C–H bond of $\mathrm{CH_{2}}$ is relatively flat, which is a characteristic feature of covalent bonds.Gritsenko, Leeuwen, and Baerends (1996)

The potential from Fig. 9 results in an $L_{1}$ error in the electron density of $O(10^{-3})$, which is a satisfactory result given the complexity of the reference WF, which is from complete CI computations in the QZ4P basis set. The CI result includes population of 0.05958 electrons in the p orbital, a significant amount that is difficult to capture in the pure-state KS representation of the density. While this example has not been studied before using RKS theory, the SlaterRKS method appears up to the task of treating strong correlation in the challenging $\mathrm{CH_{2}}$ molecule.

V Conclusions

This work found that RKS theory using a Slater basis set is a useful tool for examining exchange-correlation potentials corresponding to FCI wavefunctions. Features of the potentials match well with complete-basis-set results, compared to the finite-element inverse DFT results of Kanungo et al.Kanungo, Zimmerman, and Gavini (2019, 2021) At the same time, however, quantitative accuracy of reproducing the FCI densities using moderately-sized Slater basis sets was also moderate, with $L_{1}$ errors on the order of $10^{-3}$ per electron. At the time of this work, Slater atomic orbital basis sets do not extend beyond quadruple zeta in quality, hindering progress towards more accurate densities via SlaterRKS. Future work to build larger, more complete Slater basis sets is likely to be instrumental in improving the accuracy of this method.

We anticipate that the SlaterRKS approach will be useful in providing facile analysis of strongly correlated molecules and their exchange-correlation potentials, using highly accurate reference wavefunctions (FCI) expressed in finite basis sets. The relative ease in convergence, physical behavior of the potentials, and correct asymptotics in the density are key advantages of SlaterRKS that merit further consideration of this method.

VI Supplementary Material

See the supplementary material for the geometries of the molecules, ionization energies, details of integration grid, description of convergence of cusp correction and the even-tempered basis set used for He, asymptotic behaviour of exchange correlation potential, basis set dependence of SlaterRKS method and effect of Laplacian kinetic energy densities.

VII Acknowledgements

This project has been supported by the Department of Energy through the grant DE-SC0022241. The authors acknowledge the computing time on the Perlmutter Supercomputer from the National Energy Research Scientific Computing Center (NERSC) through allocation m4067. ST thanks support by the Eric and Wendy Schmidt AI in Science Postdoctoral Fellowship, a Schmidt Futures program.

VIII References