Monte Carlo explicitly correlated many-body Green’s function theory

Johnson, Cole M.; Doran, Alexander E.; Ten-no, Seiichiro L.; Hirata, So

doi:10.1063/1.5054610

A highly scalable stochastic algorithm is proposed and implemented for computing the basis-set-incompleteness correction to the diagonal, frequency-independent self-energy of the second-order many-body Green’s function (GF2) theory within the explicitly correlated (F12) formalism. The 6-, 9-, 12-, and 15-dimensional integrals comprising the F12 correction are directly evaluated by the Monte Carlo method using appropriate weight functions for importance sampling. The method is naturally and easily parallelized, involves minimal memory space and no disk I/O, and can use virtually any mathematical form of a correlation factor. Its computational cost to correct all ionization energies (IEs) is observed to increase as the fourth power of system size, as opposed to the fifth power in the case of the deterministic counterparts. The GF2 calculations and their F12 corrections for the first IEs of C₆₀ and C₇₀ were executed on 128 graphical processing units (GF2) and 896 central processing units (F12), respectively, to reach the results with statistical errors of 0.04 eV or less. They showed that the basis-set-incompleteness (from aug-cc-pVDZ) accounts for only 50%–60% of the deviations from experiments, suggesting the significance of higher-order perturbation corrections.

I. INTRODUCTION

One-particle many-body Green’s function (MBGF) theory^1–8 provides a converging series of perturbation approximations to electron-binding energies; the accurate determination of which is key to understanding many processes in chemistry, biology, and materials science. Green’s function also serves as the propagator in quantum field theory, underlying nearly all quantum many-body physics ranging from condensed matter physics to particle physics. Furthermore, it is the mathematical kernel of quantum transport, scattering, and embedding theories.

Two of the present authors with co-authors⁹ implemented general-order algorithms of the Feynman–Dyson perturbation series of MBGF and quantified its convergence rate and the accuracy of the diagonal and frequency-independent approximations to the self-energy. The study showed that the convergence is surprisingly slow, making higher-order corrections potentially much more important than in many-body perturbation theory (MBPT) for the ground-state correlation energy.

This slow perturbation-order convergence is particularly troublesome in view of the high-rank polynomial size dependence of the costs of these higher-order methods. The problem size is the product of the system’s spatial extent and basis-set size. This issue is, therefore, further compounded by the notoriously slow basis-set convergence^10–16 of quantities obtained from any ab initio electron-correlation theory including MBGF.

The slow basis-set convergence, however, was tackled by one of the present authors with a co-author,¹⁷ who developed an explicitly correlated (F12) extension of the second-order MBGF (GF2) in the diagonal, frequency-independent approximation to the self-energy (GF2-F12). It allows the complete-basis-set (CBS) results to be obtained with a relatively small basis set on a massively parallel computer. The diagonal and frequency-independent approximations were later lifted in GF2-F12 by Pavešević et al.¹⁸

In this article, we develop a highly scalable stochastic algorithm of the GF2-F12 method of Ohnishi and Ten-no,¹⁷ which we call the Monte Carlo GF2-F12 or MC-GF2-F12. It is naturally and easily parallelized and its computational cost increases less steeply with the number of basis functions. Its formalism and algorithm are a straightforward combination of the Monte Carlo second-order Green’s function (MC-GF2) method¹⁹ and Monte Carlo explicitly correlated second-order many-body perturbation (MC-MP2-F12) method^20,21 reported earlier. We show that MC-GF2-F12 can compute ionization energies (IEs) near their CBS limits with a relatively small basis set for large conjugated molecules with statistical errors that are comparable or smaller than typical experimental errors. The method also allows the use of virtually any correlation factor (geminal). The details are given in the following.

II. FORMALISMS

A. GF2-F12

Following the derivation and notation of Ohnishi and Ten-no,¹⁷ we write the second-order perturbation approximation of the xth occupied-orbital energy, $ϵ_{x}^{GF 2}$ ⁠, as a sum of three parts,

ϵ_{x}^{GF 2} = ϵ_{x}^{HF} + Σ_{x}^{OR} + Σ_{x}^{CD},

(1)

where $ϵ_{x}^{HF}$ is the canonical Hartree–Fock (HF) orbital energy and $Σ_{x}^{OR}$ and $Σ_{x}^{CD}$ are the so-called orbital/pair-relaxation and correlation-difference contributions, respectively, to the second-order self-energy in the diagonal, frequency-independent approximation.⁹ We call this method GF2 in this article. Its contributions are given as

Σ_{x}^{OR} = \sum_{i, j}^{occ.} \sum_{a}^{vir.} \frac{⟨ i j | x a ⟩ (2 ⟨ x a | i j ⟩ - ⟨ a x | i j ⟩)}{ϵ_{x}^{HF} + ϵ_{a}^{HF} - ϵ_{i}^{HF} - ϵ_{j}^{HF}},

(2)

Σ_{x}^{CD} = \sum_{j}^{occ.} \sum_{a, b}^{vir.} \frac{⟨ x j | a b ⟩ (2 ⟨ a b | x j ⟩ - ⟨ b a | x j ⟩)}{ϵ_{x}^{HF} + ϵ_{j}^{HF} - ϵ_{a}^{HF} - ϵ_{b}^{HF}},

(3)

where i, j, k, l, m, n, and x denote an occupied spatial orbital in a closed-shell molecule, while a and b denote a virtual spatial orbital.

As shown numerically by Ohnishi and Ten-no,¹⁷ $Σ_{x}^{OR}$ converges rapidly toward the CBS limit, while the basis-set convergence of $Σ_{x}^{CD}$ is rather slow owing to its double virtual summation. Therefore, only $Σ_{x}^{CD}$ needs to be corrected for the basis-set incompleteness. Rewriting $Σ_{x}^{CD}$ as

Σ_{x}^{CD} = \sum_{j}^{occ.} e_{x j},

(4)

with

e_{x j} = \sum_{a, b}^{vir.} \frac{⟨ x j | a b ⟩ (2 ⟨ a b | x j ⟩ - ⟨ b a | x j ⟩)}{ϵ_{x}^{HF} + ϵ_{j}^{HF} - ϵ_{a}^{HF} - ϵ_{b}^{HF}},

(5)

we notice that second-order many-body perturbation (MP2) energy and its F12 correction,^22–24 E^MP2 and E^F12, are also written as

E^{MP 2} = \sum_{i, j}^{occ.} e_{i j},

(6)

E^{F 12} = \sum_{i, j}^{occ.} d_{i j},

(7)

where d_ij is the correction to the second-order pair energy of the occupied orbital pair ij. Hence, the F12-corrected (GF2-F12) value of the GF2 orbital energy, $ϵ_{x}^{GF 2 - F 12}$ ⁠, can simply be obtained¹⁷ as

ϵ_{x}^{GF 2 - F 12} = ϵ_{x}^{HF} + Σ_{x}^{OR} + Σ_{x}^{CD} + Σ_{x}^{F 12},

(8)

Σ_{x}^{F 12} = \sum_{j}^{occ.} d_{x j},

(9)

where $Σ_{x}^{F 12}$ is the F12 correction to $Σ_{x}^{CD}$ ⁠.

The last equation suggests that, in GF2-F12, we can reuse much of the MP2-F12 formalisms and implementations. The correction, d_xj, is given by the established formalisms^{20,21,24–27} of MP2-F12 as

d_{x j} = 2 d_{x j}^{V} + d_{x j}^{B} + d_{x j}^{X},

(10)

with

d_{x j}^{V} = \sum_{m, n}^{occ.} V_{m n}^{x j} (2 t_{x j}^{m n} - t_{j x}^{m n}),

(11)

d_{x j}^{B} = \sum_{k, l, m, n}^{occ.} t_{m n}^{k l} B_{k l}^{x j} (2 t_{x j}^{m n} - t_{j x}^{m n}),

(12)

d_{x j}^{X} = - \sum_{k, l, m, n}^{occ.} (ϵ_{m} + ϵ_{n}) t_{m n}^{k l} X_{k l}^{x j} (2 t_{x j}^{m n} - t_{j x}^{m n}),

(13)

where the geminal amplitudes, $t_{j x}^{m n}$ ⁠, are held fixed at the values dictated by the cusp conditions,^22,28,29 as per Ten-no’s SP ansatz²⁷

t_{m n}^{x j} = \frac{3}{8} δ_{m x} δ_{n j} + \frac{1}{8} δ_{m j} δ_{n x} .

(14)

Other factors are designated in the standard notations,²¹ whose definitions are written as

