Direct determination of optimal real-space orbitals for correlated electronic structure of molecules

We demonstrate how to determine nearly numerically exact orthonormal orbitals that are optimal for evaluation of the energy of arbitrary (correlated) states of atoms and molecules by minimization of the energy Lagrangian. Orbitals are expressed in real space using multiresolution spectral element basis that is refined adaptively to achieve the user-specified target precision while avoiding the ill-conditioning issues that plague AO basis set expansions traditionally used for correlated models of molecular electronic structure. For light atoms the orbital solver in conjunction with a variational electronic structure model [selected Configuration Interaction (CI)] provides energies of comparable precision to a state-of-the-art atomic CI solver. The computed electronic energies of atoms and molecules are significantly more accurate than the counterparts obtained with the Gaussian AO bases of the same rank, and can be determined even when linear dependence issues preclude the use of the AO bases.


Introduction
Predictive computation of properties of molecules and materials requires highly-accurate treatment of many-electron correlation effects as well as precise numerical representation of the many-body states and operators involved therein. In molecular and, increasingly, in solid-state contexts many-body methods have utilized the Linear Combination of Atomic Orbitals (LCAO) numerical representation of 1-particle (orbital) states constructed from pre-optimized sequences of atom-centered basis sets (represented by analytic Gaussian-type or Slater-type functions, or represented numerically 1 ). However, (augmented) plane-wave (PW) representation of 1-particle states can also be used. 2 Unfortunately representation of many-particle by products of 1-particle states suffers from slow convergence in the vicinity of electron cusps, 3 thereby mandating the use of basis set extrapolation 4 or the augmentation with explicit 2-particle basis ("explicit correlation"). [5][6][7][8] Although the LCAO-based many-body methods can be pushed match experimental thermochemical 9 and spectroscopic data 10 for small molecules, extending such successes to increasingly larger systems, and to higher accuracy, is limited by the rapid onset of illconditioning, high-order rise of operator evaluation with the AO angular momenta, and the challenges of fast (reduced-scaling) reformulation of such methods due to the increasing density of the representation in the high-precision regime. LCAO representation in practice is a heavily empirical endeavour: at this moment the Basis Set Exchange reference database 11 lists 402 Gaussian AO basis sets designed for representing orbitals of a carbon atom. Albeit the large number of officially-recognized bases is simply due to the long history of the Gaussian AO technology, many recently-designed bases have rather niche uses; navigating the basis set zoo is a near-expert task, especially when treating high-density solids 12 or computing properties other than the energy. 13 Meanwhile plane wave representation struggles with description of localized features (e.g., atomic shell structure, nuclear cusps and Coulomb holes) of the electronic wave functions and thus are not naturally suitable to compact and/or fast formulations.
An attractive alternative to the LCAO and PW representations for many-body electronic structure are real-space numerical representations, a loosely-defined group including finite difference (FD), which only use grids and lack explicit basis functions, and finite element (FE) and spectral element (SE) methods, which utilize basis functions with strictly local support. Such representations combine the ability to resolve short and long lengthscale features of the many-body electronic wave functions with promise of systematic bias-free improvability without ill-conditioning issues. Due to these favorable features there has been a flurry of activity in developing real-space numerical technology for fast 1-body methods, such as Kohn-Sham density functional theory.
Unfortunately the use of real-space representation has been far more limited for manybody methods applicable to general molecules due to the exponential growth of the representation of a k-particle state with k. Even methods like MP2 and CCSD in real-space formulation involve representation of 2-particle wave operators (hence, 6-dimensional functions); while for atoms pair theories have been explored for a relatively long time, [14][15][16][17] for molecules the feasibility of such representations has been demonstrated only recently by some of us (MP2) 18,19 and our collaborators (CIS(D), CC2) 20,21 by combining volume-element-wise pair function compression and R12/F12-style regularization of electron cusps. 18 It is possible to obtain energies and properties with such methods that are competitive to the state-of-theart LCAO technology for small polyatomics, but the methodology is still relatively limited and likely cannot be made practical 3-and higher-body methods.
A more practical alternative to the direct real-space representation of 2-particle states is to expand them in products of real-space orbitals. Although such approaches, unlike the direct many-body representations, suffer from the slow convergence near electron cusps, the same remedies as used for the LCAO-based many-body methods (namely, the basis set extrapolation or R12/F12-style explicit correlation) can be deployed to fix these shortcomings. In the purely orbital-based many-body framework the key differences between spectral (LCAO, PW) and real-space representations is limited to the orbital optimization and evaluation of the physical operators. Due to the cost of both in the real-space representation real-space many-body methods have been implemented for atoms (both at the multiconfigu- The first demonstration of any real-space orbital-based many-body method applicable to general molecules was shown recently by Kottmann, Bischoff, and Valeev. 32 The key to this advance was the realization that far more compact orbital representation of the weaklyoccupied (virtual) orbitals needed for dynamical correlation is provided by the pair-natural orbitals. [33][34][35] The PNO-MP2-F12 method using the real-space representation for all orbitals was demonstrated for molecules with more than 30 atoms, with precision surpassing what is achievable by the state-of-the-art LCAO technology.
While the work showed a viable roadmap for developing many-body perturbation theory (MBPT) and related methods with real-space orbitals, a single Slater determinant often does not provide a qualitatively correct starting point. In this paper we ask whether it is possible 1. to express arbitrary (not just MBPT-based) orbital-based many-body wave functions in terms of real-space orbitals, and 2. robustly determine optimal design for such orbitals.
We conclude that the answer to these questions is affirmative.
The rest of the paper is organized as follows. Section 2 describes the correlated multiresolution orbital solver. Section 3 describes the technical details of the computational experiments that are reported and analyzed in Section 4. Section 5 summarizes our findings and outlines future directions.

Formalism
Consider the energy functional of a (unit-normalized) state Ψ in a Fock-space state defined by M (spin-)orbitals {φ p }: where γ i j and γ ij kl are the (spin-orbital) 1-and 2-RDMs (see Appendix A for the notation and standard definitions). For a fixed model defining the map from orbitals {φ p } to |Ψ , Eq. To enforce the orbital orthonormality we introduce a Lagrangian, Variation φ * m → φ * m + δ produces the following first-order change in the Lagrangian: The optimal orbitals are determined by solving The multipliers are determined by projecting Eq. (4) onto the current (orthonormal) orbitals: This choice of the multipliers makes the gradient automatically orthogonal to the current orbitals.
If |Ψ is a single Slater determinant (SD) (e.g., a Hartree-Fock determinant) the RDMs simplify, and the gradient becomes: Its projection onto an arbitrary orbital φ a orthogonal to all current orbitals is vanishing of which with the optimal orbitals is the Brillouin condition.
For multiconfigurational states Eqs. (3) and (6) for the Lagrangian gradient and multipliers take the SD-like form when the 2-RDM is rewritten in terms of its cumulant (see

Appendix A). This nominally results in replacements
Clearly, Eq. (9) follows from Eq. (11) due to the vanishing 2-RDM cumulant for the SD state.
The use of natural orbitals (NOs), which make the 1-RDM diagonal, further simplifies the gradient and makes possible to design a simple orbital update strategy similar to that used for SD 36 and correlating orbitals 32 in a real-space representation: Partitioning the Fock operatorf into its free-particle componentd (the kinetic energy in the nonrelativistic case, or the Dirac operator in the relativistic case) and the restv, the orbital stationarity condition becomes: is the effective energy of NO m. The orbital update described below exploits the fact that˜ m < 0 for all orbitals, even weakly occupied. To understand why that is consider the following decomposition of the electronic energy (Eq. (1)) into the contributions from natural orbital m and the rest: Since the electronic energy cannot decrease in magnitude as the orbital basis is expanded, E m approximates the absolute decrease in the energy due to the exclusion of orbital m from the orbital set, hence this is why m m is nonnegative.
Updated orbitals are obtained by the action of Green's operator of particle in constant potential −˜ m > 0 on the r.h.s. of Eq. (15): For a nonrelativistic particle and open boundary conditions Green's operator has the familiar Yukawa kernel: extensions to the periodic boundary conditions 37 and Dirac particles 38 are also known. Since the SD energy is invariant with respect to orbital rotations, the orbital update Eq. (18) can be optimally simplified by switching to the canonical orbitals that diagonalize the Lagrangian coefficients (i.e., the Fock matrix in this case): which coincides with the orbital update used in Ref. 36.

Technical Details
The orbital solver has been implemented in a developmental version of the Massively Parallel The CI energies obtained with fixed Gaussian AO bases and the corresponding spectralelement orbitals obtained therefrom are shown in Table 1.
Unsurprisingly, the use of optimized numerical orbitals results in significant lowerings of the CI energy. Clearly, orbital optimization alone only addresses the suboptimality of the given orbital set of a fixed finite rank ("orbital error", quantified by E AO − E MRA ); the basis set incompleteness due to the finite rank ("rank error", quantified by E MRA − E exact ) is not addressed by the orbital optimization and requires increasing the rank (ideally, in conjunction with extrapolation or explicit correlation). Nevertheless, it is somewhat surprising how substantial is the the orbital errors of the correlation consistent basis sets (that are optimized for atomic CI energies 44 ) relative to their rank errors, shown in the second-to-last column in Table 1. Even for a quadruple-zeta basis orbital relaxation lowers the energy by 60% of the residual basis set error. These findings are in line with the unexpected sensitivity of the CC singles and doubles energy to the details of the unoccupied (pair-natural) orbitals that we found in the LCAO context. 45 What is more surprising, however, is that usually significant (> 50%) fraction of the energy lowering due to the orbital optimization can be attributed to the relaxation of the weakly-occupied (correlating) orbitals. To estimate this fraction the last column shows the the energy lowering due to relaxing the occupied Hartree-Fock (1s) orbital from its finite basis form to the exact form, as a percentage of the total CI energy lowering due to the orbital optimization. For example, with the cc-pVTZ AO basis only 33 % of the CI energy lowering can be attributed to the basis set incompleteness of the Hartree-Fock occupied orbital; this is somewhat unexpected considering that the correlation consistent basis family was optimized for atomic CI computations. This finding emphasizes that even numerical deficiencies of even weakly occupied orbitals can contribute significantly to the overall basis set errors of correlated energies. As Figure 1 illustrates, the changes in the natural orbitals due to the basis set relaxation can indeed be significant, not only near the nucleus (where Gaussian orbitals lack the cusp) and far from the nucleus (where Gaussians decay too quickly) but also at the intermediate distances from the nucleus where the orbital magnitude and positions of the radial nodes can be strongly impacted by the orbital optimization. Table 2, the orbital solver can robustly optimize orbital sets with occupancies differing by many orders of magnitude. As elaborated in Section 2, the effective orbital energy,˜ , that controls the orbital update via the Greens function (Eq. (18)) is negative even for the weakly-occupied NOs; in fact, they are even more negative that the energy of the strongly occupied NO. In contrast, the corresponding diagonal elements of the Fock operator are positive, as expected. The difference between the Fock and effective orbital energies is negligible for the strongly occupied NO.

Be
For the ground state of Be the energy lowering due to the orbital relaxation was similar to what was found for He (see Table Table 3. Microhartree-level agreement was found with the atomic CI solver in the state-of-the-art GRASP package. 46 Table 3: Converged CI energies (mE h ) of the ground-state He atom evaluated with Gaussian AO ("LCAO"), optimized real-space orbitals ("MRA"), and the atomic finite-difference MCSCF solver in the GRASP2K software ("FD"). 46 The LCAO and MRA energies are converged to all digits shown.

H 10
We chose H 10 in its ground state as a model system for testing solver robustness due to the availability of extensive benchmarks for a number of methods. 48 The system offers the ability to control the nature and strength of electron correlation by varying the interatomic distance and the size of the orbital basis. In particular, at short interatomic distances (R) even the smallest nonminimal (double-zeta) Gaussian AO basis becomes nearly linearly dependent and accurate evaluation of the Hamiltonian and the method-specific iterative solvers can become untenable. At short R the importance of angular correlation greatly increases compared to the equilibrium and stretched R. Lastly, the use of a minimal AO basis (1 orbital per atom) maximizes the degree of "on-site" repulsion and increases the extent of longitudinal correlation.
The initial investigation of our orbital solver focused on the maximally-compressed geometry considered in Ref. 48, where achieving the basis set limit with the standard Gaussian basis sets is difficult due to the rapid onset of ill-conditioning, which recently motivated developments of mixed numerical representations marrying Gaussian AOs with real-space approaches. 49 As expected, the orbital optimization greatly reduces the basis set incompleteness of the CI energy, as can be seen in the rightmost column of Table 4, by ×3.4 for the STO-6G basis an by ×7.4 for the cc-pVDZ basis. The estimated residual basis set error with the MRA-optimized cc-pVDZ orbitals is < 3mE h / atom; this is already comparable to the 1mE h uncertainty in the estimated nonrelativistic energy at this distance. 48 Table 4  NOs are similar: both exhibit the same number of nodes and the same symmetry with respect to the inversion center. As expected, and just like in the atoms, the real-space NOs exhibit correct nuclear cusps and slower asymptotic decay into the vacuum. Interestingly, however, real-space NOs have lower probability densities at most nuclei than their Gaussian counterparts.

Summary
We presented a general method for constructing orbitals that are optimal for computing energies with correlated (many-body methods) in arbitrary molecules and with robust control of numerical errors. In this work we employed multiresolution adaptive discontinuous spectral-element representation for the orbitals, but the presented method can be used with other real-space representations. Although utilized with the conventional and selected configuration interaction wave functions, the orbital solver can be used with arbitrary maps from orbitals to 1-and 2-particle reduced density matrices (e.g., tensor networks, stochastic methods, or quantum circuits). Numerical performance of the method was demonstrated by computing the CI energies of small atoms (matching the numerical precision of the atomic finite-difference CI codes) and a paradigmatic compressed polyatomic molecule; in all cases the use of real-space representation for the orbitals yielded energies that were significantly lower than the Gaussian AO-based counterparts, with correct asymptotic behavior at the nuclei and in vacuum.
Further algorithmic improvements should allow to speed up the rate determining steps in the orbital updates by 1 or 2 orders of magnitude even without employing accelerated hardware (such as general-purpose graphical processing units). The formalism should also map straightforwardly to the optimization of molecular spinors for relativistic correlated methods by borrowing the advances in 4-component relativistic numerical electronic structure introduced by some of us recently. 38 Work along these lines is already underway and will be reported elsewhere.

Acknowledgement
This work was supported by the U.S. Department of Energy via award DE-SC0022327. EFV is grateful to Greg Beylkin for providing the quadratures described in Ref. 50.

Appendix A Notation and Standard Definitions
In this work we largely follow a covariant tensor notation of the many-body quantum chemistry. 51,52 Einstein summation convention is implied, unless sums over any index are shown explicitly: namely, summation is implied over every symbol that appears once in a contravariant (upper, bra) position and once as in a covariant (lower, ket) position in a given tensor product, with the summation range defined by the symbol.
Matrix elements of one-, two-, and higher-body operators are denoted as follows: It is convenient to use similar notation for one-and higher-body functions resulting from omission of one or more bra indices (and the corresponding integrations) in the matrix elements:ô All many-body operators in this manuscript are particle-symmetric, i.e. o(1, 2) = o(2, 1), hence o p 1 p 2 q 1 q 2 = o p 2 p 1 q 2 q 1 . Matrix elements of the 1-and 2-particle parts of the Hamiltonian are denoted by h p q and g pq rs , respectively.
where γ s r denotes elements of 1-RDM (see Eq. (29)), are the familiar matrix elements of the Fock operator. The overlap of two functions is denoted by s p q . δ p q denotes the Kronecker delta: 2-and 1-RDMs for state |Ψ are defined as: withâ † i andâ i the standard creators and annihilators, respectively. If |Ψ is a multideterminantal/multiconfiguration (MC) state it is convenient to represent 2-RDM as with λ ij kl the 2-RDM cumulant. 53

Appendix B Spherical projection of orbital updates
In atomic computations it is useful to be able to obtain solutions with specific target number of subshells, e.g. 3s2p1d for a second-row atom. In general such solutions are saddle points; tracking such saddle points can fail if the gradient does not preserve the relevant symmetries exactly. Unfortunately the finite-precision spectral element bases employed in this work do not preserve spherical symmetry. Thus we need to impose symmetry restrictions on the orbital updates; without such restrictions the solutions may break spherical symmetry.
Consider a subshell of 2l + 1 "noisy" orbitals. To project out the components of the orbitals corresponding to other angular momenta we want to express each atomic orbital as a product of its radial component f (r) times the angular factor, i.e. as where the angular component is written as a unitary linear combination of spherical harmonics to support arbitrary orientation of each orbital. We determine f (r) by minimizing the average deviation of the subshell orbitals: The gradient with respect to f is: This suggests the following formula for f for the given fixed coefficients u: To determine nearly-optimal values of coefficients u for a given subshell orbital i we evaluate its overlaps with the spherical harmonics S (i) m,g ≡ φ i |Y lm (r g ) on a uniform radial grid {r g = gδ r , ∀g = 1..n g }. Here we set the grid size n g to 10 and the grid extent δ g (n g + 1) to 2π 2 T i , which is the de Broglie wavelength of a free particle with kinetic energy T i ≡ φ i |T |φ i . Note that these kinetic energies are available since detection of orbital subshells also used their kinetic energies. Coefficients u m i are given by the left-hand singular vector corresponding to the largest singular value of (S (i) ) m,g ≡ S m,g . Angular integration is performed via the lowest-order (72 points) rotationally invariant quadratures of Ahrens and Beylkin. 50