Direct determination of optimal pair-natural orbitals in a real-space representation: The second-order Moller–Plesset energy

Kottmann, Jakob S.; Bischoff, Florian A.; Valeev, Edward F.

doi:10.1063/1.5141880

An efficient representation of molecular correlated wave functions is proposed, which features regularization of the Coulomb electron–electron singularities via the F12-style explicit correlation and a pair-natural orbital factorization of the correlation components of the wave function expressed in the real space. The pair-natural orbitals are expressed in an adaptive multiresolution basis and computed directly by iterative variational optimization. The approach is demonstrated by computing the second-order Moller–Plesset energies of small- and medium-sized molecules. The resulting MRA-PNO-MP2-F12 method allows for the first time to compute correlated wave functions in a real-space representation for systems with dozens of atoms (as demonstrated here by computations on alkanes as large as C₁₀H₂₂), with precision exceeding what is achievable with the conventional explicitly correlated MP2 approaches based on the atomic orbital representations.

I. INTRODUCTION

Predictive computation of properties of molecules and materials requires a highly precise numerical representation of many-body electronic wave functions/operators. Accurate electronic wave functions, such as in the coupled-cluster method, are traditionally represented in the Linear Combination of Atomic Orbitals (LCAO) representation using pre-optimized sequences of atom-centered basis sets (represented by analytic Gaussian-type or Slater-type functions, or represented numerically¹). Recently, many-body electronic structure computations have also been demonstrated in the (augmented) plane-wave (PW) representation.² Both representations suffer from slow convergence to the basis set limit, requiring the use of basis set extrapolation^3,4 or explicit correlation.^5,6 Although it is possible to closely match experimental chemical energy differences, such as atomization energies of small molecules,⁷ with the LCAO approaches, such computations suffer from high cost and poor conditioning and do not lend themselves to fast (reduced-scaling) reformulation due to the fully dense representation of electronic states and operators. The known bias of common AO basis sets toward ground-state energies makes it difficult to attain similar accuracy for excited states and response properties. Meanwhile, the plane wave representation struggles with the description of position-localized (e.g., intra-atomic) features of the electronic wave functions and thus also do not naturally lend themselves to low-order formulations.

Grid-based real-space numerical representations of many-body wave functions, such as finite difference and finite/spectral elements, provide an alternative to the LCAO and PW representations with a number of attractive features, such as the ability to resolve spatially localized features, systematic bias-free improvability, and good conditioning. The main bottleneck to routine use of such representations is the impractical size of the 3k-dimensional grid needed to represent a k-particle state. Even pair theories (k = 2) such as MP2 and CCSD become prohibitively expensive for molecular systems due to the grid size; only symmetry-based dimension reduction in atoms makes the application of pair theories possible.^8–12

Recently, we demonstrated for the first time real-space 2-body electronic structure methods (MP2, CIS(D), and CC2), in general, molecular systems without any use of symmetry by exploiting (1) multiresolution adaptive spectral element representation of the 2-electron (pair) functions, (2) re-formulation of the integrodifferential pair function equations in the integral form, (3) tensor compression of pair function coefficients, and (4) regularization of electron–electron and electron–nuclear Coulomb singularities via explicit electron–electron and electron–nuclear correlation.^13–18 Let us re-emphasize that our approach is distinct from other notable recent uses of grids in many-body electronic structures^19,20 in that the LCAO representation is not used at all. Our real-space approach achieved supremacy to the best available LCAO-F12 counterparts by producing correlation energies and other molecular properties with smaller numerical errors; this suggests that when the utmost numerical precision is necessary, the real space representation becomes the preferred choice.

Despite the promising numerical performance, real space many-body wave functions are costly to compute due to the sheer size of the 2-particle basis, even with the adaptive basis refinement, coefficient compression, and operator regularization. Here, we propose to address this challenge by leveraging a lesson of the reduced scaling LCAO many-body technology based on pair natural orbitals,^21–30 namely, that pair functions are globally separable with a low (≪100) separation rank. The initial demonstration will be limited to the ground-state MP2 method for simplicity, but extensions to higher-order methods, and to excited states and properties, will be straightforward. The key novel feature of the real-space MRA-PNO-MP2-F12 method is the direct optimization of the pair natural orbitals (PNOs) represented in a real-space representation and the adaptive determination of the PNO subspace ranks.

The rest of this paper is organized as follows. Section II briefly recaps the application of the multiresolution analysis (MRA) to the many-body electronic structure methods as well as the proposed MRA-PNO-MP2-F12 method. Section III describes the technical details of the computations whose results are reported and analyzed in Sec. IV. Our findings are summarized in Sec. V.

II. FORMALISM

A. Multiresolution spectral element representation

The multiresolution analysis (MRA) of a Hilbert space defines its recursive decomposition that translates into a straightforward adaptive construction of scale-invariant real-space spectral element bases.³¹ The multiresolution adaptive spectral element representation (or, simply, MRA representation) allows a compact representation of sufficiently smooth functions with arbitrary guaranteed precision ϵ, sparse representation of many differential and integral operators, and fast algorithms for their applications,^32–35 thus realizing a basis-free computational integrodifferential calculus. A number of electronic structure methods have been demonstrated in the multiresolution representation, including one-body (Hartree–Fock and DFT³⁶) and two-body (MP2¹³ and CC2¹⁸) methods, with energies, forces, and other molecular properties^17,38,39 of ground and excited states. Complete details can be found in Refs. 34, 36, 40, and 41.

B. MRA-MP2 and MRA-PNO-MP2-F12

1. MRA-MP2

The closed-shell second-order Moller–Plesset energy can be obtained by minimizing its Hylleraas functional⁴² with respect to the first-order pair functions τ_ij(r₁, r₂) ≡ ⟨r₁r₂|τ_ij⟩,