V_{m n}^{i j} = ⟨) i j | r_{12}^{- 1} {\hat{Q}}_{12} f_{12} | m n (⟩,

(15)

B_{m n}^{i j} = ⟨) i j | f_{12} {\hat{Q}}_{12} ({\hat{F}}_{1} + {\hat{F}}_{2}) {\hat{Q}}_{12} f_{12} | m n (⟩,

(16)

X_{m n}^{i j} = ⟨) i j | f_{12} {\hat{Q}}_{12} f_{12} | m n (⟩,

(17)

where f₁₂ is the correlation factor or geminal (an explicit function of r₁₂), ${\hat{F}}_{n}$ is the Fock operator of electron n, and ${\hat{Q}}_{12}$ is the strong-orthogonality projector,^30–33

{\hat{Q}}_{12} = (1 - Ô_{1}) (1 - Ô_{2}) - {\hat{V}}_{1} {\hat{V}}_{2} .

(18)

Here, Ô_n and ${\hat{V}}_{n}$ are the projectors onto the occupied and virtual orbital spaces of electron n, respectively.

Utilizing the generalized and extended Brillouin conditions^21,24,27,34 (which consolidate the B and X terms into the combined “BX” term) and explicitly symmetrizing each summation with respect to an interchange of x and j (a necessity in its Monte Carlo integration), we arrive at the working equation for the F12 correction

Σ_{x}^{F 12} \equiv \sum_{j}^{occ.} d_{x j} = 2 Σ_{x}^{V} + Σ_{x}^{B X},

(19)

where

Σ_{x}^{V} = \sum_{j}^{occ.} \frac{5}{16} (V_{x j}^{x j} + V_{j x}^{j x}) - \sum_{j}^{occ.} \frac{1}{16} (V_{j x}^{x j} + V_{x j}^{j x}),

(20)

\begin{matrix} Σ_{x}^{B X} & = & \sum_{j}^{occ.} \frac{7}{64} \{{(B X)}_{x j}^{x j} + {(B X)}_{j x}^{j x}\} + \sum_{j}^{occ.} \frac{1}{64} \{{(B X)}_{j x}^{x j} + {(B X)}_{x j}^{j x}\}, \end{matrix}

(21)

with

{(B X)}_{m n}^{i j} = ⟨) i j | f_{12} {\hat{Q}}_{12} [{\hat{F}}_{1} + {\hat{F}}_{2}, f_{12}] | m n (⟩ .

(22)

When the correlation factor, f₁₂, is optimized by minimizing this correction, the following identity holds (as in MP2-F12):²¹

Σ_{x}^{F 12} \equiv 2 Σ_{x}^{V} + Σ_{x}^{B X} = Σ_{x}^{V} .

(23)

The same identity still holds approximately but accurately by merely optimizing a few parameters (γ in the following examples) in the correlation factor insofar as its shape can closely mimic the correlation hole.^21,27,35 Such correlation factors include the Slater-type geminal (STG),²⁶ Yukawa–Coulomb (YC) geminal,^36,37 and Slater–Jastrow (SJ) factor³⁸

f_{12}^{STG} = {1 - \exp (- γ r_{12})} / γ,

(24)

f_{12}^{YC} = {2 \exp (- γ r_{12}) - 2} / (γ^{2} r_{12}),

(25)

f_{12}^{SJ} = \exp {r_{12} / (1 + γ r_{12})} .

(26)

We call the GF2-F12 method based on Eq. (19) the VBX approximation, whereas the one using Eq. (23), the V approximation. The magnitude of the VBX correction serves as an upper bound for the true basis-set-incompleteness correction.

The foregoing derivation was taken from Ohnishi and Ten-no.¹⁷

B. MC-GF2-F12

Here, we illustrate the key derivation steps of the working equations of the Monte Carlo algorithm of GF2-F12 (MC-GF2-F12) (see Ref. 21 for the closely related MC-MP2-F12 formalism). Expanding ${\hat{Q}}_{12}$ ⁠, i.e., substituting Eq. (18) into Eq. (20), we have

Σ_{x}^{V} = Σ_{x}^{V, 2 e} + Σ_{x}^{V, 3 e} + Σ_{x}^{V, 4 e},

(27)

with

\begin{matrix} Σ_{x}^{V, 2 e} & = & \frac{5}{16} \sum_{j}^{occ.} (⟨) x j | \frac{f_{12}}{r_{12}} | x j (⟩ + ⟨) j x | \frac{f_{12}}{r_{12}} | j x (⟩) \\ - \frac{1}{16} \sum_{j}^{occ.} (⟨) x j | \frac{f_{12}}{r_{12}} | j x (⟩ + ⟨) j x | \frac{f_{12}}{r_{12}} | x j (⟩), \end{matrix}

(28)

\begin{matrix} Σ_{x}^{V, 3 e} & = & - \frac{5}{8} \sum_{j, k}^{occ.} (⟨) x j k | \frac{f_{23}}{r_{12}} | k j x (⟩ + ⟨) j x k | \frac{f_{23}}{r_{12}} | k x j (⟩) \\ + \frac{1}{8} \sum_{j, k}^{occ.} (⟨) x j k | \frac{f_{23}}{r_{12}} | k x j (⟩ + ⟨) j x k | \frac{f_{23}}{r_{12}} | k j x (⟩), \end{matrix}

(29)

\begin{matrix} Σ_{x}^{V, 4 e} & = & \frac{5}{16} \sum_{j, k, l}^{occ.} (⟨) x j k l | \frac{f_{34}}{r_{12}} | k l x j (⟩ + ⟨) j x k l | \frac{f_{34}}{r_{12}} | k l j x (⟩) \\ - \frac{1}{16} \sum_{j, k, l}^{occ.} (⟨) x j k l | \frac{f_{34}}{r_{12}} | k l j x (⟩ + ⟨) j x k l | \frac{f_{34}}{r_{12}} | k l x j (⟩) \\ - \frac{5}{16} \sum_{j}^{occ.} \sum_{a, b}^{vir.} (⟨) x j a b | \frac{f_{34}}{r_{12}} | a b x j (⟩ + ⟨) j x a b | \frac{f_{34}}{r_{12}} | a b j x (⟩) \\ + \frac{1}{16} \sum_{j}^{occ.} \sum_{a, b}^{vir.} (⟨) x j a b | \frac{f_{34}}{r_{12}} | a b j x (⟩ + ⟨) j x a b | \frac{f_{34}}{r_{12}} | a b x j (⟩), \end{matrix}

(30)

where the superscript on Σ indicates the number of electrons involved in the integral (thus its dimension) and x denotes the occupied spatial orbital whose energy is being corrected by F12. These are then converted into an MC integrable form as

Σ_{x}^{V, 2 e} = \iint d r_{1} d r_{2} F_{x}^{V, 2 e} (r_{1}, r_{2}),

(31)

Σ_{x}^{V, 3 e} = ∭ d r_{1} d r_{2} d r_{3} F_{x}^{V, 3 e} (r_{1}, r_{2}, r_{3}),

(32)

Σ_{x}^{V, 4 e} = ⨌ d r_{1} d r_{2} d r_{3} d r_{4} F_{x}^{V, 4 e} (r_{1}, r_{2}, r_{3}, r_{4}),

(33)

where

\begin{matrix} F_{x}^{V, 2 e} (r_{1}, r_{2}) & = & \frac{5}{16} \frac{f_{12} () X_{11} O_{22} + O_{11} X_{22}))}{r_{12}} \\ - \frac{1}{16} \frac{f_{12} () X_{12} O_{21} + O_{12} X_{21}))}{r_{12}}, \end{matrix}

(34)

\begin{matrix} F_{x}^{V, 3 e} (r_{1}, r_{2}, r_{3}) & = & - \frac{5}{8} \frac{f_{23} () X_{13} O_{22} O_{31} + O_{13} X_{22} O_{31}))}{r_{12}} \\ + \frac{1}{8} \frac{f_{23} () X_{12} O_{23} O_{31} + O_{12} X_{23} O_{31}))}{r_{12}}, \end{matrix}

(35)

\begin{matrix} F_{x}^{V, 4 e} (r_{1}, r_{2}, r_{3}, r_{4}) \\ = \frac{5}{16} \frac{f_{34} () X_{13} O_{24} O_{31} O_{42} + O_{13} X_{24} O_{31} O_{42}))}{r_{12}} \\ - \frac{1}{16} \frac{f_{34} () X_{14} O_{23} O_{31} O_{42} + O_{14} X_{23} O_{31} O_{42}))}{r_{12}} \\ - \frac{5}{16} \frac{f_{34} () X_{13} O_{24} V_{31} V_{42} + O_{13} X_{24} V_{31} V_{42}))}{r_{12}} \\ + \frac{1}{16} \frac{f_{34} () X_{14} O_{23} V_{31} V_{42} + O_{14} X_{23} V_{31} V_{42}))}{r_{12}} . \end{matrix}

