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

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,^{3, 4} or represented numerically⁵). However, (augmented) plane-wave (PW) representation of 1-particle states can also be used.^{6, 7} Unfortunately representation of many-particle by products of 1-particle states suffers from slow convergence in the vicinity of electron cusps,^{8, 9} thereby mandating the use of basis set extrapolation¹⁰ or the augmentation with explicit 2-particle basis (“explicit correlation”)^{11, 12, 13, 14}.

Although the LCAO-based many-body methods can be pushed to match experimental thermochemical¹⁵ and spectroscopic data¹⁶ for small molecules, extending such successes to increasingly larger systems, and to higher accuracy, is limited by the rapid onset of ill-conditioning, 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¹⁷ 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^{18, 19} or computing properties other than the energy.²⁰ 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 (i.e., the basis functions are only nonzero in a finite domain). 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^{21, 22, 23, 24, 25, 26, 27, 28, 29, 30} 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 many-body 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 first-quantized 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,^{32, 33, 34, 35} for molecules the feasibility of such representations has been demonstrated only recently by some of us (MP2)^{36, 37} and our collaborators (CIS(D), CC2)^{38, 39} by combining volume-element-wise pair function compression and regularization of electron cusps³⁶ similar to that performed in explicitly correlated R12/F12 methods.¹² It is possible to obtain energies and properties with such methods that are competitive to the state-of-the-art 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 (also referred to as the algebraic approximation⁴⁰). 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 multiconfiguration self-consistent field (MCSCF)^{41, 42, 43} and perturbative levels^{44, 45}) and diatomic molecules (MCSCF,^{46, 47, 43} perturbation theory⁴⁸ and coupled-cluster (CC)^{49, 50}); an excellent recent review of the vast relevant literature can be found in Ref. 51.

The first demonstration of any real-space orbital-based many-body method applicable to general molecules was shown recently by Kottmann, Bischoff, and Valeev⁵². The key to this advance was the realization that far more compact orbital representation of the weakly-occupied (virtual) orbitals needed for dynamical correlation is provided by the pair-natural orbitals.^{53, 54, 55} 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. 1.

   to express arbitrary (not just MBPT-based) orbital-based many-body wave functions in terms of real-space orbitals, and

2. 2.

   robustly determine optimal design for such orbitals.

We conclude that the answer to these questions is affirmative.

Note that much the formalism for optimizing correlated real-space orbitals presented in this work is generic and can be used with any real-space representation. The specific numerical representation that we used here is based on the multiresolution analysis (MRA) of a Hilbert space that leads to adaptive construction of scale-invariant real-space discountinuous high-order spectral element bases.⁵⁶ The resulting multiresolution adaptive spectral element representation (or, simply, MRA representation)^{57, 58, 59, 60} is usable for treating describing both valence and compact (core) electronic states with robust control of numerical error $\epsilon$. The multiscale structure of the representation results in low-order computational complexity for treating large systems. A number of electronic structure methods have been demonstrated in the multiresolution representation, including one-body (Hartree-Fock and DFT²³) and 2-body (MP2^{36, 37}, CC2³⁸) methods, with energies, forces, and other molecular properties^{61, 39, 62} of ground and excited states. The complete technical details of the numerical representation can be found in Refs. ^{59, 23, 63}; for a pedagogical introduction the reader is referred to Ref. 64.

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.

Consider the energy functional of a (unit-normalized) state $\Psi$ in a Fock-space state defined by $M$ (spin-)orbitals $\{\phi_{p}\}$:

where $\gamma^{i}_{j}$ and $\gamma^{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 $\{\phi_{p}\}$ to $\ket{\Psi}$, Eq. 1 defines a functional of the orbitals. Our objective is to minimize this functional with respect to $M$ orthonormal orbitals $\{\phi_{p}\}$. Unlike the traditional representations of orbitals in terms of a fixed basis, such as an atomic orbital basis in the LCAO formalism, in this work each orbital is expressed in a basis adaptively refined to approach the orbital’s exact representation in the full Hilbert space. The lack of a fixed basis precludes automatic enforcement of the orbital orthonormality by, e.g., unitary parametrization. Therefore it will be necessary to maintain orthonormality explicitly.

Here we follow the standard Lagrangian formalism for introducing orbital orthonormality constraints. Variation of the Lagrangian

by replacement $\phi_{m}^{*}\to\phi_{m}^{*}+\delta$ produces the following first-order change:

with $\hat{g}^{p}_{q}$ defined in Eq. 28. 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 $\ket{\Psi}$ 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 $\phi_{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 $h\to f$, $\gamma\to\lambda$:

Clearly, Eq. 9 follows from Eq. 11 due to the vanishing 2-RDM cumulant for the SD state.

Note that the single determinant expression for the gradient with respect to orbital $\phi_{m}$ (Eq. 9) only involves the action of the Fock operator onto $\phi_{m}$ itself. In contrast, the $\gamma^{m}_{i}\hat{f}\ket{\phi_{i}}$ contribution to the multideterminantal expression (15) involves the action of the Fock operator on every orbital. The use of natural orbitals (NOs), which make the 1-RDM diagonal,

allows to eliminate the kinematic coupling between the gradients:

Until now the formalism has been completely generic, i.e., applicable to any real-space representation. However, the manner in which Eq. 4 must be solved differs between numerical representation. Thus it is now necessary to specialize the formalism to the MRA numerical representation used in this work that we introduced in Section 1. Orbital optimization in the MRA representation in both SD²³ and correlated⁵² contexts uses the integral reformulation of the differential equations (this is crucial for being able to solve differential equations in discontinuous spectral element basis). Partitioning the Fock operator $\hat{f}$ into the free-particle Hamiltonian $\hat{d}$ and potential $\hat{v}$ (see Eq. 33), the orbital stationarity condition becomes:

where

is the effective energy of NO $m$. Note that the natural orbital basis is essential for the integral reformulation since it eliminates all but a single instance of the differential operator $\hat{d}$; solving for the natural orbitals is thus most sound strategy in the MRA numerical representation!

The orbital update described below exploits the fact that $\tilde{\epsilon}_{m}<0$ for all orbitals, even weakly occupied. This may appear counterintuitive, since we normally think of unoccupied orbitals as corresponding to positive energies. To rationalize this observation, consider decomposition of the electronic energy, Eq. 1, into the contributions from natural orbital $m$ and the rest:

Since the electronic energy becomes more negative as the orbital basis is expanded, $E_{m}$ approximates the absolute decrease in the energy due to the addition of orbital $m$ to the orbital set. A more rigorous mathematical analysis may be warranted, but the simple argument should suffice for now.

Updated orbitals are obtained by the action of Green’s operator of particle in constant potential $-\tilde{\epsilon}_{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⁶⁵ and Dirac particles⁶⁶ are also known. The updated orbitals are not orthonormal, thus the standard symmetric (Löwdin) orthonormalization⁶⁷ is applied to restore the orbital orthonormality while minimizing the induced change to facilitate extrapolation.⁶⁸

For highly-symmetric cases, like NOs with high angular momentum in atoms, it may be also desirable to project orbitals onto the corresponding invariant subspaces of the symmetry group. This is due to the spectral element basis lacking the desired symmetry properties, thus symmetries of the orbitals are maintained to precision $\epsilon$ used to define the basis; this is further compounded by the additional approximations invoked when evaluating Eq. 18. Appendix B described the relevant projection scheme used for the case of atoms where the problem can become particularly severe when low precision ($\epsilon>10^{-4}$) is used.

The rate of convergence of Eq. 18 is rapid, but it is possible to accelerate it further by extrapolation. To this end we deploy the Krylov-subspace Accelerated Inexact Newton (KAIN) method;⁶⁹ other strategies like the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method⁷⁰ should be usable as well. Note that to be able to make KAIN extrapolation robust it is essential to choose a canonical frame for the orbitals. Detailed description of the canonicalization approach are outside of the scope of this article (canonicalization is only essential for accelerating the solver) and will be reported elsewhere.

It is useful to point out how Eq. 18 is related to the orbital update used to determine optimal single-determinant (SD) orbitals in Ref. 23. For a single determinant the cumulant vanishes,

hence the orbitals are automatically natural and the gradients (9) are kinematically decoupled. The residual coupling in the SD version of (18) is via the Lagrange multipliers. Such coupling can be removed completely by switching to the canonical orbitals which make the Fock operator diagonal ($f^{m}_{n}=\delta^{m}_{n}\epsilon_{m}$) and simplify the orbital update to

This coincides with the standard orbital update used for single determinant SCF solver in MADNESS when orbitals are not localized; with localized orbitals the SD version of Eq. 18 (with all cumulant contributions omitted) is used instead. Note that in all of these approaches the Lagrange multipliers are evaluated every iteration directly (as matrix elements of the Fock operator in the basis of current orbitals) rather than iteratively updated as in the original solver described in Ref. 23.

Specialization of Eqs. 12 and 18 for the important case of the complete active space is straightforward, namely for inactive orbitals the cumulant contributions vanish. No additional changes are needed.

The orbital solver has been implemented in a developmental version of the Massively Parallel Quantum Chemistry (MPQC) package⁷¹ using the implementation of the MRA calculus in the MADNESS library.⁶³ The solver is bootstrapped with Gaussian AOs taken from the user-provided basis and projected onto the spectral element basis. Mean-field orbitals are then obtained by optimization within the projected AO basis. RDMs were obtained from 1- and 2-electron integrals evaluated in the mean-field (SCF) orbitals expressed in real-space basis using a Heat-Bath CI (HCI)⁷² solver developed in MPQC. All computations on atoms and the minimal-basis H₁₀ computation set HCI selection thresholds $\epsilon_{1,2}$ to zero; the H₁₀ computation bootstrapped from the cc-pV{D,T}Z bases used $\epsilon_{\{1,2\}}=10^{\{-4,-6\}}$.

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 a 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.

Despite using the efficient production-grade implementation of the MRA calculus in MADNESS, the current implementation is not fully optimal. Although as we demonstrated here optimization of more than 100 correlated orbitals is already feasible with the current (CPU-only) code on a single multicore node, the implementation of the time-determining step in the orbital update, namely, the cumulant-dependent contribution to Eq. 18 (evaluated at the $\mathcal{O}(n_{\mathrm{act}}^{4})$ asymptotic cost, with $n_{\mathrm{act}}$ the number of active orbitals), does not explicitly optimize for CPU caches and can be significantly improved (by one to two orders of magnitude) even without exploiting accelerators. The optimization of $~{}140$ correlated orbitals in H₁₀ reported in Section 4.2.2 took less than 4 hours per iteration on a single node of the TinkerCliffs cluster at Virginia Tech’s Advanced Research Computing equipped with a 128-core AMD EPYC 7702 processor; this was already comparable with $~{}30$ minutes per iteration that the selected CI 2-RDM evaluation took.

The implementation can also be improved in the other asymptotic regime, namely for applications with a modest number of active orbitals ($n_{\mathrm{a}}\leq 50$) and a large number of inactive orbitals, $n_{\mathrm{i}}$. The asymptotic cost is then dominated by the evaluation of the exchange contribution to the $\hat{v}\ket{\phi_{m}}$ term on the right-hand side of Eq. 18. By orbital localization of the inactive orbitals fast application of the exchange operator will be possible, as demonstrated by the past development of the optimally-scaling MRA Hartree-Fock method⁸⁴. In the current implementation the inactive orbitals are not localized, hence the cost of the exchange operator application is at worst $\mathcal{O}(n^{3}_{\mathrm{c}})$. Nevertheless practical applications to general systems can be readily performed without these optimizations. For example, a (4,4)-CASSCF($\epsilon=10^{-5})$ evaluation of the ground state energy of the butadiene molecule took 380 seconds on an Apple MacBook Pro with a 12-core Apple M2 Max chip; for comparison, the state-of-the-art Hartree-Fock MRA solver in MADNESS took 130 seconds. We expect the described improvements to bring down the cost of the CAS computation with small active spaces close to the cost of the optimized MRA Hartree-Fock implementation. The conventional (4,4) CASSCF wave function of comparable quality requires a cc-pV5Z basis (see Table 5); such computation with BAGEL⁸⁰ (using the matching cc-pV5Z-JKFIT density fitting basis) on the same machine took 186 seconds.

After the outlined algorithmic improvements we expect the cost of orbital optimization to become largely trivial relative to the wave function/density optimization. We are thus hopeful that the real-space representations will become competitive with or even preferred to the AO representations for correlated models.

Our formalism should 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.⁶⁶ Extension to periodic systems is also straigtforward by using the existing extension of Poisson and Helmholtz Green’s operators to periodic boundary conditions.⁶⁵ Work along these lines is already underway and will be reported elsewhere.

This work was supported by the U.S. Department of Energy via award DE-SC0022327. We are grateful to Greg Beylkin for providing the quadratures described in Ref. 86. The authors acknowledge Advanced Research Computing at Virginia Tech (https://arc.vt.edu/) for providing computational resources and technical support that have contributed to the results reported within this paper.