\begin{matrix} L^{(2)} [τ] & = & \sum_{i j} ⟨) τ_{i j} | (2 - {\hat{P}}_{12}) {\hat{Q}}_{12} \\ \times \{2 g_{12} | i j (⟩ + {\hat{F}}_{12} | τ_{i j} (⟩ - \sum_{k} (| τ_{k j} (⟩ F_{k i} + | τ_{i k} (⟩ F_{k j})\}, \end{matrix}

(1)

where ${\hat{P}}_{12}$ permutes electrons 1 and 2, $F_{i j} \equiv ⟨) i | \hat{F} | j (⟩$ is the matrix element of the Fock operator $\hat{F}$ ⁠, ${\hat{F}}_{12} = {\hat{F}}_{1} + {\hat{F}}_{2}$ ⁠, and $g_{12} \equiv r_{12}^{- 1}$ is the electron–electron repulsion. The strong orthogonality projector ${\hat{Q}}_{12}$ is defined as

{\hat{Q}}_{12} \equiv {\hat{Q}}_{1} {\hat{Q}}_{2} = (1 - {\hat{O}}_{1}) (1 - {\hat{O}}_{2}),

(2)

with $\hat{O} \equiv \sum_{m} | m (⟩ ⟨) m |$ being the projector onto the subspace spanned by the occupied orbitals of the reference determinant.³⁷ For real pair functions, the stationarity condition of the Hylleraas functional can be compactly written as

\begin{matrix} 0 & = & \frac{1}{6} (2 + {\hat{P}}_{12}) \frac{δ L^{(2)} [τ]}{δ τ_{i j}} \\ = & {\hat{F}}_{12} | τ_{i j} (⟩ + {\hat{Q}}_{12} g_{12} | i j (⟩ - \sum_{k} (| τ_{i k} (⟩ F_{k j} + | τ_{k j} (⟩ F_{i k}) \end{matrix}

(3)

with the functional derivative with respect to the pair functions given by

\begin{matrix} \frac{δ L^{(2)} [τ]}{δ τ_{i j}} & = & (2 - {\hat{P}}_{12}) (2 {\hat{Q}}_{12} g_{12} | i j (⟩ + 2 {\hat{F}}_{12} | τ_{i j} (⟩) \\ - (2 \sum_{k} (| τ_{k j} (⟩ F_{k i} + | τ_{i k} (⟩ F_{k j})) . \end{matrix}

(4)

The real-space projection of Eq. (3) is an integrodifferential equation that can be solved by a number of methods; here, we follow the standard approach of Refs. 13, 14, 17, and 18 by converting it to the integral form,

| τ_{i j} (⟩ = - 2 {\hat{G}}_{μ} \{{\hat{Q}}_{12} g_{12} | i j (⟩ + {\hat{v}}^{F}_{12} | τ_{i j} (⟩ - \sum_{k \neq j} | τ_{i k} (⟩ F_{k j} - \sum_{k \neq i} | τ_{k j} (⟩ F_{i k}\} .

(5)

Here, we introduced a shorthand notation for the potential energy component of the Fock operator,

{\hat{v}}^{F}_{12} \equiv {\hat{v}}^{F}_{1} + {\hat{v}}^{F}_{2},

(6)

{\hat{v}}^{F} \equiv \hat{F} - \hat{T} = V_{nuc} + \hat{J} - \hat{K},

(7)

where $\hat{J}$ and $\hat{K}$ are the usual Coulomb and exchange operators and ${\hat{G}}_{μ} \equiv {({\hat{T}}_{12} - μ^{2})}^{- 1}$ is the 6-dimensional Helmholtz Green’s operator for free-space bound states^36,44 parameterized by

μ = \sqrt{- 2 (F_{i i} + F_{j j})} .

(8)

While for some methods the initial guess has to be constructed by nontrivial approaches, here, the trial real-space pair function can be set to zero. The coupled system of integral equations is then iterated until the fixed point is reached to the desired precision.

2. MRA-PNO-MP2-F12

In this work, we will solve the first-order MP equations by representing the pair functions as products of pair natural orbitals (PNOs). First introduced by Edmiston and Krauss²¹ and developed by Meyer, Kutzelnigg, Ahlrichs, and others,^23–25 these were recently reintroduced by Neese as a means to reduce the cost of reduced-scaling coupled-cluster methods in the LCAO representation.^26–30 Here, we will extend the PNO-based MP2 approach as a means to make the real-space representation practical for many-body methods.

The PNOs in the LCAO approach are usually computed by diagonalizing pair contributions to the unoccupied block of the correlation contribution to 1-RDM evaluated from approximate (uncoupled) pair functions. Once determined, the PNOs are traditionally fixed, although recently one of us discovered that PNO optimization can offer a substantial improvement in the context of coupled-cluster singles and doubles.⁴⁵ In a real-space representation, it is not feasible to form 1-RDM (a function of 2 particle coordinates), considering the large number of parameters needed to represent even 1-particle functions. Thus, we will instead determine the PNOs directly by minimizing the second-order Hylleraas functional expressed in terms of PNO-expanded pair functions.

Each first-order pair function |τ_ij⟩ will be expanded as a linear combination of all possible products of its PNOs {|c_ij⟩},

| τ_{i j} (⟩ \overset{PNO}{=} \sum_{c_{i j} d_{i j}} t_{c_{i j} d_{i j}} | c_{i j} d_{i j} (⟩ .

(9)

The second-order Hylleraas functional thus depends on the PNOs and the corresponding amplitudes,

\begin{align} L^{(2)} [a, t] & = \sum_{i j} \sum_{c_{i j} d_{i j}} t_{c_{i j} d_{i j}} ⟨) c_{i j} d_{i j} | (2 - {\hat{P}}_{12}) \{2 {\hat{Q}}_{12} g_{12} | i j (⟩ + \sum_{e_{i j} f_{i j}} {\hat{F}}_{12} | e_{i j} f_{i j} (⟩ t_{e_{i j} f_{i j}}) \\ (- \sum_{k} (\sum_{e_{k j} f_{k j}} | e_{k j} f_{k j} (⟩ F_{k i} t_{e_{k j} f_{k j}} + \sum_{e_{i k} f_{i k}} | e_{i k} f_{i k} (⟩ F_{k j} t_{e_{i k} f_{i k}})\} \\ = \sum_{i j} \sum_{c_{i j} d_{i j}} {\tilde{t}}_{c_{i j} d_{i j}} \{2 ⟨) c_{i j} d_{i j} | g_{12} | i j (⟩ + \sum_{e_{i j}} t_{c_{i j} e_{i j}} F_{d_{i j} e_{i j}} + \sum_{e_{i j}} t_{e_{i j} d_{i j}} F_{c_{i j} e_{i j}}) \\ (- \sum_{k} (\sum_{e_{k j} f_{k j}} S_{c_{i j} e_{k j}} S_{d_{i j} f_{k j}} t_{e_{k j} f_{k j}} F_{k i} + \sum_{e_{i k} f_{i k}} S_{c_{i j} e_{i k}} S_{d_{i j} f_{i k}} t_{e_{i k} f_{i k}} F_{k j})\}, \end{align}

(10)

where the usual contravariant amplitudes were introduced,

{\tilde{t}}_{a b} \equiv 2 t_{a b} - t_{b a} .

(11)

Note that the overlaps S_ab ≡ ⟨a|b⟩ in Eq. (10) appear since the PNOs for different pairs are not orthogonal to each other.

The stationarity conditions given with respect to both quantities can be written in compact form as

\begin{align} 0 & = \frac{1}{6} (2 + {\hat{P}}_{a_{i j} b_{i j}}) (\frac{\partial L^{(2)} [a, t]}{\partial t_{a_{i j} b_{i j}}}) \\ = ⟨) a_{i j} b_{i j} | g_{12} | i j (⟩ + \sum_{c_{i j}} (t_{a_{i j} c_{i j}} F_{c_{i j} b_{i j}} + t_{c_{i j} b_{i j}} F_{c_{i j} a_{i j}}) \\ - \sum_{k} (\sum_{c_{i k} d_{i k}} S_{a_{i j} c_{i k}} S_{b_{i j} d_{i k}} t_{c_{i k} d_{i k}} F_{k j} + \sum_{c_{k j} d_{k j}} S_{a_{i j} c_{k j}} S_{b_{i j} d_{k j}} t_{c_{k j} d_{k j}} F_{k i}), \end{align}

(12)

\begin{align} 0 & = |\frac{δ L^{(2)}}{δ a_{i j}}⟩ = 2 \hat{Q} \sum_{c_{i j}} ({\tilde{t}}_{a_{i j} c_{i j}} g_{j}^{c_{i j}} | i (⟩ + {\tilde{t}}_{c_{i j} a_{i j}} g_{i}^{c_{i j}} | j (⟩) \\ + 2 \sum_{c_{i j} d_{i j}} \hat{F} | c_{i j} (⟩ ({\tilde{t}}_{c_{i j} d_{i j}} t_{a_{i j} d_{i j}} + {\tilde{t}}_{d_{i j} c_{i j}} t_{d_{i j} a_{i j}}) \\ + 2 \sum_{c_{i j} d_{i j} e_{i j}} | c_{i j} (⟩ F_{d_{i j} e_{i j}} ({\tilde{t}}_{c_{i j} d_{i j}} t_{a_{i j} e_{i j}} + {\tilde{t}}_{d_{i j} c_{i j}} t_{e_{i j} a_{i j}}) \\ - 2 \sum_{k} \sum_{c_{i j} e_{i k} f_{i k}} | e_{i k} (⟩ S_{c_{i j} f_{i k}} F_{k j} ({\tilde{t}}_{c_{i j} a_{i j}} t_{f_{i k} e_{i k}} + {\tilde{t}}_{a_{i j} c_{i j}} t_{e_{i k} f_{i k}}) \\ - 2 \sum_{k} \sum_{c_{i j} e_{k j} f_{k j}} | e_{k j} (⟩ S_{c_{i j} f_{k j}} F_{i k} ({\tilde{t}}_{c_{i j} a_{i j}} t_{f_{k j} e_{k j}} + {\tilde{t}}_{a_{i j} c_{i j}} t_{e_{k j} f_{k j}}), \end{align}

(13)

where

g_{k}^{l} (r) \equiv \int d^{3} r \frac{1}{| r - r^{'} |} k (r^{'}) l (r^{'}) .

(14)

The amplitude gradient (12) is evaluated similar to Eq. (12) of Ref. 28 where the PNOs for a given pair were forced to be mutually orthonormal, hence leading to a slightly different equation. In this work, the mutual orthogonality of the PNOs for a given pair is optional and can be ensured by diagonalizing the pair specific Fock matrices in each iteration. Note that Ref. 28 used a fixed (pair-specific) basis of projected atomic orbitals to represent the PNOs, while in this work, the basis used to represent the PNO changes from iteration to iteration. The PNO gradient (13) can be further simplified to

\begin{align} |\frac{δ L}{δ a_{i j}}⟩ & = 2 \hat{Q} \sum_{c_{i j}} ({\tilde{t}}_{a_{i j} c_{i j}} g_{j}^{c_{i j}} | i (⟩ + {\tilde{t}}_{c_{i j} a_{i j}} g_{i}^{c_{i j}} | j (⟩) + \sum_{c_{i j}} (\hat{F} | c_{i j} (⟩ D_{c_{i j} a_{i j}} + | c_{i j} (⟩ E_{c_{i j} a_{i j}}) \\ - \sum_{k} (\sum_{e_{i k}} | e_{i k} (⟩ F_{k j} D_{e_{i k} a_{i j}} + \sum_{e_{k j}} | e_{k j} (⟩ F_{i k} D_{e_{k j} a_{i j}}) \end{align}

(15)

by introducing the PNO–PNO block of the cross-pair reduced density matrix,

D_{a_{i j} b_{k j}} \equiv 2 \sum_{c_{i j} d_{k j}} S_{c_{i j} d_{k j}} ({\tilde{t}}_{a_{i j} c_{i j}} t_{b_{k j} d_{k j}} + {\tilde{t}}_{c_{i j} a_{i j}} t_{d_{k j} b_{k j}}),

(16)

which for the same-pair case (i.e., k = i) reduces, due to the orthogonality between the PNOs of a given pair, to the usual PNO 1-RDM matrix,²⁸

D_{a_{i j} b_{i j}} = 2 \sum_{c_{i j}} ({\tilde{t}}_{a_{i j} c_{i j}} t_{b_{i j} c_{i j}} + {\tilde{t}}_{c_{i j} a_{i j}} t_{c_{i j} b_{i j}}) .

(17)

We also introduced the energy-weighted variant of the 1-pair RDM, given by

E_{a_{i j} b_{i j}} \equiv 2 \sum_{c_{i j} d_{i j}} F_{c_{i j} d_{i j}} ({\tilde{t}}_{a_{i j} c_{i j}} t_{b_{i j} d_{i j}} + {\tilde{t}}_{c_{i j} a_{i j}} t_{d_{i j} b_{i j}}) .

(18)

C. Explicitly correlated ansatz for MRA-PNO-MP2

The first-order pair functions in real space contain the 3-dimensional seam of the electron–electron cusp caused by the singular electron–electron interaction.⁴³ While this leads to the well-known slow asymptotic decay of the basis set error in LCAO,⁴⁶ the cusp seam results in an impractical number of parameters (e.g., grid size) in a real-space representation of the pair functions. Regardless of a representation, it is mandatory to describe the known analytic structure of the cusp by an explicit contribution to the pair function. We thus introduce the explicitly correlated R12/F12 ansatz for the pair function^5,47

| τ_{i j} (⟩ = | u_{i j} (⟩ + | w_{i j} (⟩

(19)

with the explicitly correlated cusp-containing term given by

| w_{i j} (⟩ = \hat{Q} f_{12} | i j (⟩

(20)

and |u_ij⟩ is a cusp-free electron pair function. The explicitly correlated ansatz regularizes the MP1 equations by canceling the singularity of the Coulomb potential with the kinetic energy operator acting on $w$ _ij.^13,16,18

With the ansatz of Eq. (19), the Hylleraas functional is split into one functional depending only on the regularized pair functions u_ij, one depending only on the explicitly correlated residues $w$ _ij, and a coupling functional,

L^{(2)} [τ] = L^{(2)} [u] + L^{(2)} [w] + {L^{(2)}}_{c} [u, w] .

(21)

The regularized pair functions u are expanded into the PNOs as in Eq. (9); henceforth, all PNOs and amplitudes will refer to the regularized pair functions u.

For the evaluation of the functionals depending on $w$ ⁠, we use the regularized potential ${\tilde{g}}_{12}$ defined by¹⁸

[{\hat{F}}_{1} + {\hat{F}}_{2}, f_{12}] = {\hat{U}}_{12} - [{\hat{K}}_{1} + {\hat{K}}_{2}, f_{12}] - g_{12} \equiv {\tilde{g}}_{12} - g_{12} .

(22)

In the coupling functional ${L^{(2)}}_{c} [u, w]$ and the Hylleraas functional that depends on $w$ only, $L^{(2)} [w]$ ⁠, all coupling terms between the individual electron pairs cancel out, leaving only

\begin{align} L^{(2)} [w] & = \sum_{i j} (4 ⟨) w_{i j} | g_{12} | i j (⟩ - 2 ⟨) w_{i j} | g_{12} | j i (⟩) \\ + (2 ⟨) w_{i j} | ({\tilde{g}}_{12} - g_{12}) | i j (⟩ - ⟨) w_{i j} | ({\tilde{g}}_{12} - g_{12}) | j i (⟩) \\ = \sum_{i j} (2 ⟨) w_{i j} | ({\tilde{g}}_{12} + g_{12}) | i j (⟩ - ⟨) w_{i j} | ({\tilde{g}}_{12} + g_{12}) | j i (⟩), \end{align}

(23)

\begin{align} {L_{c}}^{(2)} [u, w] & = 2 \sum_{i j} (⟨) u_{i j} | ({\tilde{g}}_{12} - g_{12}) | i j (⟩ - ⟨) u_{i j} | ({\tilde{g}}_{12} - g_{12}) | j i (⟩) \\ = \sum_{i j} \sum_{a_{i j} b_{i j}} {\tilde{t}}_{a_{i j} b_{i j}} ⟨) a_{i j} b_{i j} | ({\tilde{g}}_{12} - g_{12}) | i j (⟩ . \end{align}

(24)

The regularized Hylleraas functional $L^{(2)} [u]$ takes the same form as in Eq. (10). If the coupling functional is combined with the Hylleraas functional of the regularized pair functions, the Coulomb parts will cancel, resulting in the same functional as before with the Coulomb potential replaced by the regularized potential ${\tilde{g}}_{12}$ ⁠,

L^{(2)} [τ] = {\tilde{L}}^{(2)} [u] + L^{(2)} [w],

(25)

where we introduced the regularized PNO-MP2-F12 Hylleraas functional,

\begin{align} {\tilde{L}}^{(2)} [u] & = L^{(2)} [u] + {L^{(2)}}_{c} [u, w] \\ = \sum_{i j} \sum_{a_{i j} b_{i j}} (2 t_{a_{i j} b_{i j}} - t_{b_{i j} a_{i j}}) [) \sum_{c_{i j}} t_{c_{i j} b_{i j}} F_{a_{i j} c_{i j}} + \sum_{c_{i j}} t_{a_{i j} c_{i j}} F_{b_{i j} c_{i j}} \\ - \sum_{k} (\sum_{c_{i k} d_{i k}} \sum t_{c_{i k} d_{i k}} S_{a_{i j} c_{i k}} S_{b_{i j} d_{i k}} F_{k j} + \sum_{c_{k j} d_{k j}} t_{c_{k j} d_{k j}} S_{a_{i j} c_{k j}} S_{b_{i j} d_{k j}} F_{k i}) \\ + 2 ⟨) a_{i j} b_{i j} | {\tilde{g}}_{12} | i j (⟩]) . \end{align}

(26)

Since the $L^{(2)} [w]$ does not depend on PNOs or their amplitudes, the explicitly correlated equations for the PNOs and PNO amplitudes are identical to the conventional counterparts [Eqs. (12) and (15)], except the Coulomb potential g₁₂ is replaced by the regularized potential ${\tilde{g}}_{12}$ ⁠.

1. MP2-F12 energies

Having determined the regularized pair functions, the MP2-F12 energy can be evaluated by the Hylleraas functional

L^{(2)} [τ] = L^{(2)} [u] + L^{(2)} [w] + L_{c}^{(2)},

(27)

\begin{align} L^{(2)} [w] & = \sum_{i j} ⟨) w_{i j} | (2 - {\hat{P}}_{12}) ({\tilde{g}}_{12} + g_{12}) | i j (⟩ \\ = \sum_{i j} ⟨) i j | f_{12} {\hat{Q}}_{12} (2 - {\hat{P}}_{12}) ({\tilde{g}}_{12} + g_{12}) | i j (⟩, \end{align}

(28)

L_{c}^{(2)} = \sum_{i j} \sum_{a_{i j} b_{i j}} {\tilde{t}}_{a_{i j} b_{i j}} ⟨) a_{i j} b_{i j} | ({\tilde{g}}_{12} - g_{12}) | i j (⟩ .

(29)

The pure f₁₂ dependent functional $L^{(2)} [w]$ of Eq. (28) will contain terms such as

\begin{align} ⟨) i j | f_{12} {\hat{Q}}_{12} {\hat{K}}_{1} f_{12} | i j (⟩ \\ = ⟨) i j | f_{12} {\hat{K}}_{1} f_{12} | i j (⟩ + ⟨) i j | f_{12} ({\hat{O}}_{1} {\hat{O}}_{2} - {\hat{O}}_{1} - {\hat{O}}_{2}) {\hat{K}}_{1} f_{12} | i j (⟩ \end{align}

(30)

emerging from the regularized potential $\tilde{g}$ ⁠. The second term on the right-hand side of Eq. (30) can be straightforwardly evaluated without the need to represent 2-particle quantities. This is not the case for the first term,

\begin{align} ⟨) i j | f_{12} {\hat{K}}_{1} f_{12} | i j (⟩ & = \sum_{k} \int d^{3} r_{1} ϕ_{i}^{*} (r_{1}) ϕ_{k} (r_{1}) \int d^{3} r_{2} ϕ_{j}^{*} (r_{2}) f_{12} ϕ_{j} (r_{2}) \\ \times \int d^{3} r_{3} ϕ_{k}^{*} (r_{3}) g_{13} f_{32} ϕ_{i} (r_{3}) . \end{align}

(31)

One way to evaluate such terms within an MRA framework is to compute a two-particle intermediate function

f (r_{1}, r_{2}) = \int d^{3} r_{3} ϕ_{k}^{*} (r_{3}) g_{13} f_{32} ϕ_{i} (r_{3}),

(32)

making the costs of the MRA-PNO-MP2-F12 computation comparable to the full (non-PNO) MRA-MP2-F12 computation and therefore shares its disadvantages, namely, the need to compute and store a 6-dimensional function. Another possibility is to introduce an auxiliary virtual basis set $\{| v_{k} (⟩\}$ in order to approximate $\hat{Q} \approx \sum_{k} | v_{k} (⟩ ⟨) v_{k} |$ and reduce Eq. (30) to an expectation value over one-particle functions

\begin{align} ⟨) i j | f_{12} {\hat{Q}}_{12} {\hat{K}}_{1} f_{12} | i j (⟩ = & \sum_{k} ⟨) i j | f_{12} | v_{k} (1) (⟩ {\hat{Q}}_{2} ⟨) v_{k} (1) | {\hat{K}}_{1} f_{12} | i j (⟩ \\ = & \sum_{k} ⟨) j | f_{v_{k}}^{i} \hat{Q} f_{i}^{{\tilde{v}}_{k}} | j (⟩ with | {\tilde{v}}_{k} (⟩ \equiv \hat{K} | v_{k} (⟩, f_{v}^{i} \equiv ⟨) i | f_{12} | v (⟩ . \end{align}

(33)

Note that it is not necessary to approximate the full two-particle projector ${\hat{Q}}_{12}$ ⁠. The introduction of such an auxiliary basis set is possible and can be performed in reasonable time compared to the optimization of the PNOs, but it introduces an additional approximation.

To avoid these challenges, we compute the energy by projection

E^{MP 2} [τ] = \sum_{i j} ⟨) i j | (2 - {\hat{P}}_{12}) g_{12} | τ_{i j} (⟩

(34)

such that terms such as Eq. (28) do not appear anymore and the energy can be calculated as

E^{MP 2 - F 12} [u, w] = E^{MP 2} [u] + E^{MP 2} [w],

(35)

E^{MP 2} [u] = \sum_{i j} \sum_{a_{i j} b_{i j}} {\tilde{t}}_{a_{i j} b_{i j}} ⟨) i j | g_{12} | a_{i j} b_{i j} (⟩,

(36)

E^{MP 2} [w] = \sum_{i j} ⟨) i j | (2 - {\hat{P}}_{12}) g_{12} {\hat{Q}}_{12} f_{12} | i j (⟩ .

(37)

In LCAO-MP2 without explicit correlation, the projection formula and Hylleraas functional yield identical working equations for the energy. The two approaches differ if the MP1 wave function is not solved exactly, as it is the case in LCAO-MP2-F12. Since MRA-MP2-F12 and MRA-PNO-MP2-F12 compute the MP1 (residual) wave function numerically exactly, the two expressions again yield the same energy once the wave function is converged with respect to the number of PNOs. Another argument for using projected MP2 energies is that they are similar to the projected coupled-cluster energies so that the quality of the optimized PNOs when used in the coupled-cluster context can be better estimated by the projected MP2 energies. Note, however, that the projected MP2 energies are not variational (in the sense of the Hylleraas functional) so that the final MRA-PNO-MP2-F12 energy can be lower than the exact (basis set limit) MP2 energy. In addition, if the PNOs are far from being converged, the Hylleraas functional and the projected energies can differ quite significantly. If not explicitly stated otherwise, we will always show projected energies in this article.

D. Optimizing PNOs and their amplitudes

1. Optimization of the PNO amplitudes

At each step of the PNO optimization procedure (see Fig. 2), the amplitudes are optimized for the current PNOs. The straightforward way to do so is by solving Eq. (12), but this might cause convergence problems, especially for larger systems with localized orbitals. A solution to this is the canonicalization of the PNOs, i.e., rotating them in a way that the PNO–PNO Fock matrix $F_{a_{i j} b_{i j}}$ for each pair becomes diagonal and the coupling between the PNOs of the pair vanishes. Note that if localized orbitals are used, the coupling between the different pairs is still present so that the amplitudes still need to be solved in an iterative way. It is also possible to use a quadratic optimization scheme to improve the convergence behavior with or without reduced coupling; here, we used the KAIN solver introduced in Ref. 48.

2. Optimization of the PNOs

The real-space PNOs are optimized in a manner similar to how the Hartree–Fock/Kohn–Sham orbitals are optimized by using the BSH Green’s operator in an iterative fashion.^36,44 However, compared to the mean-field methods, the number of PNOs is not fixed; this makes the PNO optimization more complicated compared to the standard LCAO-based PNO methods in which a set of PNOs is determined once and fixed. The current procedure is closer in spirit to the optimized PNO-CCSD approach of Ref. 45, where the number of PNOs was varied.

Before optimization, a given set of PNOs is truncated by diagonalizing the reduced density matrix given in Eq. (17) and using the eigenvalues as a truncation criterium, i.e., for each pair, the PNOs that correspond to eigenvalues of the reduced density matrix smaller than δ are removed. This will be referred to as compression of the PNOs; this step is completely analogous how the PNOs are determined in the LCAO-based MP2 and CC methods. In the representation that diagonalizes the reduced density matrix, the PNO gradient of Eq. (15) is

\begin{align} |\frac{δ L}{δ a_{i j}}⟩ = & 2 \sum_{c_{i j}} ({\tilde{t}}_{a_{i j} c_{i j}} g_{j}^{c_{i j}} | i (⟩ + \sum_{c_{i j}} {\tilde{t}}_{c_{i j} a_{i j}} g_{i}^{c_{i j}} | j (⟩) + \hat{F} | a_{i j} (⟩ d_{a_{i j}} + \sum_{c_{i j}} | c_{i j} (⟩ E_{c_{i j} a_{i j}} \\ - \sum_{k} (\sum_{e_{i k}} | e_{i k} (⟩ F_{k j} D_{e_{i k} a_{i j}} + \sum_{e_{k j}} | e_{k j} (⟩ F_{i k} D_{e_{k j} a_{i j}}), \end{align}

(38)

where $d_{a_{i j}}$ is the corresponding eigenvalue of the PNO RDM (note that the corresponding inter-pair RDMs are not diagonal). The gradient can be rearranged to give

(\hat{T} + \frac{E_{a_{i j} a_{i j}}}{d_{a_{i j}}} - F_{i i} - F_{i j}) | a_{i j} (⟩ = - \frac{2}{d_{a_{i j}}} (({\hat{v}}_{n e} + 2 \hat{J} - \hat{K}) | a_{i j} (⟩ + V_{a_{i j}})

(39)

with the potential

\begin{align} V_{a_{i j}} \equiv & 2 \sum_{c_{i j}} ({\tilde{t}}_{a_{i j} c_{i j}} g_{j}^{c_{i j}} | i (⟩ + {\tilde{t}}_{c_{i j} a_{i j}} g_{i}^{c_{i j}} | j (⟩) + \sum_{c_{i j} \neq a_{i j}} | c_{i j} (⟩ E_{c_{i j} a_{i j}} \\ - \sum_{k \neq j} \sum_{e_{i k}} | e_{i k} (⟩ F_{k j} D_{e_{i k} a_{i j}} - \sum_{k \neq i} \sum_{e_{k j}} | e_{k j} (⟩ F_{i k} D_{e_{k j} a_{i j}} . \end{align}

(40)

The optimized PNOs are then given as a solution to

| a_{i j} (⟩ = - \frac{2}{d_{a_{i j}}} {\hat{G}}_{ϵ_{a_{i j}}} (({\hat{v}}_{n e} + 2 \hat{J} - \hat{K}) | a_{i j} (⟩ + V),

(41)

where ${\hat{G}}_{ϵ_{a_{i j}}}$ is the pair-specific 3D BSH Green’s operator parameterized by

ϵ_{a_{i j} a_{i j}} \equiv \frac{E_{a_{i j} a_{i j}}}{d_{a_{i j}}} - F_{i i} - F_{j j} .

(42)

This operator is constructed and applied as described in Ref. 36. Due to the $\frac{2}{d_{a_{i j}}}$ factor in Eq. (41), the MRA precision ϵ should be in a similar range as the RDM cutoff threshold δ if all PNOs should be optimized to convergence. Since the RDM threshold is usually chosen to be 10⁻⁶ or smaller, this is, in general, not practical, but since the PNOs corresponding to the small RDM eigenvalues will also contribute less to the overall energy, it seems unnecessary to fully optimize those PNOs.

The PNO optimization as described below does not increase their number; the PNO update keeps the number of PNOs unchanged, whereas the compression step can only decrease their number. Thus, it is important to consider how to adaptively grow the set of PNOs to be able to guarantee the target precision of the MP2 energy. The simplest possibility is just to read in a sufficiently large set of guess PNOs, defined by some predefined atomic basis set and use this as first guess for all pairs. Here, we instead adaptively generate the guess PNOs from scratch in a pair-specific way as follows: The procedure is similar to how guess functions were generated in MRA-CIS, namely, the orbital-like functions in the linear-response equations were created by multiplying “excitation operators” onto the occupied reference orbitals.⁴⁹ In Ref. 49, monomials (i.e., x^ay^bz^c) were used as “excitation operators” in order to introduce desired symmetries into the guess functions. Due to the regularization with the explicitly correlated ansatz, we expect that the low order of the monomials used to generate the guess PNOs will be sufficient for an accurate representation. In this work, we will follow this approach where we wrap the roots of the monomials with a sin function (i.e., $x^{a} y^{b} z^{c} \to \sin^{a} (x) \sin^{b} (y) \sin^{c} (z)$ ⁠) in order to remove the growth of the functions toward the boundary, which can lead to large undesired numerical noise at the boundaries after multiplication with the reference orbital (see Table I). Reference 50 also used a similar strategy to create natural orbital-like functions from the occupied reference orbitals called “product plane waves.” The main difference is that the functions here only serve as guess to provide correct initial symmetries and are optimized afterward.

TABLE I.

Multipole operators to create the initial PNOs of each macroiteration and the used abbreviations.

dipole+ = {x, y, z, r²}

quadrupole = {x, y, z, x², y², z², xy, xz, yz}

octupole = quadrupole ∪ {x³, y³, z³, x²y, x²z, xy², xz², y²z, yz², xyz}

multipole = \{\begin{matrix} macroiteration = 1 : dipole + \\ macroiteration = 2 : quadrupole \\ macroiteration > 2 : octopole \end{matrix})

The guess functions for the PNOs are created pair-specifically as

\{| {a_{i j}}^{(0)} (⟩\} = \{e x o p | i (⟩ \cup e x o p | j (⟩\},

(43)

where exop is a set of PNO generators given in Table I. After initialization, the PNOs are iterated until convergence, which defines one macroiteration. A macroiteration consist of several microiterations where the PNOs are compressed and optimized. Due to compression, the number of PNOs is usually reduced. In the next macroiteration, a new set of PNOs is added to the already converged ones and the procedure is repeated until the overall convergence of the energy is reached. In order to save computational time, it is possible to grow the PNO ranks with a lower MRA threshold and reiterate the final pre-optimized PNOs again with a larger MRA threshold. The two phases of PNO optimization are illustrated in Figs. 1 and 2.

FIG. 1.

View large Download slide

Macroiteration cycles of the adaptive MRA-PNO-MP2-F12 solver.

FIG. 2.

View large Download slide

Microiteration cycles of the adaptive MRA-PNO-MP2-F12 solver.

III. COMPUTATIONAL DETAILS

The MRA-PNO-MP2-F12 method was implemented in the open-source package MADNESS (revision c6bdfa1fd762080ce9754d17c0b88f26b366b529).⁴⁰

MRA-PNO-MP2-F12 computations with MRA target precision ϵ and PNO truncation threshold δ are abbreviated as MRA(ϵ, δ) (Table II). In this work, we restrict ourselves to 7 initial polynomials per dimension (i.e., k = 7) for 7³ = 343 grid points per volume element.

TABLE II.

Overview over notation and conventions.

Projector onto the occupied space

\hat{O} = \sum_{i \in occ} | i (⟩ ⟨) i |

Projector onto the virtual space

\hat{Q} = 1 - \hat{O}

Two-body notation

{\hat{X}}_{1} \equiv \hat{X} (1), {\hat{X}}_{12} \equiv \hat{X} (1,2)

Two-body Fock operator (same for

\hat{K}

and

\hat{J}

⁠)

{\hat{F}}_{12} = {\hat{F}}_{1} + {\hat{F}}_{2}

Fock matrix

F_{p q} = ⟨) p | \hat{F} | q (⟩

Fock operator

\hat{F} = \hat{T} + \hat{J} - \hat{K}

Potential part

{\hat{v}}^{F} \equiv \hat{F} - \hat{T}

Particle and orbital permutation

{\hat{P}}_{12}, {\hat{P}}_{p q},

Explicit correlation factor

f_{12} \equiv f (1,2) = 1 - \exp (- r_{12})

Coulomb potential

g_{12} = r_{12}^{- 1}

Hylleraas functional

L^{(2)}

see Eqs. , , and

BSH Green’s function

{\hat{G}}_{μ}

see Ref. 36

MP2 pair function

|τ_ij⟩

MRA with accuracy 10^−x and PNO threshold of 10^−δ

MRA(x, δ)

aug-cc-pVXZ basis-set

aXZ

In general, a lower k value reduces memory demand in any vertex of the MRA tree but also leads to deeper trees and therefore more vertices.^41,51,52 It is, therefore, hard to state which k values should be used in general. A rule of thumb is to choose k ≥ 5 and increase it with target precision where k = 5 is usually only used in six dimensional representations where higher values have significant impact on memory consumption.¹⁴ We chose k = 7 because it usually performs well within the heuristics of madness for three dimensional representations and for the accuracy range chosen for this work.⁴⁰ We performed some tests with k = 8 and k = 6 and did not note significant changes.

Unless otherwise mentioned, the MRA computations are performed with reference orbitals localized via the standard Foster–Boys⁵³ procedure. All final correlation energies were obtained by the adaptive PNO generation algorithm as described in Secs. II D 2 and IV B.

For comparison with the LCAO-based methods, we used the MP2-F12 method with ansatz 2B, F+K approximation, and fixed-amplitudes (sp) ansatz⁵⁴ implemented in the ricc2⁵⁵ program of the TURBOMOLE⁵⁶ package with the basis sets of Dunning et al.^57–60 with the matching auxiliary (density-fitting) basis sets^61,62 of the TURBOMOLE library. In most plots and figures, the Dunning basis sets aug-cc-pVXZ are abbreviated as aVXZ. For second row elements, we used the corrected aug-cc-pV(X+d)Z⁶³ basis sets abbreviated as aV(Z+d)Z and the corresponding auxiliary basis sets.⁶⁴ All computations are performed using the frozen-core approximation (i.e., neglecting the core orbitals in the electron correlation treatment).

IV. RESULTS

A. Saturating the PNO subspace: The helium atom

In Fig. 3, the convergence of MRA-PNO-MP2-F12 computations with different sets of virtual orbitals and MRA thresholds is shown. The smallest set of initial PNOs is created by the “dipole+” operator of Table I, which for the case of the ground state He atom generates an s-type and three p-type orbitals from the reference orbital (crosses and open circles). The energy after 5 iterations (at point M1) differs by several mE_h from the MP2-F12/aV6Z reference, indicating an incomplete basis. Adding d-type functions (“quadrupole,” filled circles) decreases the basis set error. Note that since the energy is computed by projection, the variational principle does not apply, and the energy can be overestimated, resulting in positive errors with respect to the reference.

FIG. 3.

View large Download slide

Variation of the MRA-PNO-MP2-F12 energy of a helium atom with respect to the PNO generator set. The horizontal axis refers to the iterations of the two-phase PNO optimization algorithm (see Sec. II D 2). The most important PNOs at the end of each macroiteration are shown in the lower half of the figure.

In the second macroiteration (converged at point M2), additional PNOs of the same symmetry are generated and the energy changes significantly, while in the last macroiteration, no such changes are observed. At that point, the partial wave sets are saturated, and it can be seen that the 2s2p set (open circles) still misses the 2d functions (filled circles), but no f-functions are needed (filled circles and triangles are on top of each other).

The final macroiteration adds another set of functions of the respective symmetry, which does not change the energies much. However, the two sp-sets start to differ from the third iteration on, where the additional s- and p-functions are removed from the set in the low-threshold calculation (crosses), while they are kept in the tight-threshold calculation (open circles).

In order to reach the sub-millihartree range, it is necessary to not only have the partial wave expansion saturated but also increase the MRA and PNO thresholds. If the d- and f-functions are used, the difference to the MP2-F12/aV6Z reference is below 0.1 mE_h.

B. Saturating the PNO subspace: The methane molecule

A similar analysis was performed for methane, a highly symmetric system with more than one electron pair. Figure 4 illustrates the convergence of the MRA-PNO-MP2-F12 energy of a methane molecule with respect to the PNO generators used to generate the initial PNOs.

FIG. 4.

View large Download slide

Variation of the MRA-PNO-MP2-F12 energy of a methane molecule with respect to the PNO generator set. The horizontal axis refers to the iterations of the two-phase PNO optimization algorithm (see Sec. II D 2).

Just as in the He atom, the error with respect to the MP2-F12/aV6Z reference decreases as the number of PNO generators is increased. The choice of the MRA threshold is investigated in Fig. 5 where the same computation is performed with different MRA thresholds. As in the case of a helium atom, the 10⁻⁴ and 10⁻⁵ thresholds do not lead to a significantly different result if accuracy in the millihartree regime is desired, while an MRA threshold of 10⁻³ underestimates the correlation energy by several mE_h. If the threshold is, however, tightened to 10⁻⁴ in the last macrocycle, the result is nearly the same as for the computation with MRA(4,6). This indicates that the MRA threshold does not have significant effects on the quality of the generated PNOs and the error arises dominantly from the accumulation of small errors in the individual pair energies.

FIG. 5.

View large Download slide

Variation of the MRA-PNO-MP2-F12 energy of a methane molecule with respect to the MRA threshold. In the last iteration of MRA(3-4,6), the threshold was increased from 10⁻³ to 10⁻⁴.

The first step in Fig. 5 of the third macrocycle (first filled blue dot of M3) of the MRA(3-4,6) computation shows the result for the PNOs generated and optimized with MRA(3,6) but the projected energy calculated with MRA(4,6). This steps corrects the majority of the error of the lower threshold MRA(3,4) result compared to the MP2-F12/aV6Z reference. Since there are still significant changes after the optimization with MRA(4,6), it is not enough to just recompute the energy with tightened precision; the orbitals need to be also relaxed. It is, however, in this case, sufficient to generate and pre-optimize the PNOs with the lower MRA precision so that only a few iterations with the tightened precision are necessary for an accurate result.

In Fig. 6, an adaptive protocol is used where the PNOs are generated and pre-optimized with MRA(3,6) and where the PNO generator sets also increase in each macrocycle. After convergence, the MRA threshold is increased to (4,6) and the pre-optimized orbitals are re-iterated but not further expanded. The computation is performed for the methane and the butane molecule resulting in sub-millihartree agreement with the MP2-F12/aV6Z reference.

FIG. 6.

View large Download slide

The MRA-PNO-MP2-F12 energy of methane and butane molecules. In the last iteration of MRA(3-4,6), the threshold was increased from 10⁻³ to 10⁻⁴. The horizontal axis refers to the iterations of the two-phase PNO optimization algorithm (see Sec. II D 2).

C. Energy computation: The projection formula vs the Hylleraas functional

In Fig. 7, the Hylleraas functional and the projected energies are compared for the water molecule where the $L^{(2)} [w_{i j}]$ part of the functional was computed using a quintuple-ζ auxiliary basis set. We note that the projected energies typically converge from below (i.e., the correlation energies are too large in magnitude), while the energies computed with the Hylleraas functions converge from above, as they should according to the variational principle. Energies evaluated with the two approaches converge to the same value but differ significantly in the beginning. The advantage of the projected energies is that there is no need for a 6D computation or the introduction of an auxiliary basis for the computation of the $f_{12} [{\hat{K}}_{12}, f_{12}]$ term.

FIG. 7.

View large Download slide

Comparison of the projected MRA-PNO-MP2-F12 energies with the value of the Hylleraas functional for water with localized orbitals. The PNOs are generated with the adaptive multipole procedure.

D. Reference orbitals: Canonical vs localized

The effects of the localization are illustrated in Fig. 8. The MRA calculation with canonical orbitals differs significantly from the MP2-F12/aV6Z reference and the error increases with the size of the linear alkane chain, while both are not the case for the localized computation.

FIG. 8.

View large Download slide

Difference in the total MRA-PNO-MP2-F12 correlation energy compared to MP2-F12/aV6Z for the linear alkanes (methane to pentane). Comparison of canonical and localized (Foster–Boys) orbitals. The insets show the average and maximal PNO ranks.

The more an orbital becomes delocalized, the smaller the absolute values of the function on any grid point is. Due to the truncation scheme of MRA, this leads to loss of precision. In the extreme case, when an orbital is delocalized enough and its function values fall below the MRA threshold, it would be truncated completely. Tightening the thresholds helps pushing this limiting case further, but loss of precision will occur eventually. Localized orbitals, on the other hand, have magnitude that is asymptotically independent of the system size, and this allows us to maintain high precision.

E. Error growth with the system size: n-alkanes

We computed MRA-PNO-MP2-F12 energies for n-alkanes of varying sizes (methane to decane), comparing the results to MP2-F12/aVXZ (X = D, T, Q, 5, 6) in Figs. 9 and 10. Note that the conventional (LCAO) quintuple- and hextuple-ζ calculations could not be performed for the larger alkanes, and the quadruple-ζ MP2-F12 is not a sufficiently precise reference (e.g., the MP2-F12/aVQZ energy differs from the MP2-F12/aV6Z reference by more than a mE_h already for pentane). Hence, in Fig. 10, the MRA values were used as the reference.

FIG. 9.

View large Download slide

Errors in the MRA-PNO-MP2-F12 and LCAO MP2-F12 energies relative to the MP2-F12/aV6Z energies for linear n-alkanes (methane → pentane). The inset shows the mE_h regime in greater detail.

FIG. 10.

View large Download slide

Errors in the LCAO MP2-F12 energies relative to the MRA(3-4,6)-multipole energies for linear n-alkanes (methane → decane). On the right, the average pair energy and the maximum and average PNO ranks are plotted.

For the smaller basis sets, growth of the absolute error with respect to the system size is observable, as expected due to the accumulation of errors in the individual pair energies. The dependence of the error on the system size is strongly nonlinear since these small systems are too small to reach the asymptotic $O (N)$ behavior of the MP2 energy (see also the inset in Fig. 10).

The MRA results are in good agreement with the corresponding best available MP2-F12/aVXZ energies. The error of MRA-PNO-MP2-F12 energy does not show artifactual dependence on the system size as was observed for the standard MRA(3,6)-MP2-F12 in Ref. 18. We believe that this is due to the use of localized reference orbitals, which prevents the number of significant pairs from growing quadratically (see the inset in Fig. 10), and due to the increased MRA precision threshold. Note that, in Fig. 6, this system size dependent error is observable for the computations with the MRA threshold of 10⁻³.

F. Error growth with the system size: n-silanes

A similar analysis was performed for n-silanes to test the ability to describe heavier-element core orbitals in the MRA-PNO-MP2-F12 method, which, despite the frozen core approximation, enter the expressions via the strong orthogonality projectors. The results presented in Fig. 11 suggest that, for heavier elements, tighter MRA thresholds are needed, presumably to be able to ensure strict orthogonality of the unoccupied states to the core one-particle states; the unoccupied states need to be resolved on the lengthscale of the tightest occupied orbitals to maintain strict orthogonality. The use of pseudopotentials could alleviate the need for the high resolution of the core regions of the unoccupied orbitals. In all-electron computations, no such simplifications appear possible.

FIG. 11.

View large Download slide

Difference in the total MRA-PNO-MP2-F12 and MP2-F12 correlation energies compared to MP2-F12/aug-cc-pw(5+d)Z for the first four linear silanes. (Top) MRA-PNO-MP2-F12 and (bottom) canonical LCAO MP2-F12.

G. Conventional LCAO PNOs vs optimized MRA PNOs: Comparison of ranks and precision

MRA PNOs can be viewed as essentially exact PNOs (modulo the finite MRA precision parameter). In principle, these should be equivalent to the PNOs obtained from LCAO MP1-F12 amplitudes [with the strong orthogonality projector of Eq. (2) rather than the traditional F12 projector] determined in a sufficiently complete basis. However, in the LCAO representation, the PNOs are traditionally formed from conventional (non-F12) MP1 amplitudes. Thus, a fair comparison of our approach is to the F12-optimized PNOs of Tew et al.,³⁰ which result in more compact PNO^30,66 sets that approximate the exact PNOs determined from MP1-F12. To understand the differences between the exact and approximate LCAO (F12-)PNOs, we compared PNO ranks and the corresponding PNO truncation errors in MRA and LCAO approaches for the HCl–H₂O complex (see Table III). For comparable accuracy of about −547mE_h, the LCAO requires average PNO rank 84, where the optimized PNOs of MRA require only 51 PNOs per pair. This suggests that the LCAO F12-PNOs are far from being optimal even for high ranks.

TABLE III.

PNO ranks for H₂O–HCl: Comparison between optimized MRA and standard LCAO orbitals.

Method	Max rank	Av. rank	Energy (mE_h)
PNO-MP2-F12, aV(Q+d)Z, t = 6	41	22	−542.025
PNO-MP2-F12, aV(Q+d)Z, t = 7	68	42	−543.863
PNO-MP2-F12, aV(Q+d)Z, t = 8	105	70	−544.365
MP2-F12, aV(Q+d)Z	307	307	−544.520
PNO-MP2-F12, aV(5+d)Z, t = 6	45	23	−542.961
PNO-MP2-F12, aV(5+d)Z, t = 7	83	46	−545.485
PNO-MP2-F12, aV(5+d)Z, t = 8	127	84	−546.138
MP2-F12, aV(5+d)Z	503	503	−546.350
MP2-F12, aV6Z	768	768	−547.056
MRA(3-4,6)	37	17	−548.729
MRA(3-4,7)	69	32	−547.849
MRA(3-4,8)	118	51	−547.058

Method	Max rank	Av. rank	Energy (mE_h)
PNO-MP2-F12, aV(Q+d)Z, t = 6	41	22	−542.025
PNO-MP2-F12, aV(Q+d)Z, t = 7	68	42	−543.863
PNO-MP2-F12, aV(Q+d)Z, t = 8	105	70	−544.365
MP2-F12, aV(Q+d)Z	307	307	−544.520
PNO-MP2-F12, aV(5+d)Z, t = 6	45	23	−542.961
PNO-MP2-F12, aV(5+d)Z, t = 7	83	46	−545.485
PNO-MP2-F12, aV(5+d)Z, t = 8	127	84	−546.138
MP2-F12, aV(5+d)Z	503	503	−546.350
MP2-F12, aV6Z	768	768	−547.056
MRA(3-4,6)	37	17	−548.729
MRA(3-4,7)	69	32	−547.849
MRA(3-4,8)	118	51	−547.058

H. Performance for relative energies: Water clusters

The LCAO-based computations often achieve remarkable accuracy for relative energies compared to their performance on total energies due to robust cancellation of “intra-atomic” errors as the geometry varies; the LCAO-based PNO relative energies are also known to be very precise. The MRA-based approaches, in contrast, do not generally profit from the intra-atomic error cancellation by construction, since the structure of the basis near each nucleus is not fixed. Thus, it is vitally important to examine the performance of MRA-PNO-MP2-F12 for relative energetics. To this end, we computed binding energies of small water clusters up to tetramers. The MP2 binding energies of the water clusters are just a few millihartrees; thus, generally very high precision in absolute energies, even for the LCAO approaches, is needed to compute the binding energies precisely.

The errors due to the MP2-F12/aV6Z reference for the total and binding energies are shown in Figs. 12 and 13, respectively. All MRA-PNO-MP2-F12 results give precise results in the binding energies, except MRA(3-4,6); we hypothesized that this is caused by premature truncation of PNOs in the optimization. To investigate this further, we performed computations with a tight PNO truncation threshold (δ = 10⁻⁸), but the final PNO ranks were restricted to a fixed small value [see the “MRA(3-4,8)-maxrank-X” data in Figs. 12 and 13, with X = 20, 30]. The absolute energies obtained with the final PNO ranks were restricted to 20 and 30 PNOs, which differed significantly from those with the non-truncated counterpart (Fig. 12); the average and maximum PNO ranks for the single H₂O molecule in the latter were 60 and 90, i.e., much larger. In contrast, the effect of the final PNO truncation on the relative energies (Fig. 13) is very small due to the apparent systematic convergence of the PNO truncation error with the cluster size. This suggests that the PNO construction procedure can be improved further by gradually reducing PNO ranks during their construction.

FIG. 12.

View large Download slide

Total MRA-PNO-MP2-F12 and MP2-F12 correlation energies of small water clusters compared to the CBS MP2 (estimated by MP2-F12/aV6Z). (Top) MRA-PNO-MP2-F12 and (bottom) canonical LCAO MP2-F12.

FIG. 13.

View large Download slide

MRA-PNO-MP2-F12 and MP2-F12 correlation part of the binding energies of small water clusters compared to MP2-F12/aV6Z. The correlation part of the binding energies is computed as $n \cdot E^{MP 2 - F 12} (H_{2} O) - E^{MP 2 - F 12} ({(H_{2} O)}_{n})$ ⁠, and the total MP2-F12/aV6Z binding energies are 1.7, 2.4, and 5.7 mE_h, respectively. (Top) MRA-PNO-MP2-F12 and (bottom) canonical LCAO MP2-F12.

I. Formal scaling and expected computational performance

The formal computational complexity of the presented MRA-PNO-MP2-F12 approach is identical to that of the corresponding LCAO counterpart. The cost of computing the MP1 wave function in the LCAO representation grows with the fifth power of the system size or more precisely as $N_{occ}^{3} N_{virt}^{2}$ without any approximation. With localized orbitals, the cost reduces to $N_{occ}^{3}$ ⁠; using the real-space representation for the pair functions with localized occupied orbitals results in the same complexity.¹⁸ The complexity is further reduced to quadratic due to the exponential decay of the Fock operator with a finite HOMO–LUMO gap. Thus, the MRA-PNO-MP2-F12 method as implemented here has quadratic complexity with the system size. Further reduction of complexity to linear could be achieved by pair screening but was not implemented in our work. Although the complexity of our implementation is formally quadratic, the cost is dominated by the short-range pairs, which have the highest PNO ranks, and thus, the apparent scaling of the cost with the system size is subquadratic.

A rigorous and fair comparison of the computational cost of MRA vs LCAO PNO-MP2-F12 is complicated due to the need to compare costs of computations that achieve same precision. An additional practical complication was that we performed the two on different computers due to technical and licensing reasons. After accounting for the hardware factors, some MRA-PNO calculations were faster than the LCAO counterparts, namely, aug-cc-pV6Z PNO-MP2-F12 computations for propane and larger systems. For the LCAO calculations, the wall-times are dominated by the Hartree–Fock part, while for PNO-MRA-MP2 calculations, the wall-times are dominated by the MP2 part. In MRA based computations, Hartree–Fock solves for the optimized occupied orbitals only, and the most time consuming step is the actual optimization of the PNOs, while in LCAO based calculations, Hartree–Fock is solved by rotating in the full orbital space defined by the basis-set. In other words, in both calculations, the construction of the basis in which the MP2 equations are solved is the expensive part.⁶⁵

A potentially interesting question in comparing LCAO with MRA-PNO is also the performance of truncated schemes with regard to relative energies. We illustrated that truncated MRA-PNO schemes can, in principle, achieve consistent and accurate results in cases where small and medium LCAO based calculations fail even with the F12 correction (see Fig. 13). A rigorous investigation is beyond the scope of this work.

V. SUMMARY AND PERSPECTIVE

We presented a real-space implementation of the simplest many-body electronic structure method, namely, the second-order Moller–Plesset perturbation theory, with the pair functions represented as a linear combination of products of pair-natural orbitals represented in real space on an adaptive multiresolution spectral-element basis. Similarly to the previously demonstrated6-dimensional” MP2 approach^13,14 in which pair functions were represented directly in real space, regularization of the electron–electron Coulomb singularity is mandatory to be able to solve equations robustly and reach high precision. However, in contrast to the original approach in which each 6D pair function was written as a sum of single-particle products only locally (within each element of the adaptive volume decomposition), employing global separation of pair functions into products of orbitals allows us to greatly extend the size of systems that can be treated with highly precise real-space MP2 computations on up to 32 atoms demonstrated for the first time.

Compared to the PNO methods in the LCAO representation in which the PNOs are determined once and fixed (for an exception, see Ref. 45 in which PNOs were optimized), in our approach, the PNOs must be determined directly by iterative minimization of an appropriate Lagrangian (namely, the Hylleraas functional for the MP2 energy). Similarly, the rank of the PNO basis is varied adaptively to reach a given precision by employing PNO generating operators (which in atomic systems increase principal and orbital quanta of the argument functions) to repeatedly augment PNO subspaces for each pair. The PNO subspaces generated by this procedure converge rapidly and produce more compact and precise sets than those obtained in the LCAO representation.

Benchmark calculations on n-alkanes (up to decane) shows that a size-independent maximum rank of only 25–30 is needed to achieve chemical accuracy of the absolute correlation energy for the whole molecule with perhaps a few dozen atoms and can be applied to systems too large to be treatable by conventional LCAO MP2-F12 approaches of similar precision. Further calculations on a linear chain of alkanes show that the basis set limit can be reached in a systematic way, where, in this case, higher thresholds were necessary, presumably due to the increased number of core orbitals still present in the projector. Initial assessment for relative energies was performed using weakly interacting water clusters up to four molecules where the PNO representation adaptively grown to a low fixed rank achieved precise interaction energies despite having large absolute errors.

The proposed approach can be generalized straightforwardly to higher-order methods, such as coupled-cluster. Of particular interest would be application to systems in electronically excited states and to molecular properties other than energy, where the weaknesses of the LCAO representation become more dramatic. Another prospective application would be to atoms and molecules in strong electromagnetic fields where the atom-centered representations such as LCAO are not appropriate due to their inability to describe large spatial distortions of the electronic wavepacket.

ACKNOWLEDGMENTS

J.S.K. acknowledges the support from Fonds der Chemischen Industrie and Deutsche Forschungsgemeinschaft (Grant No. BI-1432/2-1). The work of F.A.B. was supported by the Deutsche Forschungsgemeinschaft (Grant No. BI-1432/2-1). The work of E.F.V. was supported by the U.S. National Science Foundation (Award Nos. 1550456 and 1800348).

REFERENCES

1.

V.

Blum

,

R.

Gehrke

,

F.

Hanke

,

P.

Havu

,

V.

Havu

,

X.

Ren

,

K.

Reuter

, and

M.

Scheffler

, “

Ab initio molecular simulations with numeric atom-centered orbitals

,”

Comput. Phys. Commun.

180

,

2175

–

2196

(

2009

).

https://doi.org/10.1016/j.cpc.2009.06.022

Google Scholar

Crossref

2.

T.

Schäfer

,

B.

Ramberger

, and

G.

Kresse

, “

Quartic scaling MP2 for solids: A highly parallelized algorithm in the plane wave basis

,”

J. Chem. Phys.

146

,

104101

(

2017

).

https://doi.org/10.1063/1.4976937

Google Scholar

Crossref

PubMed

3.

T.

Helgaker

,

W.

Klopper

, and

H.

Koch

, “

Basis-set convergence of correlated calculations on water

,”

J. Chem. Phys.

106

,

9639

–

9646

(

1997

).

https://doi.org/10.1063/1.473863

Google Scholar

Crossref

4.

J. J.

Shepherd

,

A.

Grüneis

,

G. H.

Booth

,

G.

Kresse

, and

A.

Alavi

, “

Convergence of many-body wave-function expansions using a plane-wave basis: From homogeneous electron gas to solid state systems

,”

Phys. Rev. B

86

,

035111

(

2012

).

https://doi.org/10.1103/physrevb.86.035111

Google Scholar

Crossref

5.

L.

Kong

,

F. A.

Bischoff

, and

E. F.

Valeev

, “

Explicitly correlated R12/F12 methods for electronic structure

,”

Chem. Rev.

112

,

75

–

107

(

2012

).

https://doi.org/10.1021/cr200204r

Google Scholar

Crossref

PubMed

6.

A.

Grüneis

,

J. J.

Shepherd

,

A.

Alavi

,

D. P.

Tew

, and

G. H.

Booth

, “

Explicitly correlated plane waves: Accelerating convergence in periodic wavefunction expansions

,”

J. Chem. Phys.

139

,

084112

(

2013

).

https://doi.org/10.1063/1.4818753

Google Scholar

Crossref

PubMed

7.

A.

Tajti

,

P. G.

Szalay

,

A. G.

Császár

,

M.

Kállay

,

J.

Gauss

,

E. F.

Valeev

,

B. A.

Flowers

,

J.

Vázquez

, and

J. F.

Stanton

, “

HEAT: High accuracy extrapolated ab initio thermochemistry

,”

J. Chem. Phys.

121

,

11599

–

11613

(

2004

).

https://doi.org/10.1063/1.1811608

Google Scholar

Crossref

PubMed

8.

J. R.

Flores

, “

High-precision atomic computations from finite-element techniques: Second-order correlation energies of rare-gas atoms

,”

J. Chem. Phys.

98

,

5642

–

5647

(

1993

).

https://doi.org/10.1063/1.464908

Google Scholar

Crossref

9.

H. M.

Schmidt

and

J.

Linderberg

, “

Finite-element computation of perturbation energies for the two-electron atom

,”

Phys. Rev. A

49

,

4406

–

4410

(

1994

).

https://doi.org/10.1103/physreva.49.4404

Google Scholar

Crossref

10.

J.

Ackermann

, “

Finite-element-method expectation values for correlated two-electron wave functions

,”

Phys. Rev. A

52

,

1968

–

1975

(

1995

).

https://doi.org/10.1103/physreva.52.1968

Google Scholar

Crossref

PubMed

11.

J. R.

Flores

and

D.

Kolb

, “

Atomic MP2 correlation energies fast and accurately calculated by FEM extrapolations

,”

J. Phys. B: At., Mol. Opt. Phys.

32

,

779

–

790

(

1999

).

https://doi.org/10.1088/0953-4075/32/3/019

Google Scholar

Crossref

12.

J. R.

Flores

, “

New benchmarks for the second-order correlation energies of Ne and Ar through the finite element MP2 method

,”

Int. J. Quantum Chem.

108

,

2172

–

2177

(

2008

).

https://doi.org/10.1002/qua.21742

Google Scholar

Crossref

13.

F. A.

Bischoff

,

R. J.

Harrison

, and

E. F.

Valeev

, “

Computing many-body wave functions with guaranteed precision: The first-order Møller-Plesset wave function for the ground state of helium atom

,”

J. Chem. Phys.

137

,

104103

(

2012

).

https://doi.org/10.1063/1.4747538

Google Scholar

Crossref

PubMed

14.

F. A.

Bischoff

and

E. F.

Valeev

, “

Computing molecular correlation energies with guaranteed precision

,”

J. Chem. Phys.

139

,

114106

(

2013

).

https://doi.org/10.1063/1.4820404

Google Scholar

Crossref

PubMed

15.

F. A.

Bischoff

, “

Regularizing the molecular potential in electronic structure calculations. I. SCF methods

,”

J. Chem. Phys.

141

,

184105

(

2014

).

https://doi.org/10.1063/1.4901021

Google Scholar

Crossref

PubMed

16.

F. A.

Bischoff

, “

Regularizing the molecular potential in electronic structure calculations. II. Many-body methods

,”

J. Chem. Phys.

141

,

184106

(

2014

).

https://doi.org/10.1063/1.4901022

Google Scholar

Crossref

PubMed

17.

J. S.

Kottmann

and

F. A.

Bischoff

, “

Coupled-cluster in real space. 2. cc2 excited states using multiresolution analysis

,”

J. Chem. Theor. Comput.

13

,

5956

–

5965

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00695

Google Scholar

Crossref

18.

J. S.

Kottmann

and

F. A.

Bischoff

, “

Coupled-cluster in real space. 1. CC2 ground state energies using multiresolution analysis

,”

J. Chem. Theor. Comput.

13

,

5945

–

5955

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00694

Google Scholar

Crossref

19.

E. G.

Hohenstein

,

R. M.

Parrish

,

C. D.

Sherrill

, and

T. J.

Martínez

, “

Communication: Tensor hypercontraction. III. Least-squares tensor hypercontraction for the determination of correlated wavefunctions

,”

J. Chem. Phys.

137

,

221101

(

2012

).

https://doi.org/10.1063/1.4768241

Google Scholar

Crossref

PubMed

20.

N.

Mardirossian

,

J. D.

McClain

, and

G. K.-L.

Chan

, “

Lowering of the complexity of quantum chemistry methods by choice of representation

,”

J. Chem. Phys.

148

,

044106

(

2018

).

https://doi.org/10.1063/1.5007779

Google Scholar

Crossref

PubMed

21.

C.

Edmiston

and

M.

Krauss

, “

Configuration-interaction calculation of H₃ and H₂

,”

J. Chem. Phys.

42

,

1119

–

1120

(

1965

).

https://doi.org/10.1063/1.1696050

Google Scholar

Crossref

22.

C.

Edmiston

and

M.

Krauss

, “

Pseudonatural orbitals as a basis for the superposition of configurations. II. Energy surface for linear H₃

,”

J. Comput. Phys.

49

,

192

–

205

(

1968

).

https://doi.org/10.1063/1.1669809

Google Scholar

Crossref

23.

W.

Meyer

, “

Ionization energies of water from PNO-CI calculations

,”

Int. J. Quantum Chem.

5

,

341

–

348

(

1971

).

https://doi.org/10.1002/qua.560050839

Google Scholar

Crossref

24.

W.

Meyer

, “

PNO-CI studies of electron correlation effects. I. Configuration expansion by means of nonorthogonal orbitals, and application to the ground state and ionized states of methane

,”

J. Chem. Phys.

58

,

1017

–

1035

(

1973

).

https://doi.org/10.1063/1.1679283

Google Scholar

Crossref

25.

R.

Ahlrichs

,

H.

Lischka

,

V.

Staemmler

, and

W.

Kutzelnigg

, “

PNO–CI (pair natural orbital configuration interaction) and CEPA–PNO (coupled electron pair approximation with pair natural orbitals) calculations of molecular systems. I. Outline of the method for closed-shell states

,”

J. Chem. Phys.

62

,

1225

–

1234

(

1975

).

https://doi.org/10.1063/1.430637

Google Scholar

Crossref

26.

F.

Neese

,

F.

Wennmohs

, and

A.

Hansen

, “

Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method

,”

J. Chem. Phys.

130

,

114108

(

2009

).

https://doi.org/10.1063/1.3086717

Google Scholar

Crossref

PubMed

27.

F.

Neese

,

A.

Hansen

, and

D. G.

Liakos

, “

Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis

,”

J. Chem. Phys.

131

,

064103

(

2009

).

https://doi.org/10.1063/1.3173827

Google Scholar

Crossref

PubMed

28.

P.

Pinski

,

C.

Riplinger

,

E. F.

Valeev

, and

F.

Neese

, “

Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. I. An efficient and simple linear scaling local MP2 method that uses an intermediate basis of pair natural orbitals

,”

J. Chem. Phys.

143

,

034108

(

2015

).

https://doi.org/10.1063/1.4926879

Google Scholar

Crossref

PubMed

29.

H.-J.

Werner

,

G.

Knizia

,

C.

Krause

,

M.

Schwilk

, and

M.

Dornbach

, “

Scalable electron correlation Methods I: PNO-LMP2 with linear scaling in the molecular size and near-inverse-linear scaling in the number of processors

,”

J. Chem. Theory Comput.

11

,

484

–

507

(

2015

).

https://doi.org/10.1021/ct500725e

Google Scholar

Crossref

PubMed

30.

D. P.

Tew

and

C.

Hättig

, “

Pair natural orbitals in explicitly correlated second-order Møller-Plesset theory

,”

Int. J. Quantum Chem.

113

,

224

–

229

(

2012

).

https://doi.org/10.1002/qua.24098

Google Scholar

Crossref

31.

B. K.

Alpert

, “

Sparse representation of smooth linear operators

,” Ph.D. thesis,

Yale University

,

1990

.

Google Scholar

32.

G.

Beylkin

,

R.

Coifman

, and

V.

Rokhlin

, “

Fast wavelet transforms and numerical algorithms I

,”

Commun. Pure Appl. Math.

44

,

141

–

183

(

1991

).

https://doi.org/10.1002/cpa.3160440202

Google Scholar

Crossref

33.

B. K.

Alpert

, “

A class of bases in L² for the sparse representation of integral operators

,”

SIAM J. Math. Anal.

24

,

246

–

262

(

1993

).

https://doi.org/10.1137/0524016

Google Scholar

Crossref

34.

B. K.

Alpert

,

G.

Beylkin

,

D.

Gines

, and

L.

Vozovoi

, “

Adaptive solution of partial differential equations in multiwavelet bases

,”

J. Comput. Phys.

182

,

149

–

190

(

2002

).

https://doi.org/10.1006/jcph.2002.7160

Google Scholar

Crossref

35.

G.

Beylkin

,

V.

Cheruvu

, and

F.

Perez

, “

Fast adaptive algorithms in the non-standard form for multidimensional problems

,”

Appl. Comput Harmonic Anal.

24

,

354

–

377

(

2008

).

https://doi.org/10.1016/j.acha.2007.08.001

Google Scholar

Crossref

36.

R. J.

Harrison

,

G. I.

Fann

,

T.

Yanai

,

Z.

Gan

, and

G.

Beylkin

, “

Multiresolution quantum chemistry: Basic theory and initial applications

,”

J. Chem. Phys.

121

,

11587

–

11598

(

2004

).

https://doi.org/10.1063/1.1791051

Google Scholar

Crossref

PubMed

37.

We follow the notation used in previous works denoting the particles onto which an operator acts as subscripts, i.e., ${\hat{O}}_{12} \equiv \hat{O} (1,2)$ ⁠. For one-particle equations, the subscripts are dropped for brevity.

38.

T.

Yanai

,

G. I.

Fann

,

Z.

Gan

,

R. J.

Harrison

, and

G.

Beylkin

, “

Multiresolution quantum chemistry in multiwavelet bases: Analytic derivatives for Hartree–Fock and density functional theory

,”

J. Chem. Phys.

121

,

2866

–

2876

(

2004

).

https://doi.org/10.1063/1.1768161

Google Scholar

Crossref

PubMed

39.

F. A.

Bischoff

, “

Analytic second nuclear derivatives of Hartree-Fock and DFT using multi-resolution analysis

,”

J. Chem. Phys.

146

,

124126

(

2017

).

https://doi.org/10.1063/1.4978957

Google Scholar

Crossref

PubMed

40.

R. J.

Harrison

,

G.

Beylkin

,

F. A.

Bischoff

,

J. A.

Calvin

,

G. I.

Fann

,

J.

Fosso-Tande

,

D.

Galindo

,

J. R.

Hammond

,

R.

Hartman-Baker

,

J. C.

Hill

,

J.

Jia

,

J. S.

Kottmann

,

M.-J.

Yvonne Ou

,

J.

Pei

,

L. E.

Ratcliff

,

M. G.

Reuter

,

A. C.

Richie-Halford

,

N. A.

Romero

,

H.

Sekino

,

W. A.

Shelton

,

B. E.

Sundahl

,

W. S.

Thornton

,

E. F.

Valeev

,

A.

Vázquez-Mayagoitia

,

N.

Vence

,

T.

Yanai

, and

Y.

Yokoi

, “

Madness: A multiresolution, adaptive numerical environment for scientific simulation

,”

SIAM J. Sci. Comput.

38

,

S123

–

S142

(

2016

).

https://doi.org/10.1137/15m1026171

Google Scholar

Crossref

41.

F. A.

Bischoff

, “

Computing accurate molecular properties in real space using multiresolution analysis

,” in

Advances in Quantum Chemistry

(

Elsevier

,

2019

).

Google Scholar

Crossref

42.

E. A.

Hylleraas

, “

Über den grundterm der zweielektronenprobleme von H⁻, He, Li⁺, Be⁺⁺ usw

,”

Z. Phys.

65

,

209

–

225

(

1930

).

https://doi.org/10.1007/bf01397032

Google Scholar

Crossref

43.

In this work, we only focus on the electron–electron cusp. The electron–nuclear cusp also plays a role in the real-space many-body methods since its seam has the same dimension as that of the electron–electron cusp. Although it is possible to regularize the electron–nuclei cusps via the NEMO approach described in Refs. 15 and 16, this is less crucial for the PNO-based approach we pursue here. It is because the cusps in PNOs are pointwise, and just like the cusps in the Hartree–Fock orbitals, these can be effectively described by the adaptive multiresolution refinement of the basis.

44.

M.

Kalos

, “

Monte Carlo calculations of the ground state of three- and four-body nuclei

,”

Phys. Rev.

128

,

1791

–

1795

(

1962

).

https://doi.org/10.1103/physrev.128.1791

Google Scholar

Crossref

45.

M. C.

Clement

,

J.

Zhang

,

C. A.

Lewis

,

C.

Yang

, and

E. F.

Valeev

, “

Optimized pair natural orbitals for the coupled cluster methods

,”

J. Chem. Theory Comput.

14

,

4581

–

4589

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00294

Google Scholar

Crossref

PubMed

46.

W.

Kutzelnigg

and

J. D.

Morgan

, “

Rates of convergence of the partial-wave expansions of atomic correlation energies

,”

J. Chem. Phys.

96

,

4484

(

1992

).

https://doi.org/10.1063/1.462811

Google Scholar

Crossref

47.

W.

Kutzelnigg

, “

r₁₂-Dependent terms in the wave function as closed sums of partial wave amplitudes for large l

,”

Theor. Chim. Acta

68

,

445

–

469

(

1985

).

https://doi.org/10.1007/bf00527669

Google Scholar

Crossref

48.

R. J.

Harrison

, “

Krylov subspace accelerated inexact Newton method for linear and nonlinear equations

,”

J. Comput. Chem.

25

,

328

–

334

(

2004

).

https://doi.org/10.1002/jcc.10108

Google Scholar

Crossref

PubMed

49.

J. S.

Kottmann

,

S.

Höfener

, and

F. A.

Bischoff

, “

Numerically accurate linear response-properties in the configuration-interaction singles (CIS) approximation

,”

Phys. Chem. Chem. Phys.

17

,

31453

–

31462

(

2015

).

https://doi.org/10.1039/c5cp00345h

Google Scholar

Crossref

PubMed

50.

T. E.

Baker

,

K.

Burke

, and

S. R.

White

, “

Accurate correlation energies in one-dimensional systems from small system-adapted basis functions

,”

Phys. Rev. B

97

,

085139

(

2018

).

https://doi.org/10.1103/physrevb.97.085139

Google Scholar

Crossref

51.

J. S.

Kottmann

, “

Coupled-cluster in real space

,” Ph.D. thesis,

Humboldt-Universität zu Berlin; Mathematisch-Naturwissenschaftliche Fakultät

,

2018

.

Google Scholar

52.

S. R.

Jensen

, “

Real-space all-electron density functional theory with multiwavelets

,” Ph.D. thesis,

UiT The Arctic University of Norway

,

2014

.

Google Scholar

53.

J.

Foster

and

S.

Boys

, “

Canonical configurational interaction procedure

,”

Rev. Mod. Phys.

32

,

300

(

1960

).

https://doi.org/10.1103/revmodphys.32.300

Google Scholar

Crossref

54.

R. A.

Bachorz

,

F. A.

Bischoff

,

A.

Glöß

,

C.

Hättig

,

S.

Höfener

,

W.

Klopper

, and

D. P.

Tew

, “

The MP2-F12 method in the TURBOMOLE program package

,”

J. Comput. Chem.

32

,

2492

–

2513

(

2011

).

https://doi.org/10.1002/jcc.21825

Google Scholar

Crossref

PubMed

55.

O.

Christiansen

,

H.

Koch

,

P.

Jørgensen

, and

J.

Olsen

, “

Excitation energies of H₂O, N₂ and C₂ in full configuration interaction and coupled cluster theory

,”

Chem. Phys. Lett.

256

,

185

–

194

(

1996

).

https://doi.org/10.1016/0009-2614(96)00394-6

Google Scholar

Crossref

56.

TURBOMOLE V7.0 2015, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com.

57.

T. H.

Dunning

, “

Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen

,”

J. Chem. Phys.

90

,

1007

–

1023

(

1989

).

https://doi.org/10.1063/1.456153

Google Scholar

Crossref

58.

A. K.

Wilson

,

T.

van Mourik

, and

T. H.

Dunning

, Jr., “

Gaussian basis sets for use in correlated molecular calculations. VI. Sextuple zeta correlation consistent basis sets for boron through neon

,”

J. Mol. Struct.: THEOCHEM

388

,

339

–

349

(

1996

).

https://doi.org/10.1016/s0166-1280(96)04689-1

Google Scholar

Crossref

59.

P. J.

Knowles

and

H.-J.

Werner

, “

An efficient method for the evaluation of coupling coefficients in configuration interaction calculations

,”

Chem. Phys. Lett.

145

,

514

–

522

(

1988

).

https://doi.org/10.1016/0009-2614(88)87412-8

Google Scholar

Crossref

60.

R. A.

Kendall

,

T. H.

Dunning

, and

R. J.

Harrison

, “

Electron affinities of the first-row atoms revisited. systematic basis sets and wave functions

,”

J. Chem. Phys.

96

,

6796

–

6806

(

1992

).

https://doi.org/10.1063/1.462569

Google Scholar

Crossref

61.

A.

Halkier

,

W.

Klopper

,

T.

Helgaker

, and

P.

Jørgensen

, “

Basis-set convergence of the molecular electric dipole moment

,”

J. Chem. Phys.

111

,

4424

–

4430

(

1999

).

https://doi.org/10.1063/1.480036

Google Scholar

Crossref

62.

C.

Hättig

, “

Optimization of auxiliary basis sets for RI-MP2 and RI-CC2 calculations: Core–valence and quintuple-ζ basis sets for H to Ar and QZVPP basis sets for Li to Kr

,”

Phys. Chem. Chem. Phys.

7

,

59

–

66

(

2005

).

https://doi.org/10.1039/b415208e

Google Scholar

Crossref

63.

T. H.

Dunning

,

K. A.

Peterson

, and

A. K.

Wilson

, “

Gaussian basis sets for use in correlated molecular calculations. X. The atoms aluminum through argon revisited

,”

J. Chem. Phys.

114

,

9244

–

9253

(

2001

).

https://doi.org/10.1063/1.1367373

Google Scholar

Crossref

64.

K. E.

Yousaf

and

K. A.

Peterson

, “

Optimized complementary auxiliary basis sets for explicitly correlated methods: aug-cc-pvnz orbital basis sets

,”

Chem. Phys. Lett.

476

,

303

–

307

(

2009

).

https://doi.org/10.1016/j.cplett.2009.06.003

Google Scholar

Crossref

65.

The number of iterations in the Hartree–Fock calculation did not change significantly.

66.

F.

Pavošević

,

F.

Neese

, and

E. F.

Valeev

, “

Geminal-spanning orbitals make explicitly correlated reduced-scaling coupled-cluster methods robust, yet simple

,”

J. Chem. Phys.

141

,

054106

(

2014

).

https://doi.org/10.1063/1.4890002

Google Scholar

Crossref

PubMed

2020

Author(s)

Direct determination of optimal pair-natural orbitals in a real-space representation: The second-order Moller–Plesset energy

I. INTRODUCTION

II. FORMALISM

A. Multiresolution spectral element representation

B. MRA-MP2 and MRA-PNO-MP2-F12

1. MRA-MP2

2. MRA-PNO-MP2-F12

C. Explicitly correlated ansatz for MRA-PNO-MP2

1. MP2-F12 energies

D. Optimizing PNOs and their amplitudes

1. Optimization of the PNO amplitudes

2. Optimization of the PNOs

III. COMPUTATIONAL DETAILS

IV. RESULTS

A. Saturating the PNO subspace: The helium atom

B. Saturating the PNO subspace: The methane molecule

C. Energy computation: The projection formula vs the Hylleraas functional

D. Reference orbitals: Canonical vs localized

E. Error growth with the system size: n-alkanes

F. Error growth with the system size: n-silanes

G. Conventional LCAO PNOs vs optimized MRA PNOs: Comparison of ranks and precision

H. Performance for relative energies: Water clusters

I. Formal scaling and expected computational performance

V. SUMMARY AND PERSPECTIVE

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Projector onto the occupied space	$\hat{O} = \sum_{i \in occ} \| i (⟩ ⟨) i \|$
Projector onto the virtual space	$\hat{Q} = 1 - \hat{O}$
Two-body notation	${\hat{X}}_{1} \equiv \hat{X} (1), {\hat{X}}_{12} \equiv \hat{X} (1,2)$
Two-body Fock operator (same for $\hat{K}$ and $\hat{J}$ ⁠)	${\hat{F}}_{12} = {\hat{F}}_{1} + {\hat{F}}_{2}$
Fock matrix	$F_{p q} = ⟨) p \| \hat{F} \| q (⟩$
Fock operator	$\hat{F} = \hat{T} + \hat{J} - \hat{K}$
Potential part	${\hat{v}}^{F} \equiv \hat{F} - \hat{T}$
Particle and orbital permutation	${\hat{P}}_{12}, {\hat{P}}_{p q},$
Explicit correlation factor	$f_{12} \equiv f (1,2) = 1 - \exp (- r_{12})$
Coulomb potential	$g_{12} = r_{12}^{- 1}$
Hylleraas functional	$L^{(2)}$ see Eqs. (1), (28), and (10)
BSH Green’s function	${\hat{G}}_{μ}$ see Ref. 36
MP2 pair function	\|τ_ij⟩
MRA with accuracy 10^−x and PNO threshold of 10^−δ	MRA(x, δ)
aug-cc-pVXZ basis-set	aXZ

Direct determination of optimal pair-natural orbitals in a real-space representation: The second-order Moller–Plesset energy

I. INTRODUCTION

II. FORMALISM

A. Multiresolution spectral element representation

B. MRA-MP2 and MRA-PNO-MP2-F12

1. MRA-MP2

2. MRA-PNO-MP2-F12

C. Explicitly correlated ansatz for MRA-PNO-MP2

1. MP2-F12 energies

D. Optimizing PNOs and their amplitudes

1. Optimization of the PNO amplitudes

2. Optimization of the PNOs

III. COMPUTATIONAL DETAILS

IV. RESULTS

A. Saturating the PNO subspace: The helium atom

B. Saturating the PNO subspace: The methane molecule

C. Energy computation: The projection formula vs the Hylleraas functional

D. Reference orbitals: Canonical vs localized

E. Error growth with the system size: n-alkanes

F. Error growth with the system size: n-silanes

G. Conventional LCAO PNOs vs optimized MRA PNOs: Comparison of ranks and precision

H. Performance for relative energies: Water clusters

I. Formal scaling and expected computational performance

V. SUMMARY AND PERSPECTIVE

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only