(36)

Here, O, V, and X are two-electron, six-dimensional, real-space pair functions given by

O_{p q} = \sum_{i}^{o c c .} φ_{i}^{*} (r_{p}) φ_{i} (r_{q}),

(37)

V_{p q} = \sum_{a}^{v i r .} φ_{a}^{*} (r_{p}) φ_{a} (r_{q}),

(38)

and

X_{p q} = φ_{x}^{*} (r_{p}) φ_{x} (r_{q}),

(39)

where φ_p(r_q) is the pth spatial orbital of electron q.

Likewise, Eq. (21) becomes a sum of six terms,

Σ_{x}^{B X} = Σ_{x}^{T, 2 e} + Σ_{x}^{T, 3 e} + Σ_{x}^{T, 4 e} + Σ_{x}^{K, 3 e} + Σ_{x}^{K, 4 e} + Σ_{x}^{K, 5 e},

(40)

where “T” and “K” stand for the kinetic and exchange contributions, respectively, of the Fock operator in the BX integral, Eq. (22),

[{\hat{F}}_{1} + {\hat{F}}_{2}, f_{12}] = [{\hat{T}}_{1} + {\hat{T}}_{2}, f_{12}] - [{\hat{K}}_{1} + {\hat{K}}_{2}, f_{12}] .

(41)

Here, ${\hat{T}}_{n}$ and ${\hat{K}}_{n}$ stand for the kinetic and exchange operators, respectively, for electron n (the other operators in $\hat{F}$ commute with f₁₂). Note that the dimensions of the exchange terms are one-electron higher than those of the kinetic terms. This is owing to the presence of a projector in the exchange operator, which, when expanded, increases the integral’s dimension (like ${\hat{Q}}_{12}$ ⁠). The MC-integrable expression then reads

\begin{matrix} Σ_{x}^{T, 2 e} & = & \iint d r_{1} d r_{2} F_{x}^{T 1,2 e} (r_{1}, r_{2}) \\ + \iint d r_{1} d r_{2} F_{x}^{T 2,2 e} (r_{1}, r_{2}), \end{matrix}

(42)

\begin{matrix} Σ_{x}^{T, 3 e} & = & ∭ d r_{1} d r_{2} d r_{3} F_{x}^{T 1,3 e} (r_{1}, r_{2}, r_{3}) \\ + ∭ d r_{1} d r_{2} d r_{3} F_{x}^{T 2,3 e} (r_{1}, r_{2}, r_{3}), \end{matrix}

(43)

\begin{matrix} Σ_{x}^{T, 4 e} & = & ⨌ d r_{1} d r_{2} d r_{3} d r_{4} F_{x}^{T 1,4 e} (r_{1}, r_{2}, r_{3}, r_{4}) \\ + ⨌ d r_{1} d r_{2} d r_{3} d r_{4} F_{x}^{T 2,4 e} (r_{1}, r_{2}, r_{3}, r_{4}), \end{matrix}

(44)

Σ_{x}^{K, 3 e} = ∭ d r_{1} d r_{2} d r_{3} F_{x}^{K, 3 e} (r_{1}, r_{2}, r_{3}),

(45)

Σ_{x}^{K, 4 e} = ⨌ d r_{1} d r_{2} d r_{3} d r_{4} F_{x}^{K, 4 e} (r_{1}, r_{2}, r_{3}, r_{4}),

(46)

Σ_{x}^{K, 5 e} = \int \dots \int d r_{1} \dots d r_{5} F_{x}^{K, 5 e} (r_{1}, r_{2}, r_{3}, r_{4}, r_{5}) .

(47)

The two-electron kinetic integrands are given by

\begin{matrix} F_{x}^{T 1,2 e} (r_{1}, r_{2}) \\ = \frac{7}{64} \frac{f_{12} f_{12}^{(a)} (X_{11} O_{22} + O_{11} X_{22})}{r_{12}} \\ + \frac{7}{64} \frac{f_{12} f_{12}^{(c)} (X_{11} O_{22}^{'} + O_{11} X_{22}^{'} - X_{11}^{'} O_{22} - O_{11}^{'} X_{22})}{r_{12}} \\ + \frac{1}{64} \frac{f_{12} f_{12}^{(a)} (X_{12} O_{21} + O_{12} X_{21})}{r_{12}} \\ + \frac{1}{64} \frac{f_{12} f_{12}^{(c)} (X_{21} O_{12}^{'} + O_{21} X_{12}^{'} - X_{21}^{'} O_{12} - O_{21}^{'} X_{12})}{r_{12}}, \end{matrix}

(48)

\begin{matrix} F_{x}^{T 2,2 e} (r_{1}, r_{2}) \\ = \frac{7}{64} f_{12} f_{12}^{(b)} (X_{11} O_{22} + O_{11} X_{22}) \\ + \frac{7}{64} f_{12} f_{12}^{(d)} (X_{11} O_{22}^{'} + O_{11} X_{22}^{'} - X_{11}^{'} O_{22} - O_{11}^{'} X_{22}) \\ + \frac{1}{64} f_{12} f_{12}^{(b)} (X_{12} O_{21} + O_{12} X_{21}) \\ + \frac{1}{64} f_{12} f_{12}^{(d)} (X_{21} O_{12}^{'} + O_{21} X_{12}^{'} - X_{21}^{'} O_{12} - O_{21}^{'} X_{12}), \end{matrix}

(49)

the three-electron kinetic integrands by

\begin{matrix} F_{x}^{T 1,3 e} (r_{1}, r_{2}, r_{3}) \\ = - \frac{7}{32} \frac{f_{23} f_{12}^{(a)} (X_{31} O_{22} + O_{31} X_{22}) O_{13}}{r_{12}} \\ - \frac{7}{32} \frac{f_{23} f_{12}^{(c)} (X_{22}^{'} O_{31} + O_{22}^{'} X_{31} - X_{31}^{'} O_{22} + O_{31}^{'} X_{22}) O_{13}}{r_{12}} \\ - \frac{1}{32} \frac{f_{23} f_{12}^{(a)} (X_{32} O_{21} + O_{32} X_{21}) O_{13}}{r_{12}} \\ - \frac{1}{32} \frac{f_{23} f_{12}^{(c)} (X_{32}^{'} O_{21} + O_{32}^{'} X_{21} - X_{21}^{'} O_{32} - O_{21}^{'} X_{32}) O_{13}}{r_{12}}, \end{matrix}

(50)

\begin{matrix} F_{x}^{T 2,3 e} (r_{1}, r_{2}, r_{3}) & = & - \frac{7}{32} f_{23} f_{12}^{(b)} (X_{31} O_{22} + O_{31} X_{22}) O_{13} - \frac{7}{32} f_{23} f_{12}^{(d)} (X_{22}^{'} O_{31} + O_{22}^{'} X_{31} - X_{31}^{'} O_{22} - O_{31}^{'} X_{22}) O_{13} \\ - \frac{1}{32} f_{23} f_{12}^{(b)} (X_{32} O_{21} + O_{32} X_{21}) O_{13} - \frac{1}{32} f_{23} f_{12}^{(d)} (X_{32}^{'} O_{21} + O_{32}^{'} X_{21} - X_{21}^{'} O_{32} - O_{21}^{'} X_{32}) O_{13}, \end{matrix}

(51)

and the four-electron kinetic integrands by

\begin{matrix} F_{x}^{T 1,4 e} (r_{1}, r_{2}, r_{3}, r_{4}) & = & \frac{7}{64} \frac{f_{34} f_{12}^{(a)} (X_{31} O_{42} + O_{31} X_{42}) O_{13} O_{24}}{r_{12}} + \frac{7}{64} \frac{f_{34} f_{12}^{(c)} (X_{42}^{'} O_{31} + O_{42}^{'} X_{31} - X_{31}^{'} O_{42} - O_{31}^{'} X_{42}) O_{13} O_{24}}{r_{12}} \\ + \frac{1}{64} \frac{f_{34} f_{12}^{(a)} (X_{32} O_{41} + O_{32} X_{41}) O_{13} O_{24}}{r_{12}} + \frac{1}{64} \frac{f_{34} f_{12}^{(c)} (X_{32}^{'} O_{41} + O_{32}^{'} X_{41} - X_{41}^{'} O_{32} - O_{41}^{'} X_{32}) O_{13} O_{24}}{r_{12}} \\ - \frac{7}{64} \frac{f_{34} f_{12}^{(a)} (X_{31} O_{42} + O_{31} X_{42}) V_{13} V_{24}}{r_{12}} - \frac{7}{64} \frac{f_{34} f_{12}^{(c)} (X_{42}^{'} O_{31} + O_{42}^{'} X_{31} - X_{31}^{'} O_{42} - O_{31}^{'} X_{42}) V_{13} V_{24}}{r_{12}} \\ - \frac{1}{64} \frac{f_{34} f_{12}^{(a)} (X_{32} O_{41} + O_{32} X_{41}) V_{13} V_{24}}{r_{12}} - \frac{1}{64} \frac{f_{34} f_{12}^{(c)} (X_{32}^{'} O_{41} + O_{32}^{'} X_{41} - X_{41}^{'} O_{32} - O_{41}^{'} X_{32}) V_{13} V_{24}}{r_{12}}, \end{matrix}

(52)

and

\begin{matrix} F_{x}^{T 2,4 e} (r_{1}, r_{2}, r_{3}, r_{4}) & = & \frac{7}{64} f_{34} f_{12}^{(b)} (X_{31} O_{42} + O_{31} X_{42}) O_{13} O_{24} + \frac{7}{64} f_{34} f_{12}^{(d)} (X_{42}^{'} O_{31} + O_{42}^{'} X_{31} - X_{31}^{'} O_{42} - O_{31}^{'} X_{42}) O_{13} O_{24} \\ + \frac{1}{64} f_{34} f_{12}^{(b)} (X_{32} O_{41} + O_{32} X_{41}) O_{13} O_{24} + \frac{1}{64} f_{34} f_{12}^{(d)} (X_{32}^{'} O_{41} + O_{32}^{'} X_{41} - X_{41}^{'} O_{32} - O_{41}^{'} X_{32}) O_{13} O_{24} \\ - \frac{7}{64} f_{34} f_{12}^{(b)} (X_{31} O_{42} + O_{31} X_{42}) V_{13} V_{24} - \frac{7}{64} f_{34} f_{12}^{(d)} (X_{42}^{'} O_{31} + O_{42}^{'} X_{31} - X_{31}^{'} O_{42} - O_{31}^{'} X_{42}) V_{13} V_{24} \\ - \frac{1}{64} f_{34} f_{12}^{(b)} (X_{32} O_{41} + O_{32} X_{41}) V_{13} V_{24} - \frac{1}{64} f_{34} f_{12}^{(d)} (X_{32}^{'} O_{41} + O_{32}^{'} X_{41} - X_{41}^{'} O_{32} - O_{41}^{'} X_{32}) V_{13} V_{24} . \end{matrix}

(53)

The primed real-space pair functions are defined by

X_{p q}^{'} = φ_{x}^{*} (r_{p}) r_{12} \cdot \nabla_{q} φ_{x} (r_{q}),

(54)

O_{p q}^{'} = \sum_{i}^{o c c .} φ_{i}^{*} (r_{p}) r_{12} \cdot \nabla_{q} φ_{i} (r_{q}),

(55)

where r₁₂ = r₁ − r₂. The kinetic-energy integrands such as of Eq. (42) are divided into two terms, $F_{x}^{T 1, n e}$ and $F_{x}^{T 2, n e}$ (n = 2, 3, 4), with the first having a singularity at r₁₂ = |r₁₂| = 0 and the other not, thus requiring different weight functions in MC integrations (see below). They are defined with the quantities derivable from the correlation factor

[{\hat{T}}_{1} + {\hat{T}}_{2}, f_{12}] = \frac{f_{12}^{(a)}}{r_{12}} + f_{12}^{(b)} - \{\frac{f_{12}^{(c)}}{r_{12}} + f_{12}^{(d)}\} r_{12} \cdot (\nabla_{1} - \nabla_{2}) .

(56)

See Appendix A of Ref. 21 for a detailed derivation of equations analogous to Eqs. (48) and (49).

The exchange integrands are

\begin{matrix} F_{x}^{K, 3 e} (r_{1}, r_{2}, r_{3}) \\ = \frac{7}{32} \frac{f_{23} (f_{23} - f_{13})}{r_{12}} (X_{12} O_{33} + O_{12} X_{33}) O_{21} \\ + \frac{1}{32} \frac{f_{23} (f_{23} - f_{13})}{r_{12}} (X_{13} O_{32} + O_{13} X_{32}) O_{21}, \end{matrix}

(57)

\begin{matrix} F_{x}^{K, 4 e} (r_{1}, r_{2}, r_{3}, r_{4}) \\ = - \frac{7}{32} \frac{f_{24} (f_{23} - f_{13})}{r_{12}} (X_{43} O_{21} + O_{43} X_{21}) O_{34} O_{12} \\ - \frac{1}{32} \frac{f_{24} (f_{23} - f_{13})}{r_{12}} (X_{41} O_{23} + O_{41} X_{23}) O_{34} O_{12} \\ - \frac{7}{32} \frac{f_{34} (f_{14} - f_{24})}{r_{12}} (X_{32} O_{44} + O_{32} X_{44}) O_{13} O_{21} \\ - \frac{1}{32} \frac{f_{34} (f_{14} - f_{24})}{r_{12}} (X_{34} O_{42} + O_{34} X_{42}) O_{13} O_{21}, \end{matrix}

(58)

\begin{matrix} F_{x}^{K, 5 e} (r_{1}, r_{2}, r_{3}, r_{4}, r_{5}) \\ = \frac{7}{32} \frac{f_{35} (f_{14} - f_{24})}{r_{12}} (X_{32} O_{54} + O_{32} X_{54}) O_{31} O_{45} O_{21} \\ + \frac{1}{32} \frac{f_{35} (f_{14} - f_{24})}{r_{12}} (X_{34} O_{52} + O_{34} X_{52}) O_{31} O_{45} O_{21} \\ - \frac{7}{32} \frac{f_{35} (f_{14} - f_{24})}{r_{12}} (X_{32} O_{54} + O_{32} X_{54}) V_{31} V_{45} O_{21} \\ - \frac{1}{32} \frac{f_{35} (f_{14} - f_{24})}{r_{12}} (X_{34} O_{52} + O_{34} X_{52}) V_{31} V_{45} O_{21} . \end{matrix}

(59)

III. MONTE CARLO ALGORITHM

Just like in MC-MP2-F12,^19–21 the GF2-F12 corrections are evaluated by the Monte Carlo (MC) method³⁹ utilizing the redundant-walker algorithm.⁴⁰ In each MC step, the ratio between the integrand and a judiciously chosen weight function is calculated and summed at coordinates distributed randomly but according to the weight function. The weight function must cancel the singularities of the integrand, be always positive, be analytically integrable, and should generally behave like the integrand.^41,42 The random distribution according to the weight function is attained by the Metropolis algorithm.⁴³

The V integrals [Eqs. (31)–(33)] are evaluated²⁰ as

Σ_{x}^{V, 2 e} = \iint d r_{1} d r_{2} F_{x}^{V, 2 e} (r_{1}, r_{2}) \approx \frac{1}{N} \sum_{n = 1}^{N} \frac{F_{x}^{V, 2 e} (r_{1}^{[n]}, r_{2}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]})},

(60)

\begin{matrix} Σ_{x}^{V, 3 e} & = & ∭ d r_{1} d r_{2} d r_{3} F_{x}^{V, 3 e} (r_{1}, r_{2}, r_{3}) \\ \approx & \frac{1}{N} \sum_{n = 1}^{N} \frac{1}{m} \sum_{k = 1}^{m} \frac{F_{x}^{V, 3 e} (r_{1}^{[n]}, r_{2}^{[n]}, r_{3 k}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]}) w_{1 e} (r_{3 k}^{[n]})}, \end{matrix}

(61)

\begin{matrix} Σ_{x}^{V, 4 e} & = & ⨌ d r_{1} d r_{2} d r_{3} d r_{4} F_{x}^{V, 4 e} (r_{1}, r_{2}, r_{3}, r_{4}) \\ \approx & \frac{1}{N} \sum_{n = 1}^{N} \frac{2!}{m (m - 1)} \\ \times \sum_{k = 1}^{m - 1} \sum_{l = k + 1}^{m} \frac{F_{x}^{V, 4 e} (r_{1}^{[n]}, r_{2}^{[n]}, r_{3 k}^{[n]}, r_{4 l}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]}) w_{1 e} (r_{3 k}^{[n]}) w_{1 e} (r_{4 l}^{[n]})}, \end{matrix}

(62)

where N is the total number of MC steps and m is the number of “redundant walkers” (see below for the definition).

There are two different types of walkers used for these integrands. For strongly correlated coordinates, $(r_{1}^{[n]}, r_{2}^{[n]})$ ⁠, the weight function that guides their propagation is given by

w_{2 e} (r_{1}, r_{2}) = \frac{1}{N_{2 e}} \frac{g (r_{1}) g (r_{2})}{r_{12}},

(63)

with g(r) being a sum of atom-centered s-type Gaussian-type orbitals,

g (r) = \sum_{α = 1}^{n_{λ}} \sum_{A = 1}^{n_{atom}} c_{A}^{(α)} e^{- ζ_{A}^{(α)} | r - r_{A} |^{2}},

(64)

where n_atom is the number of atoms, r_A is the position of the Ath atom, and n_λ is the number of Gaussian functions per atom. For uncoupled coordinates, $r_{3 k}^{[n]}$ or $r_{4 l}^{[n]}$ ⁠, the weight function is

w_{1 e} (r_{3}) = \frac{g (r_{3})}{N_{1 e}} .

(65)

With this choice of g(r), normalization coefficients, N_2e and N_1e, are determined analytically.⁴⁴ Usually two Gaussians per atom (n_λ = 2) suffice; one tight and one diffuse to expand the electron density. In the applications to the IEs of C₆₀ and C₇₀ (see below), a third is added to account for the large size of the salient molecular orbitals.

The reason for this separation of the walker coordinates (between r₁, r₂ and r₃ and r₄) can be understood by inspecting the structure of the $F_{3 e}^{V}$ [Eq. (35)] and $F_{4 e}^{V}$ [Eq. (36)] integrands. The coordinates of electrons 1 and 2 are strongly coupled through the singular $r_{12}^{- 1}$ factor in these integrands. This factor should, therefore, be canceled by the same factor in the weight function. Electrons 3 and 4 have no such correlation and the corresponding weight function can be an independent function of each electron’s coordinate.

In the redundant-walker algorithm,^20,40 m “redundant” one-electron walkers are propagated according to the weight function, $w_{1 e}$ (r), at a cost that is m times that of propagating the minimal two one-electron walkers in integrating $Σ_{x}^{V, 4 e}$ [Eq. (62)]. With m walkers instead of two, there are m(m − 1)/2 distinct ways of choosing the two walker coordinates that can be substituted in the integrand of Eq. (62). This increases the sampling number by a factor of m(m − 1)/2 at an m-fold increase in the sampling cost, a net increase in the sampling efficiency by (m − 1)/2. However, Eqs. (60) and (61) are unaffected by the algorithm, making the overall efficiency boost difficult to predict.

Similarly, the T integrals [Eqs. (42)–(44)] are evaluated as

\begin{matrix} Σ_{x}^{T, 2 e} & \approx & \frac{1}{N} \sum_{n = 1}^{N} \frac{F_{x}^{T 1,2 e} (r_{1}^{[n]}, r_{2}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]})} \\ + \frac{1}{N} \sum_{n = 1}^{N} \frac{2!}{m (m - 1)} \sum_{k = 1}^{m - 1} \sum_{l = k + 1}^{m} \frac{F_{x}^{T 2,2 e} (r_{1 k}^{[n]}, r_{2 l}^{[n]})}{w_{1 e} (r_{1 k}^{[n]}) w_{1 e} (r_{2 l}^{[n]})}, \end{matrix}

(66)

\begin{matrix} Σ_{x}^{T, 3 e} & \approx & \frac{1}{N} \sum_{n = 1}^{N} \frac{1}{m} \sum_{k = 1}^{m} \frac{F_{x}^{T 1,3 e} (r_{1}^{[n]}, r_{2}^{[n]}, r_{3 k}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]}) w_{1 e} (r_{3 k}^{[n]})} \\ + \frac{1}{N} \sum_{n = 1}^{N} \frac{3!}{m (m - 1) (m - 2)} \\ \times \sum_{k = 1}^{m - 2} \sum_{l = k + 1}^{m - 1} \sum_{h = l + 1}^{m} \frac{F_{x}^{T 2,3 e} (r_{1 k}^{[n]}, r_{2 l}^{[n]}, r_{3 h}^{[n]})}{w_{1 e} (r_{1 k}^{[n]}) w_{1 e} (r_{2 l}^{[n]}) w_{1 e} (r_{3 h}^{[n]})}, \end{matrix}

(67)

and

\begin{matrix} Σ_{x}^{T, 4 e} & \approx & \frac{1}{N} \sum_{n = 1}^{N} \frac{2!}{m (m - 1)} \\ \times \sum_{k = 1}^{m - 1} \sum_{l = k + 1}^{m} \frac{F_{x}^{T 1,4 e} (r_{1}^{[n]}, r_{2}^{[n]}, r_{3 k}^{[n]}, r_{4 l}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]}) w_{1 e} (r_{3 k}^{[n]}) w_{1 e} (r_{4 l}^{[n]})} \\ + \frac{1}{N} \sum_{n = 1}^{N} \frac{4!}{m (m - 1) (m - 2) (m - 3)} \sum_{k = 1}^{m - 3} \sum_{l = k + 1}^{m - 2} \sum_{h = l + 1}^{m - 1} \sum_{i = h + 1}^{m} \\ \times \frac{F_{x}^{T 2,4 e} (r_{1 k}^{[n]}, r_{2 l}^{[n]}, r_{3 h}^{[n]}, r_{4 i}^{[n]})}{w_{1 e} (r_{1 k}^{[n]}) w_{1 e} (r_{2 l}^{[n]}) w_{1 e} (r_{3 h}^{[n]}) w_{1 e} (r_{4 i}^{[n]})}, \end{matrix}

(68)

of which all factors are already defined.

The K integrals [Eqs. (45)–(47)] are computed as

Σ_{x}^{K, 3 e} \approx \frac{1}{N} \sum_{n = 1}^{N} \frac{1}{m} \sum_{k = 1}^{m} \frac{F_{x}^{K, 3 e} (r_{1}^{[n]}, r_{2}^{[n]}, r_{3 k}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]}) w_{1 e} (r_{3 k}^{[n]})},

(69)

\begin{matrix} Σ_{x}^{K, 4 e} & \approx & \frac{1}{N} \sum_{n = 1}^{N} \frac{2!}{m (m - 1)} \\ \times \sum_{k = 1}^{m - 1} \sum_{l = k + 1}^{m} \frac{F_{x}^{K, 4 e} (r_{1}^{[n]}, r_{2}^{[n]}, r_{3 k}^{[n]}, r_{4 l}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]}) w_{1 e} (r_{3 k}^{[n]}) w_{1 e} (r_{4 l}^{[n]})}, \end{matrix}

(70)

and

\begin{matrix} Σ_{x}^{K, 5 e} & \approx & \frac{1}{N} \sum_{n = 1}^{N} \frac{3!}{m (m - 1) (m - 2)} \sum_{k = 1}^{m - 2} \sum_{l = k + 1}^{m - 1} \sum_{h = l + 1}^{m} \\ \times \frac{F_{x}^{K, 5 e} (r_{1}^{[n]}, r_{2}^{[n]}, r_{3 k}^{[n]}, r_{4 l}^{[n]}, r_{5 h}^{[n]})}{w_{2 e} (r_{1}^{[n]}, r_{2}^{[n]}) w_{1 e} (r_{3 k}^{[n]}) w_{1 e} (r_{4 l}^{[n]}) w_{1 e} (r_{5 h}^{[n]})} . \end{matrix}

(71)

IV. RESULTS AND DISCUSSION

MC-GF2 (Refs. 19 and 45) and MC-GF2-F12 calculations were performed for the first (and occasionally the second) IEs of the molecules listed in Fig. 1 using the aug-cc-pVDZ basis set in the frozen core approximation. The structures of molecules A–O were taken from Ref. 17, whereas those of molecules P and Q came from Ref. 46. The MC-GF2 and F12 calculations were executed simultaneously; for every 8 processors, one was a graphical processing unit (GPU) running MC-GF2 and the other 7 were central processing units (CPUs) performing MC-GF2-F12. For molecules A–O, 64 GPUs (GF2) and 448 CPUs (F12) of the Blue Waters supercomputer were employed; for molecules P and Q, 128 GPUs (GF2) and 896 CPUs (F12) were used. 512 redundant walkers were propagated for MC-GF2, while 40 redundant one- and two-electron walkers were used for MC-GF2-F12. The statistical uncertainty was taken after 7 blocking transformations.⁴⁷

FIG. 1.

View large Download slide

A: benzene; B: naphthalene; C: anthracene; D: pyrene; E: coronene; F: ovalene; G: porphyrin; H: tetraphenylporphyrin; I: thiophene; J: bithiophene; K: terthiophene; L: quaterthiophene; M: quinquethiophene; N: sexithiophene; O: septithiophene; P: C₆₀; Q: C₇₀.

Table I compares the MC-GF2-F12 results for molecules A–O in the V and VBX approximations with the deterministic GF2-F12 results.¹⁷ The MC results were converged to a statistical error of 0.04 eV or less after the number of MC steps indicated in the table, which may be comparable to a typical experimental error (or “chemical accuracy” of 0.043 eV). The VBX data reproduce the CBS limits obtained by the deterministic GF2-F12 method¹⁷ even though the two methods differ slightly in detail as well as in the γ value and basis set used. The mean absolute deviation (MAE) of 0.018 eV between these two methods is, therefore, largely (if not entirely) due to the statistical error (0.017 eV) in the MC-GF2-F12 method. The MAE between the V approximation and deterministic GF2-F12 is distinctly greater (0.049 eV) and is likely a bias caused by the nonvariational nature of the F12 correction in the V approximation exacerbated by the small aug-cc-pVDZ basis set used [recall that the fidelity of the V approximation or the degree to which Eq. (23) is satisfied depends rather sensitively on the mathematical form of the correlation factor and on the value of γ]. It can be said that the MC algorithm works correctly with almost no bias, reproducing the deterministic counterparts within a sufficiently small statistical error, as long as the VBX formalism and an appropriate value of γ are used.

TABLE I.

The first (and second) IEs (in eV) calculated by the deterministic GF2-F12/aug-cc-pVTZ (Ref. 17) and MC-GF2-F12/aug-cc-pVDZ (in either the V or VBX approximation). The values in parentheses are the statistical errors from both MC-GF2 and MC-GF2-F12 steps. A Slater-type geminal with γ = 1.1 was used in the MC-GF2-F12 calculations. N stands for the number of MC steps in the MC-GF2 or MC-GF2-F12 calculation running on 64 GPUs or on 448 CPUs, respectively, of the Blue Waters supercomputer.

		MC-GF2-F12		N/10⁹
Molecule^a	GF2-F12^b	V	VBX	GF2	F12
A	9.235	9.174(01)	9.217(01)	0.51	0.43
B	8.045	7.981(02)	8.025(03)	0.41	0.33
C	7.289	7.231(04)	7.278(05)	0.33	0.22
D	7.314	7.246(05)	7.291(06)	0.31	0.19
E	7.185	7.131(10)	7.177(13)	0.21	0.08
F	6.563	6.496(21)	6.542(26)	0.21	0.06
G (²B_3u)	6.718	6.671(19)	6.713(22)	0.20	0.07
G (²A_u)	7.095	7.019(09)	7.066(11)	0.39	0.13
H (²B_u)	6.155	6.125(27)	6.172(38)	1.48	0.23
H (²A_u)	6.754	6.653(21)	6.691(30)	1.86	0.29
I	8.914	8.857(02)	8.891(04)	0.57	0.45
J	7.747	7.720(07)	7.748(11)	0.42	0.33
K	7.257	7.238(12)	7.269(18)	0.33	0.19
L	7.026	7.001(13)	7.039(21)	0.96	0.44
M	6.825	6.804(16)	6.840(26)	0.96	0.30
N	6.767	6.732(18)	6.775(29)	2.00	0.53
O	6.695	6.677(19)	6.715(29)	1.81	0.39
MAE^c	0	0.049(12)	0.018(17)	⋯	⋯

		MC-GF2-F12		N/10⁹
Molecule^a	GF2-F12^b	V	VBX	GF2	F12
A	9.235	9.174(01)	9.217(01)	0.51	0.43
B	8.045	7.981(02)	8.025(03)	0.41	0.33
C	7.289	7.231(04)	7.278(05)	0.33	0.22
D	7.314	7.246(05)	7.291(06)	0.31	0.19
E	7.185	7.131(10)	7.177(13)	0.21	0.08
F	6.563	6.496(21)	6.542(26)	0.21	0.06
G (²B_3u)	6.718	6.671(19)	6.713(22)	0.20	0.07
G (²A_u)	7.095	7.019(09)	7.066(11)	0.39	0.13
H (²B_u)	6.155	6.125(27)	6.172(38)	1.48	0.23
H (²A_u)	6.754	6.653(21)	6.691(30)	1.86	0.29
I	8.914	8.857(02)	8.891(04)	0.57	0.45
J	7.747	7.720(07)	7.748(11)	0.42	0.33
K	7.257	7.238(12)	7.269(18)	0.33	0.19
L	7.026	7.001(13)	7.039(21)	0.96	0.44
M	6.825	6.804(16)	6.840(26)	0.96	0.30
N	6.767	6.732(18)	6.775(29)	2.00	0.53
O	6.695	6.677(19)	6.715(29)	1.81	0.39
MAE^c	0	0.049(12)	0.018(17)	⋯	⋯

^a

See Fig. 1 and also Ref. 17.

^b

Reference 17.

^c

Mean absolute error from the GF2-F12 values.

In our previous study,²¹ we established the O(n²) scaling of the computational cost (wall time) per MC step with the number of basis functions (n) in an MC-MP2-F12 calculation. This combined with the observations that the relative statistical error increases as O(n) and that the same falls off as N^−1/2 with the number of MC steps (N) led to the overall O(n⁴) scaling of the MC-MP2-F12 method, which is one-rank lower than the O(n⁵) scaling of deterministic MP2-F12 methods with the SP ansatz [the scaling would be even worse and O(n⁶) for deterministic MP2-F12 methods that do not use the SP ansatz]. In this sense, the MC algorithm is guaranteed to outperform deterministic ones in the large n limit, though the crossover point seems rather far.¹⁵

If we determine O(n) IEs simultaneously from a single GF2-F12 calculation (which is realistic for large molecules and band structures in solids in particular), the scaling of the cost per MC step in MC-GF2-F12 would be the same O(n²) as in MC-MP2-F12 because of the algorithmic similarity illustrated by Eqs. (6)–(9). To know the scaling of MC-GF2-F12, we also need to know how the statistical error in an IE behaves with n. Unlike the total (correlation) energy computed in MC-MP2-F12, we demand that the bare statistical error (not the relative error) fall below a certain threshold (such as 0.04 eV) in the case of IEs. Figure 2 clearly shows the O(n) scaling of the bare statistical error, regardless of the details (V versus VBX) of the formalism. This and the O(n²) scaling of the cost per MC step lead to the overall O(n⁴) scaling of MC-GF2-F12 for calculating all IEs. This is one-rank lower than the O(n⁵) scaling of deterministic GF2-F12 calculations.⁴⁸

FIG. 2.

View large Download slide

The statistical error, σ (in E_h), of IEs calculated by MC-GF2-F12 (in the V and VBX approximations) as a function of the number of basis functions n. A linear function of n is overlaid in blue.

A unique advantage of MC-GF2-F12 is its ability to use virtually any correlation factor.^21,35 Table II compares the performance of three of high-performing correlation factors,^21,35 STG, YC, and SJ, as defined by Eqs. (24)–(26), for the first IEs of the water, methane, and benzene molecules. These three correlation factors yield nearly interchangeable F12 corrections (within 0.03 eV of one another) but have vastly different statistical errors. The SJ factor, in particular, suffers from large statistical errors likely because of the greater complexity and number of derivatives appearing in the kinetic-energy operator commutator of Eq. (56). This leaves STG, which happens to be convenient for deterministic F12 algorithms,²⁶ the most attractive choice for the MC algorithms also.

TABLE II.

The F12 corrections (in eV) to the first IEs of the water, methane, and benzene molecules obtained by MC-GF2-F12 in the VBX approximation using the Slater-type geminal (STG), Yukawa–Coulomb (YC) geminal, or Slater–Jastrow (SJ) factor with the γ value given in the corresponding row. The statistical errors taken after 4.59 × 10⁷ MC steps are indicated in parentheses.

Correlation factor	γ	Water	Methane	Benzene
STG [Eq. (24)]	1.1	0.578(01)	0.331(01)	0.402(03)
YC [Eq. (25)]	1.6	0.569(02)	0.327(01)	0.411(10)
SJ [Eq. (26)]	1.1	0.569(07)	0.335(06)	0.384(78)

Correlation factor	γ	Water	Methane	Benzene
STG [Eq. (24)]	1.1	0.578(01)	0.331(01)	0.402(03)
YC [Eq. (25)]	1.6	0.569(02)	0.327(01)	0.411(10)
SJ [Eq. (26)]	1.1	0.569(07)	0.335(06)	0.384(78)

Table III lists the results of the MC-GF2 and MC-GF2-F12 calculations for the first IEs of C₆₀ and C₇₀ using the aug-cc-pVDZ basis set. It should be stressed that MC-GF2 and MC-GF2-F12 calculations are always executable no matter how large a system is insofar as the preceding HF calculation can be completed. This is because the MC algorithms use tiny memory footprint and no disk I/O, while being more demanding of CPU time, although the latter issue is mitigated by the high degree of parallelism and the effective use of GPUs.⁴⁵ The MC algorithms are, therefore, an example of space-time trade-off⁵⁰ (somewhat like Almlöf’s direct self-consistent field algorithm⁵¹), which minimizes memory space at the expense of increased operation cost, thereby lifting a hardware-made limitation in applicability to large problems. Whether such calculations can be converged to within a reasonable statistical error with available resources is a different question, but the answer is affirmative for C₆₀ and C₇₀ using a small fraction of the Blue Waters supercomputer.

TABLE III.

The first IEs (in eV) of C₆₀ and C₇₀ calculated by MC-GF2 and MC-GF2-F12 (in the VBX approximation) with the aug-cc-pVDZ basis set. The values in parentheses are the statistical errors. A Slater-type geminal with γ = 1.1 was used in the MC-GF2-F12 calculations. Calculations were performed on 128 GPUs (GF2) or on 896 CPUs (F12) of the Blue Waters supercomputer. The numbers of MC steps taken were 3.67 × 10⁸ (GF2) and 6.15 × 10⁷ (F12) for C₆₀ and 3.14 × 10⁸ (GF2) and 4.54 × 10⁷ (F12) for C₇₀.

Molecule^a	MC-GF2	MC-GF2-F12	Experiment^b
P (C₆₀)	6.70(2)	7.17(3)	7.64(2)
Q (C₇₀)	6.73(2)	7.17(4)	7.47(2)

^a

See Fig. 1 and also Ref. 46.

^b

Reference 49.

After about 4 × 10⁸ (GF2) and 6 × 10⁷ (F12) MC steps, the statistical errors in both molecules fall below 0.04 eV. The F12 part of the calculation is more rapidly convergent than the GF2 part likely because of the smallness of the integral value and the partial removal of singularities in the integrand in the former. This fact combined with the flexibility in the mathematical form of the correlation factor make it most advantageous to evaluate the F12 corrections stochastically. Figures 3 and 4 show their convergence with MC steps for C₆₀ and C₇₀, respectively. These plots (displaying only small initial portions of the whole calculations) suggest that the statistical errors in MC-GF2-F12 largely originate from those in MC-GF2, which nevertheless shows stable convergence at reasonably reliable values already after 5 × 10⁶ steps.

FIG. 3.

View large Download slide

The MC-GF2 and MC-GF2-F12 results for the first IE of C₆₀ as a function of the number of MC steps.

FIG. 4.

View large Download slide

The MC-GF2 and MC-GF2-F12 results for the first IE of C₇₀ as a function of the number of MC steps.

The MC-GF2/aug-cc-pVDZ values of the IEs (6.70 eV and 6.73 eV for C₆₀ and C₇₀, respectively) have substantial errors of 0.94 and 0.74 eV from the experimental values. The F12 method improves these results by reducing the errors by 0.47 eV and 0.44 eV, respectively, to 0.47 eV and 0.30 eV from the experimental values. Therefore, the basis-set incompleteness (in GF2/aug-cc-pVDZ) is responsible for 50%–60% of the errors. The remaining errors must be ascribed to the deficiency of the GF2 method itself such as the lack of higher-order perturbation corrections, the diagonal approximation, and/or the frequency-independent approximation, with the first-listed likely being the primary source.⁹ Work to extend MC-GF2 to higher perturbation orders is underway in our laboratory.

V. CONCLUSIONS

We have implemented the MC-GF2-F12 method capable of computing correlated IEs near the CBS limit in a highly scalable manner. It is naturally and easily parallelized with the operation cost for computing all IEs scaling as O(n⁴) with the number (n) of basis functions as opposed to O(n⁵) of its deterministic counterparts. The largest calculation performed in this study (C₇₀) executed the GF2 part on 128 GPUs and the F12 part on 896 CPUs, yielding its first IE near the CBS limit with a statistical error of 0.04 eV. The applicability of these MC calculations is presently only limited by the feasibility of the preceding HF calculations. The three correlation factors that were shown to perform well in MC-MP2-F12 also work equally well in MC-GF2-F12 except that the relative simplicity of STG seems to minimize the statistical errors and is preferred over the YC or SJ factor.

ACKNOWLEDGMENTS

This work was supported by the Center for Scalable, Predictive methods for Excitation and Correlated phenomena (SPEC), which is funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, as a part of the Computational Chemical Sciences Program. It is also part of the Blue Waters sustained-petascale computing project, which is supported by the U.S. National Science Foundation (Award Nos. OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This work was also supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) as Priority Issue on Post-K computer (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use) and JSPS Grant-in-Aids for Scientific Research (A) (Grant No. JP18H03900).

REFERENCES

1.

J.

Linderberg

and

Y.

Öhrn

,

Proc. R. Soc. A

285

,

445

(

1965

).

https://doi.org/10.1098/rspa.1965.0116

Google Scholar

Crossref

2.

L. S.

Cederbaum

and

W.

Domcke

,

Adv. Chem. Phys.

36

,

205

(

1977

).

https://doi.org/10.1002/9780470142554.ch4

Google Scholar

Crossref

3.

P.

Jørgensen

and

J.

Simons

,

Second Quantization-Based Methods in Quantum Chemistry

(

Academic Press

,

New York

,

1981

).

Google Scholar

4.

J.

Oddershede

,

Adv. Chem. Phys.

69

,

201

(

1987

).

https://doi.org/10.1002/9780470142943.ch3

Google Scholar

Crossref

5.

F.

Aryasetiawan

and

O.

Gunnarsson

,

Rep. Prog. Phys.

61

,

237

(

1998

).

https://doi.org/10.1088/0034-4885/61/3/002

Google Scholar

Crossref

6.

J. V.

Ortiz

,

Adv. Quant. Chem.

35

,

33

(

1999

).

https://doi.org/10.1016/s0065-3276(08)60454-2

Google Scholar

Crossref

7.

G.

Onida

,

L.

Reining

, and

A.

Rubio

,

Rev. Mod. Phys.

74

,

601

(

2002

).

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

Google Scholar

Crossref

8.

J. V.

Ortiz

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

3

,

123

(

2013

).

https://doi.org/10.1002/wcms.1116

Google Scholar

Crossref

9.

S.

Hirata

,

A. E.

Doran

,

P. J.

Knowles

, and

J. V.

Ortiz

,

J. Chem. Phys.

147

,

044108

(

2017

).

https://doi.org/10.1063/1.4994837

Google Scholar

Crossref

PubMed

10.

W.

Klopper

,

F. R.

Manby

,

S.

Ten-no

, and

E. F.

Valeev

,

Int. Rev. Phys. Chem.

25

,

427

(

2006

).

https://doi.org/10.1080/01442350600799921

Google Scholar

Crossref

11.

T.

Shiozaki

,

E. F.

Valeev

, and

S.

Hirata

,

Annu. Rep. Comput. Chem.

5

,

131

(

2009

).

https://doi.org/10.1016/s1574-1400(09)00506-4

Google Scholar

Crossref

12.

C.

Hättig

,

W.

Klopper

,

A.

Köhn

, and

D. P.

Tew

,

Chem. Rev.

112

,

4

(

2012

).

https://doi.org/10.1021/cr200168z

Google Scholar

Crossref

PubMed

13.

S.

Ten-no

and

J.

Noga

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

2

,

114

(

2012

).

https://doi.org/10.1002/wcms.68

Google Scholar

Crossref

14.

L.

Kong

,

F. A.

Bischoff

, and

E. F.

Valeev

,

Chem. Rev.

112

,

75

(

2012

).

https://doi.org/10.1021/cr200204r

Google Scholar

Crossref

PubMed

15.

A.

Grüneis

,

S.

Hirata

,

Y. Y.

Ohnishi

, and

S.

Ten-no

,

J. Chem. Phys.

146

,

080901

(

2017

).

https://doi.org/10.1063/1.4976974

Google Scholar

Crossref

PubMed

16.

Q. L.

Ma

and

H.-J.

Werner

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

8

,

e1371

(

2018

).

https://doi.org/10.1002/wcms.1371

Google Scholar

Crossref

17.

Y.-Y.

Ohnishi

and

S.

Ten-no

,

J. Comput. Chem.

37

,

2447

(

2016

).

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

Google Scholar

Crossref

PubMed

18.

F.

Pavošević

,

C.

Peng

,

J. V.

Ortiz

, and

E. F.

Valeev

,

J. Chem. Phys.

147

,

121101

(

2017

).

https://doi.org/10.1063/1.5000916

Google Scholar

Crossref

PubMed

19.

S. Y.

Willow

,

K. S.

Kim

, and

S.

Hirata

,

J. Chem. Phys.

138

,

164111

(

2013

).

https://doi.org/10.1063/1.4801862

Google Scholar

Crossref

PubMed

20.

S. Y.

Willow

,

J.

Zhang

,

E. F.

Valeev

, and

S.

Hirata

,

J. Chem. Phys.

140

,

031101

(

2014

).

https://doi.org/10.1063/1.4862255

Google Scholar

Crossref

PubMed

21.

C. M.

Johnson

,

A. E.

Doran

,

J.

Zhang

,

E. F.

Valeev

, and

S.

Hirata

,

J. Chem. Phys.

145

,

154115

(

2016

).

https://doi.org/10.1063/1.4964854

Google Scholar

Crossref

PubMed

22.

W.

Kutzelnigg

and

J. D.

Morgan

,

J. Chem. Phys.

96

,

4484

(

1992

).

https://doi.org/10.1063/1.462811

Google Scholar

Crossref

23.

E. F.

Valeev

and

C. L.

Janssen

,

J. Chem. Phys.

121

,

1214

(

2004

).

https://doi.org/10.1063/1.1759319

Google Scholar

Crossref

PubMed

24.

H.-J.

Werner

,

T. B.

Adler

, and

F. R.

Manby

,

J. Chem. Phys.

126

,

164102

(

2007

).

https://doi.org/10.1063/1.2712434

Google Scholar

Crossref

PubMed

25.

W.

Kutzelnigg

,

Theor. Chim. Acta

68

,

445

(

1985

).

https://doi.org/10.1007/bf00527669

Google Scholar

Crossref

26.

S.

Ten-no

,

Chem. Phys. Lett.

398

,

56

(

2004

).

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

Google Scholar

Crossref

27.

S.

Ten-no

,

J. Chem. Phys.

121

,

117

(

2004

).

https://doi.org/10.1063/1.1757439

Google Scholar

Crossref

PubMed

28.

T.

Kato

,

Commun. Pure Appl. Math.

10

,

151

(

1957

).

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

Google Scholar

Crossref

29.

R. T.

Pack

and

W.

Byers-Brown

,

J. Chem. Phys.

45

,

556

(

1966

).

https://doi.org/10.1063/1.1727605

Google Scholar

Crossref

30.

K.-C.

Pan

and

H. F.

King

,

J. Chem. Phys.

53

,

4397

(

1970

).

https://doi.org/10.1063/1.1673955

Google Scholar

Crossref

31.

P.

Wind

,

W.

Klopper

, and

T.

Helgaker

,

Theor. Chem. Acc.

107

,

173

(

2002

).

https://doi.org/10.1007/s00214-001-0318-6

Google Scholar

Crossref

32.

A. J.

May

and

F. R.

Manby

,

J. Chem. Phys.

121

,

4479

(

2004

).

https://doi.org/10.1063/1.1780891

Google Scholar

Crossref

PubMed

33.

M.

Torheyden

and

E. F.

Valeev

,

Phys. Chem. Chem. Phys.

10

,

3410

(

2008

).

https://doi.org/10.1039/b803620a

Google Scholar

Crossref

PubMed

34.

A. J.

May

,

E.

Valeev

,

F. R.

Manby

, and

F.

Karlsruhe

,

Phys. Chem. Chem. Phys.

7

,

2710

(

2005

).

https://doi.org/10.1039/b507781h

Google Scholar

Crossref

PubMed

35.

C. M.

Johnson

,

S.

Hirata

, and

S.

Ten-no

,

Chem. Phys. Lett.

683

,

247

(

2017

).

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

Google Scholar

Crossref

36.

A.

Grüneis

,

J. J.

Shepherd

,

A.

Alavi

,

D. P.

Tew

, and

G. H.

Booth

,

J. Chem. Phys.

139

,

084112

(

2013

).

https://doi.org/10.1063/1.4818753

Google Scholar

Crossref

PubMed

37.

A.

Grüneis

,

Phys. Rev. Lett.

115

,

066402

(

2015

).

https://doi.org/10.1103/physrevlett.115.066402

Google Scholar

Crossref

PubMed

38.

R.

Jastrow

,

Phys. Rev.

98

,

1479

(

1955

).

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

Google Scholar

Crossref

39.

M. H.

Kalos

and

P. A.

Whitlock

,

Monte Carlo Methods

(

Wiley-VCH

,

Weinheim

,

2008

).

Google Scholar

Crossref

40.

S. Y.

Willow

,

M. R.

Hermes

,

K. S.

Kim

, and

S.

Hirata

,

J. Chem. Theory Comput.

9

,

4396

(

2013

).

https://doi.org/10.1021/ct400557z

Google Scholar

Crossref

PubMed

41.

S. Y.

Willow

,

K. S.

Kim

, and

S.

Hirata

,

J. Chem. Phys.

137

,

204122

(

2012

).

https://doi.org/10.1063/1.4768697

Google Scholar

Crossref

PubMed

42.

S.

Hirata

,

X.

He

,

M. R.

Hermes

, and

S. Y.

Willow

,

J. Phys. Chem. A

118

,

655

(

2014

).

https://doi.org/10.1021/jp410587b

Google Scholar

Crossref

PubMed

43.

N.

Metropolis

,

A. W.

Rosenbluth

,

M. N.

Rosenbluth

,

A. N.

Teller

, and

E.

Teller

,

J. Chem. Phys.

21

,

1087

(

1953

).

https://doi.org/10.1063/1.1699114

Google Scholar

Crossref

44.

S.

Obara

and

A.

Saika

,

J. Chem. Phys.

84

,

3963

(

1986

).

https://doi.org/10.1063/1.450106

Google Scholar

Crossref

45.

A. E.

Doran

and

S.

Hirata

,

J. Chem. Theory Comput.

12

,

4821

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00588

Google Scholar

Crossref

PubMed

46.

B.

Chan

,

Y.

Kawashima

,

M.

Katouda

,

T.

Nakajima

, and

K.

Hirao

,

J. Am. Chem. Soc.

138

,

1420

(

2016

).

https://doi.org/10.1021/jacs.5b12518

Google Scholar

Crossref

PubMed

47.

H.

Flyvbjerg

and

H. G.

Petersen

,

J. Chem. Phys.

91

,

461

(

1989

).

https://doi.org/10.1063/1.457480

Google Scholar

Crossref

48.

C.

Steffen

,

K.

Thomas

,

U.

Huniar

,

A.

Hellweg

,

O.

Rubner

, and

A.

Schroer

,

J. Comput. Chem.

31

,

2967

(

2010

).

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

Google Scholar

Crossref

PubMed

49.

D. L.

Lichtenberger

,

M. E.

Rempe

, and

S. B.

Gogosha

,

Chem. Phys. Lett.

198

,

454

(

1992

).

https://doi.org/10.1016/0009-2614(92)80027-9

Google Scholar

Crossref

50.

G.

Baumgartner

,

A. A.

Auer

,

D. E.

Bernholdt

,

A.

Bibireata

,

V.

Choppella

,

D.

Cociorva

,

X. Y.

Gao

,

R. J.

Harrison

,

S.

Hirata

,

S.

Krishnamoorthy

,

S.

Krishnan

,

C. C.

Lam

,

Q. D.

Lu

,

M.

Nooijen

,

R. M.

Pitzer

,

J.

Ramanujam

,

P.

Sadayappan

, and

A.

Sibiryakov

,

Proc. IEEE

93

,

276

(

2005

).

https://doi.org/10.1109/jproc.2004.840311

Google Scholar

Crossref

51.

J.

Almlöf

,

K.

Faegri

, and

K.

Korsell

,

J. Comput. Chem.

3

,

385

(

1982

).

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

Google Scholar

Crossref

2018

Author(s)

Monte Carlo explicitly correlated many-body Green’s function theory

I. INTRODUCTION

II. FORMALISMS

A. GF2-F12

B. MC-GF2-F12

III. MONTE CARLO ALGORITHM

IV. RESULTS AND DISCUSSION

V. CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Monte Carlo explicitly correlated many-body Green’s function theory

I. INTRODUCTION

II. FORMALISMS

A. GF2-F12

B. MC-GF2-F12

III. MONTE CARLO ALGORITHM

IV. RESULTS AND DISCUSSION

V. CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only