Hybrid stochastic-deterministic calculation of the second-order perturbative contribution of multireference perturbation theory

TL;DR

This paper introduces a hybrid stochastic-deterministic method for calculating the second-order perturbative energy in multireference perturbation theory, achieving faster error reduction and high accuracy even for large molecular systems.

Contribution

The paper presents a novel hybrid approach that combines stochastic and deterministic calculations to efficiently compute $E^{(2)}$ in MRPT, improving convergence and accuracy over purely stochastic methods.

Findings

Polynomial decay of statistical error ($t^{-3}$ to $t^{-4}$)

Exact results with no error bar possible without stochastic part

Accurate $E^{(2)}$ for large Cr$_2$ molecule with extensive determinants

Abstract

A hybrid stochastic-deterministic approach for computing the second-order perturbative contribution $E^{(2)}$ within multireference perturbation theory (MRPT) is presented. The idea at the heart of our hybrid scheme --- based on a reformulation of $E^{(2)}$ as a sum of elementary contributions associated with each determinant of the MR wave function --- is to split $E^{(2)}$ into a stochastic and a deterministic part. During the simulation, the stochastic part is gradually reduced by dynamically increasing the deterministic part until one reaches the desired accuracy. In sharp contrast with a purely stochastic MC scheme where the error decreases indefinitely as $t^{-1/2}$ (where $t$ is the computational time), the statistical error in our hybrid algorithm displays a polynomial decay $\sim t^{-n}$ with $n=3-4$ in the examples considered here. If desired, the calculation can be carried on until the stochastic part entirely vanishes. In that case, the exact result is obtained with no error bar and no noticeable computational overhead compared to the fully-deterministic calculation. The method is illustrated on the F$_2$ and Cr$_2$ molecules. Even for the largest case corresponding to the Cr$_2$ molecule treated with the cc-pVQZ basis set, very accurate results are obtained for $E^{(2)}$ for an active space of (28e,176o) and a MR wave function including up to $2 \times 10^7$ determinants.

Tables

Table 1. Table 1: Convergence of E ( 2 ) superscript 𝐸 2 E^{(2)} for the \ce Cr2 molecule with bond length 1.68 1.68 1.68  Å  as a function of the wall-clock time for various basis sets (800 CPU cores).

Basis	$E^{\left(2\right)}$	Wall-clock time
cc-pVDZ	$-\mathrm{0.068\hspace{0.17em}3}\left(1\right)$	14 min
	$-\mathrm{0.068\hspace{0.17em}36}\left(1\right)$	55 min
	$-\mathrm{0.068\hspace{0.17em}361}\left(1\right)$	2.4 hr
	$-\mathrm{0.068\hspace{0.17em}360\hspace{0.17em}604}$	3 hr
cc-pVTZ	$-\mathrm{0.124\hspace{0.17em}4}\left(5\right)$	19 min
	$-\mathrm{0.124\hspace{0.17em}7}\left(1\right)$	58 min
	$-\mathrm{0.124\hspace{0.17em}63}\left(1\right)$	3.5 hr
	$-\mathrm{0.124\hspace{0.17em}642}\left(1\right)$	8.7 hr
	—	$\sim$ 15 hr (estimated)
cc-pVQZ	$-\mathrm{0.155\hspace{0.17em}8}\left(5\right)$	56 min
	$-\mathrm{0.155\hspace{0.17em}9}\left(1\right)$	2.5 hr
	$-\mathrm{0.155\hspace{0.17em}95}\left(1\right)$	9.0 hr
	$-\mathrm{0.155\hspace{0.17em}952}\left(1\right)$	18.5 hr
	—	$\sim$ 29 hr (estimated)

Table 2. Table 2: Variational ground-state energy E ( 0 ) superscript 𝐸 0 E^{(0)} and second-order contribution E ( 2 ) superscript 𝐸 2 E^{(2)} of the \ce Cr2 molecule with bond length 1.68 1.68 1.68  Å computed with various basis sets. For all basis sets, the reference is composed of 2 × 10 7 2 superscript 10 7 2\times 10^{7} determinants selected in the valence FCI space (28 electrons).

Reference	Basis	Active space	$E^{\left(0\right)}$	$E^{\left(2\right)}$	$E^{\left(0\right)}+E^{\left(2\right)}$
CIPSI	cc-pVDZ	(28e,76o)	$-\mathrm{2087.227\hspace{0.17em}883\hspace{0.17em}3}$	$-\mathrm{0.068\hspace{0.17em}334}\left(1\right)$	$-\mathrm{2087.296\hspace{0.17em}217}\left(1\right)$
	cc-pVTZ	(28e,126o)	$-\mathrm{2087.449\hspace{0.17em}781\hspace{0.17em}7}$	$-\mathrm{0.124\hspace{0.17em}676}\left(1\right)$	$-\mathrm{2087.574\hspace{0.17em}423}\left(1\right)$
	cc-pVQZ	(28e,176o)	$-\mathrm{2087.513\hspace{0.17em}373\hspace{0.17em}3}$	$-\mathrm{0.155\hspace{0.17em}957}\left(1\right)$	$-\mathrm{2087.669\hspace{0.17em}330}\left(1\right)$

Equations

$$\hat{H}={\hat{H}}^{(0)}+\hat{V},$$

$$\hat{V}=\hat{H}-{\hat{H}}^{(0)}$$

$${\hat{H}}^{(0)}=E^{(0)}\text{outerproduct}\Psi \Psi +\sum _{\alpha \in \mathcal{A}}{H}_{\alpha \alpha }\text{outerproduct}\alpha \alpha ,$$

$$|\Psi ⟩=\sum _{I\in \mathcal{D}}{c}_I|I⟩$$

$$\mathcal{D}=\text{quantity}|I⟩,I=1,\dots ,N_{\text{det}},$$

$$E^{(0)}=\frac{\text{expectationvalue}\hat{H}\Psi }{\text{innerproduct}\Psi \Psi }$$

$$\mathcal{A}=\text{quantity}|\alpha ⟩\notin \mathcal{D}\wedge \text{quantity}(\exists |I⟩\in \mathcal{D}\mid H_{\alpha I}\ne 0).$$

$$E^{(2)}=\sum _{\alpha \in \mathcal{A}}\frac{\text{absolutevalue}\ast {\text{matrixelement}\alpha \hat{H}\Psi }^2}{\Delta E_{\alpha }},$$

$$\mathcal{A}=\bigcup _{I=1}^{N_{\text{det}}}{\mathcal{A}}_I\text{mbox}with{\mathcal{A}}_I\cap {\mathcal{A}}_J=\varnothing \text{mbox}ifI\ne J.$$

$$w_I=\frac{c_I^2}{\text{innerproduct}\Psi \Psi }.$$

$${\mathcal{A}}_I=\text{quantity}|\alpha ⟩\in \mathcal{A}\mid H_{\alpha I}\ne 0\wedge (\forall J<I,|\alpha ⟩\notin {\mathcal{A}}_J).$$

$$E^{(2)}=\sum _{I=1}^{N_{\text{det}}}{e}_I,$$

$$e_I=\sum _{\alpha \in {\mathcal{A}}_I}\frac{\text{absolutevalue}\ast {\text{matrixelement}\alpha \hat{H}\Psi }^2}{\Delta E_{\alpha }}.$$

$$e_I=\sum _{\alpha \in {\mathcal{A}}_I}\frac{\text{absolutevalue}\ast {\text{matrixelement}\alpha \hat{H}{\Psi }_I}^2}{\Delta E_{\alpha }},$$

$$|{\Psi }_I⟩=\sum _{J=I}^{N_{\text{det}}}{c}_J|J⟩$$

$$E^{(2)}=\sum _{I=1}^{N_{\text{det}}}{e}_I=\sum _{I=1}^{N_{\text{det}}}\sum _{\alpha \in {\mathcal{A}}_I}\frac{\text{absolutevalue}\ast {\text{matrixelement}\alpha \hat{H}{\Psi }_I}^2}{\Delta E_{\alpha }}.$$

$$E^{(2)}=\sum _{I=1}^{N_{\text{det}}}{p}_I\text{quantity}(\frac{e_I}{p_I}),$$

$$E^{(2)}=\text{expectationvalue}{\frac{e_I}{p_I}}_{p_I}.$$

$$p_I^{\text{opt}}=\frac{e_I}{E^{(2)}}.$$

$$p_I=\frac{\text{innerproduct}{\Psi }_I{\Psi }_I}{{\sum }_{J=1}^{N_{\text{det}}}\text{innerproduct}{\Psi }_J{\Psi }_J}=\frac{{\sum }_{J=I}^{N_{\text{det}}}{c}_J^2}{{\sum }_{J=1}^{N_{\text{det}}}{\sum }_{K=J}^{N_{\text{det}}}{c}_K^2}.$$

$$p_I=\frac{c_I^2}{{\sum }_{J=1}^{N_{\text{det}}}{c}_J^2}=w_I.$$

$$R(I-1)\le u\le R(I),$$

$$R(I)=\sum _{J=1}^I{p}_J,$$

$$u_k=\frac{k-1+u}{M}\text{mbox}fork=1,\dots ,M,$$

$$R(I_k-1)\le u_k\le R(I_k)\text{mbox}fork=1,\dots ,M$$

$$\mathcal{D}=\bigcup _{k=1}^M{\mathcal{D}}_k\text{mbox}with{\mathcal{D}}_k\cap {\mathcal{D}}_l=\varnothing ,\forall k\ne l,$$

$$E^{(2)}=\sum _{k=1}^M\sum _{I_k\in {\mathcal{D}}_k}{e}_{I_k}.$$

$$E^{(2)}=\text{expectationvalue}{\frac{1}{M}\sum _{k=1}^M\frac{e_{I_k}}{p_{I_k}}}_{p(I_1,\dots ,I_M)},$$

$$\sum _{I_1}\dots \sum _{I_{k-1}}\sum _{I_{k+1}}\dots \sum _{I_M}p(I_1,\dots ,I_M)=M{p}_{I_k},$$

$$\sum _{I_k\in {\mathcal{D}}_k}{p}_{I_k}=\frac{1}{M}.$$

Full text

I Introduction

Multireference (MR) approaches are based upon the distinction between non-dynamical (or static) and dynamical correlation effects. Though such a clear-cut distinction is questionable, it is convenient to discriminate between the so-called static correlation effects emerging whenever the description of the molecular system using a single configuration breaks down (excited-states, transition-metal compounds, systems far from their equilibrium geometry, …)Szalay et al. (2012) and the dynamical correlation effects resulting from the short-range part of the electron-electron repulsion.Hättig et al. (2012)

To quantitatively establish this distinction, the Hamiltonian is decomposed as

[TABLE]

where the zeroth-order Hamiltonian $\hat{H}^{(0)}$ is chosen in conjunction with some MR wave function including the most chemically-relevant configurations at the origin of static correlation effects, and

[TABLE]

is the residual part describing the bulk of dynamical correlation effects. The plethora of MR methods found in the literature results from the large freedom in choosing $\hat{H}^{(0)}$, and the fact that $\hat{V}$ may or may not be treated perturbatively. Among the non-perturbative approaches, let us cite the two most common ones, namely the MR configuration interaction (MRCI)Lischka et al. (1981); Szalay et al. (2011); Lischka et al. (2001) and the MR coupled cluster (MRCC)Paldus and Li (1999); Kong et al. (2009); Jeziorski (2010); Lyakh et al. (2012) approaches. However, because of their high computational cost, these methods are usually limited to systems of moderate size.

To overcome the computational burden associated with these methods — yet still capturing the main physical effects — a natural idea is to treat the potential as a perturbation, entering the realm of MR perturbation theories (MRPT). Several flavors of MRPT exist depending on the choice of $\hat{H}^{(0)}$ (Epstein-Nesbet decomposition,Epstein (1926); Nesbet (1955) Dyall Hamiltonian,Angeli et al. (2001); Angeli, Cimiraglia, and Malrieu (2002) Fink’s partitioning,Fink (2006, 2009) etc). Among the most commonly-used approaches, we have the CASPT2Andersson et al. (1990); Andersson, Malmqvist, and Roos (1992) and NEVPT2Angeli et al. (2001); Angeli, Cimiraglia, and Malrieu (2002) methods. Regarding the construction of the zeroth-order part, CASSCF-type approaches are the most widely-used schemes,Siegbahn et al. (1980); Roos, Taylor, and Siegbahn (1980); Siegbahn et al. (1981) but other methods, such as CASCI, selected CI (see Refs. Huron, Rancurel, and Malrieu, 1973; Caffarel et al., 2016 and references therein), FCIQMCBooth, Thom, and Alavi (2009); Cleland et al. (2012); Petruzielo et al. (2012) or DMRG-type approachesWhite (1992, 1993); Sharma and Chan (2012) can also be employed.

In this work, we shall consider MRPT theories limited to the second order in perturbation (MRPT2).Andersson et al. (1990) We address the important problem of calculating efficiently the second-order perturbative contribution $E^{(2)}$ in situations where standard calculations become challenging. Here, we suppose that the MR wave function has already been constructed by any method of choice.

Although the present method can be easily generalized to any externally-decontracted MRPT approach (such as the recently-introduced JM-MRPT2 method Giner et al. (2017)), for the sake of simplicity, we shall restrict ourselves here to MR Epstein-Nesbet perturbation theory. Extension to externally-contracted methods, such as CASPT2 or NEVPT2, is less obvious — although not impossible — since the excited contracted wave functions are non-orthogonal.

The computational cost of MRPT2 can rapidly become unbearable when the number of electrons $N_{\text{el}}$ and the number of one-electron basis functions $N_{\text{bas}}$ become large. The cost is indeed proportional to the number of reference determinants $N_{\text{det}}$ times the total number of singly- and doubly-excited determinants (scaling as $N_{\text{el}}^{2}N_{\text{bas}}^{2}$). Because our main goal is to treat large, chemically-relevant systems, the development of fast and accurate schemes for computing $E^{(2)}$ becomes paramount. Of course, in actual calculations, a trade-off must be found between the price to pay to build the MR wave function and the effort needed to evaluate $E^{(2)}$. Increasing $N_{\text{det}}$ (i.e. improving the MR wave function) may appear as the natural thing to do as the magnitude of $E^{(2)}$ decreases and the contribution of the neglected higher orders is made smaller. However, its computational price (proportional to $N_{\text{det}}$) increases stiffly and the calculation becomes rapidly unfeasible. Unfortunately, this balance is strongly dependent on the method used to generate the MR wave function and on the ability to compute rapidly and accurately $E^{(2)}$.

In this work, we present a simple and efficient Monte Carlo (MC) method for computing the second-order perturbative contribution $E^{(2)}$. For all the systems reported here, the reference space is constructed using the CIPSI method,Huron, Rancurel, and Malrieu (1973); Evangelisti, Daudey, and Malrieu (1983); Caffarel et al. (2016) a selected CI approach where important determinants are selected perturbatively. However, other variants of selected CI approaches or any other method for constructing the reference wave function may, of course, be used. Note that, in this study, the reported wall-clock times only refer to the computation of $E^{(2)}$, i.e. do not take into account the preliminary calculation of the reference wave function.

A natural idea to evaluate $E^{(2)}$ with some targeted accuracy is to truncate the perturbational sum over excited determinants. However, since all the terms of the second-order sum have the same (negative) sign, the truncation will inevitably introduce a bias which is difficult to control. A way to circumvent this problem is to resort to a stochastic sampling of the various contributions. In this case, the systematic bias is removed at the price of introducing a statistical error. The key property is that this error can now be controlled thanks to the central-limit theorem. However, in practice, to make the statistical average converge rapidly and to get statistical error small enough, care has to be taken in the way the statistical estimator is built and how the sampling is performed. The purpose of the present work is to propose a practical solution to this problem.

Note that the proposal of computing stochastically perturbative contributions is not new. In the context of second-order Møller-Plesset (MP2) theory, where the reference Hamiltonian reduces to the Hartree-Fock Hamiltonian, Hirata and coworkers have proposed a MC scheme for calculating the MP2 correlation energy.Willow, Kim, and Hirata (2012); Willow et al. (2014) However, we point out that this approach, based on a single-reference wave function, samples a 13-dimensional integral (in time and space) and has no direct relation with the present method. In a recent work, Sharma et al.Sharma et al. (2017) address the very same problem of computing stochastically the second-order perturbative contribution of Epstein-Nesbet MRPT. Similarly to what is proposed here, $E^{(2)}$ is recast as a sum over contributions associated with each reference determinant and contributions are stochastically sampled. However, the definition of the quantities to be averaged and the way the sampling is performed are totally different. Finally, let us mention the recent work of Jeanmairet et al.Jeanmairet, Sharma, and Alavi (2017) addressing a similar problem in a different way. Within the framework of the recently proposed linear CC MRPT, it is shown that both the zeroth-order and first-order wave functions can be sampled using a generalization of the FCIQMC approach. Here also, $E^{(2)}$ can be expressed as a stochastic average.

The present paper is organized as follows. In Sec. II, we report notations and basic definitions for MRPT2. Section III proposes an original reformulation of the second-order contribution allowing an efficient MC sampling. The expression of the MC estimator is given and a hybrid stochastic-deterministic approach greatly reducing the statistical fluctuations is presented. In Sec. IV, some illustrative applications for the \ceF2 and \ceCr2 molecules are discussed. Finally, some concluding remarks are given in Sec. V.

II

Second-order multireference perturbation theory

II.1 Second-order energy contribution

In MR Epstein-Nesbet perturbation theory, the reference Hamiltonian is chosen to be

[TABLE]

where $H_{\alpha\alpha}=\expectationvalue{\hat{H}}{\alpha}$ and

[TABLE]

is the reference wave function expressed as a sum of $N_{\text{det}}$ determinants belonging to the reference space

[TABLE]

and

[TABLE]

is the corresponding (variational) energy. The sum in Eq. (3) is over the set of determinants $\ket{\alpha}$ that do not belong to $\mathcal{D}$ but are connected to $\mathcal{D}$ via $\hat{H}$:

[TABLE]

Due to the two-body character of the interaction, the determinants $\ket{\alpha}$ are either singly- or doubly-excited with respect to (at least) one reference determinant.Szabo and Ostlund (1989) However, several reference determinants can be connected to the same $\ket{\alpha}$.

Using such notations, the second-order perturbative contribution is written as

[TABLE]

with $\Delta E_{\alpha}=E^{(0)}-H_{\alpha\alpha}$.

II.2 Partition of $\mathcal{A}$

The first step of the method — instrumental in the MC algorithm efficiency — is the partition of $\mathcal{A}$ into $N_{\text{det}}$ subsets $\mathcal{A}_{I}$ associated with each reference determinant $\ket{I}$:

[TABLE]

To define ${\cal A}_{I}$, the determinants $\ket{I}$ are first sorted in descending order according to the weight

[TABLE]

The partition of $\mathcal{A}$ starts with $\mathcal{A}_{1}$ defined as the set of determinants $\ket{\alpha}\in\mathcal{A}$ connected to the first reference determinant (i.e. $I=1$). Then, $\mathcal{A}_{2}$ is constructed as the set of determinants of $\mathcal{A}$ connected to the determinant corresponding to $I=2$, but not belonging to $\mathcal{A}_{1}$. The process is carried on up to the last determinant. This partition is schematically illustrated in Fig. 1. Mathematically, it can be written as

[TABLE]

Because of the way they are constructed, the size of $\mathcal{A}_{I}$ is expected to decrease rapidly as a function of $I$, except for a possible transient regime for very small $I$.

A key point in the construction of the partition of $\mathcal{A}$ is to avoid both the computation of redundant contributions and the storage of unnecessary intermediates. First, when a determinant $\ket{\alpha}$ is generated by applying a single or double excitation operator to a reference determinant $\ket{I}$, one has to check that $\ket{\alpha}$ does not belong to $\mathcal{D}$. If the reference determinants are stored in a hash table, the presence of $\ket{\alpha}$ in $\mathcal{D}$ can be checked in constant time. Next, one has to know if $\ket{\alpha}$ has already been generated via another reference determinant $\ket{J}$. To do so, one must compute the number of holes and particles between $\ket{\alpha}$ and each determinant preceding $\ket{I}$ in $\mathcal{D}$. As soon as an excitation degree lower than 3 is found, the search can be aborted since the contribution is known to have been considered before. In the worst-case scenario, this step scales as $\order{N_{\text{det}}}$, and the prefactor is very small since finding the excitation degree between two determinants can be performed in less than 20 CPU cyclesScemama and Giner (2013) (comparable to a floating-point division). Furthermore, the asymptotic scaling can be further reduced by sorting the determinants in groups with same spin string. Indeed, one only has to probe determinants $\ket{J}$ that are no more than quadruply excited with respect to $\ket{I}$, and, if the search is restricted to groups with spin-up string, the asymptotic scaling reduces to $\order{\sqrt{N_{\text{det}}}}$. To provide a quantitative illustration of the computational effort associated with the construction of the partitioning, using $2\times 10^{7}$ determinants (as in the case of \ceCr2 presented below), this preliminary step is negligible: on a single 2.7 GHz core, the calculation takes $20\text{ cycles }\times N_{\text{det}}^{3/2}/(2.7\times 10^{9}\text{ cycles/s})\sim 663\text{ s}$ (CPU time), while the total execution time (wall-clock time) of the entire run ranges from $14$ minutes to $18.5$ hr using 800 cores (see Table 1).

II.3 Partition of $E^{(2)}$

Thanks to the partition of $\mathcal{A}$ (see Eq. (11)), the sum (8) can be decomposed into a sum over the reference determinants $\ket{I}$:

[TABLE]

where

[TABLE]

Moreover, noticing that by construction the determinants $\ket{\alpha}$ belonging to $\mathcal{A}_{I}$ are not connected to the part of the reference function expanded over the preceding reference determinants, we have

[TABLE]

where

[TABLE]

is a truncated reference wave function. Our final working expression for the second-order contribution $E^{(2)}$ is thus written as

[TABLE]

A key property at the origin of the efficiency of the MC simulations presented below is that the $e_{I}$’s take their largest values at very small $I$. Then, they decay very rapidly as $I$ increases.

This important property is illustrated in Fig. 2. The data have been obtained for the \ceF2 molecule at the equilibrium bond length of $R_{\text{F-F}}=1.4119$ Å  using Dunning’s cc-pVQZ basis set.Dunning (1989) The multideterminant reference space is built by selecting determinants using the CIPSI algorithm. Figure 2 displays the $e_{I}$’s for the first 20 000 selected determinants. As one can see, the $e_{I}$’s decay very rapidly with $I$. Of course, at the scale of the individual determinant, there is no guarantee of a strictly monotonic decay, and it is indeed what we observe. By averaging groups of successive $e_{I}$’s, the curve can be smoothed out. The two data sets presented in Fig. 2 have been obtained by averaging either by groups of 20 (point cloud) or 100 (solid line) values.

It is important to note that the very rapid decay of the $e_{I}$’s is a direct consequence of the way we have chosen to decompose $\mathcal{A}$. To be more precise, we note that in Eq. (14), the decay has three different origins:

• •

the number of determinants involved in the sum over $\ket{\alpha}$ decreases as a function of $I$;

• •

the excitation energies $\Delta E_{\alpha}$ increase with $I$;

• •

the norm of the truncated wave function $\Psi_{I}$ decreases rapidly (as $c^{2}_{I}$) when $I$ increases.

In addition, as a consequence of the first point, we note that the computation of $e_{I}$ becomes faster when $I$ increases.

III

Monte Carlo method

III.1 Monte Carlo estimator

To get an expression of $E^{(2)}$ suitable for MC simulations, the second-order contribution is recast as

[TABLE]

and is thus rewritten as the following MC estimator

[TABLE]

Here, $p_{I}$ is an arbitrary probability distribution. The optimal choice for $p_{I}$ is given by the zero-variance condition, i.e.

[TABLE]

Note that $e_{I}$ and $E^{(2)}$ being both negative, the probability distribution $p_{I}$ is positive, as it should be.

To build a reasonable approximation of $p_{I}$, we note that the magnitude of $e_{I}$, as expressed in Eq. (14), is essentially given by the norm of the truncated wave function $\Psi_{I}$ (see Eq. (15)). Thus, a natural choice for the probability distribution is

[TABLE]

In our simulations, we have observed that summing totally or partially the squared coefficients in the numerator does not change significantly the statistical fluctuations. As a consequence, we restrict the summation in Eq. (20) to the leading term, i.e.

[TABLE]

Let us emphasize that performing a MC simulation in the $e_{I}$ space is highly beneficial since the number of $e_{I}$ is always small enough to make them all fit in memory. Hence, one can follow the so-called lazy evaluation strategy:Hudak (1989) the value of $e_{I}$ is computed only once when needed for the first time, and its value is then stored. If the same $e_{I}$ is requested later, the stored value will be returned.

III.2 Improved Monte Carlo sampling

The stochastic calculation of $E^{(2)}$, Eq. (18), can be done in a standard way by sampling the probability distribution and averaging the successive values of $e_{I}/p_{I}$. In practice, the sampling can be realized by drawing, at each MC step, a uniform random number $u\in\quantity[0,1]$ and selecting the determinant $\ket{I}$ verifying

[TABLE]

where $R$ is the cumulative distribution function of the probability distribution defined as

[TABLE]

with $R(0)=0$.

At this stage, it is useful to take advantage of the fact that, thanks to the way the $e_{I}$’s have been constructed, the quantity to be averaged, $e_{I}/p_{I}$, is a slowly-varying function of $I$ (providing that the small-scale fluctuations present at the level of individual determinants have been averaged out). This property, which is well illustrated by Fig. 2, is shared by $p_{I}\sim c^{2}_{I}$, hence by the ratio $e_{I}/p_{I}$. Thus, an efficient way to reduce the statistical fluctuations consists in sampling piece-wisely $\mathcal{D}$ by decomposing it into subdomains where the integrand is a slowly-varying function (see the justification of this statement after Eq. (30)).

To implement this idea, the interval $[0,1]$ is divided into $M$ equally-spaced intervals $\mathcal{U}_{k}$ and a “comb” of correlated random numbers

[TABLE]

covering uniformly $[0,1]$ is created (where $u$ is a single uniform random number). At each MC step, a $M$-tuple of determinants $(I_{1},I_{2},\ldots,I_{M})$ verifying

[TABLE]

is drawn.

Defining $\mathcal{D}_{k}$ as the subset of determinants $\ket{I_{k}}$ satisfying $R(I_{k})\in\mathcal{U}_{k}$, we introduce the following partition

[TABLE]

and express $E^{(2)}$ as a sum of $M$ contributions associated with each $\mathcal{D}_{k}$:

[TABLE]

Using the process described above (Eqs. (24) and  (25)), the second-order energy can be rewritten as the following MC estimator

[TABLE]

where $p(I_{1},\ldots,I_{M})$ denotes the normalized probability distribution corresponding to Eqs. (24) and (25). Equation (28) follows from the fact that, by construction, $p_{I_{k}}$ is the $k$-th marginal distribution of $p(I_{1},\ldots,I_{M})$

[TABLE]

with

[TABLE]

By drawing determinants on separate subsets $\mathcal{D}_{k}$, the sum to be averaged in Eq. (28), is expected to fluctuate less than the very same sum computed by independently drawing determinants over $\mathcal{D}$. This remarkable property can be explained as follows. For large $M$, the fluctuations of the sum based on independent drawings behave as in any MC scheme, i.e. as $M^{-1/2}$. Using a comb covering evenly (with weight $p_{I}$) the determinant space, the situation is different since the sum can now be seen as a Riemann sum over $\mathcal{D}$ with a residual error behaving as $M^{-1}$. As a consequence, the overall reduction in statistical noise resulting from the use of the comb is expected to be of the order of $\sqrt{M}$. We emphasize that such an attractive feature is only observed because $e_{I}/p_{I}$ is a slowly-varying function of $I$ (as mentioned above). In the opposite case, the gain would vanish. In the application on the \ceF2 molecule presented below (see Fig. 5), the numerical results confirm this: a decrease of about one order of magnitude in statistical error is obtained when using $M$=100. Note that using a comb reduces the estimator’s variance, but does not change the typical inverse square root behavior of the statistical error with respect to the number of MC steps.

Note that Eq. (26) is actually not correct when some determinants (first and/or last determinant of a given subset) belong to more than one subset. Thus, special care has to be taken for determinants at the boundary of two subsets, but this difficulty can be easily circumvented by formally duplicating each of these determinants into copies with suitable weights.

III.3 Hybrid stochastic-deterministic scheme

In practice, because the first few determinants are responsible for the most significant contribution in Eq. (17), it is advantageous not to sample the entire reference space but to remove from the stochastic sampling the leading determinants. Consequently, $E^{(2)}$ is split into a deterministic $E_{\text{D}}^{(2)}$ and a stochastic $E_{\text{S}}^{(2)}$ component, such as

[TABLE]

where $\mathcal{D}_{\text{D}}$ is the set of determinants in the deterministic space, and $\mathcal{D}_{\text{S}}=\mathcal{D}\setminus\mathcal{D}_{\text{D}}$ is its stochastic counterpart.

At a given point of the simulation, some determinants have been drawn, and some have not. If we keep track of the list of the drawn determinants, we can check periodically, for each $\mathcal{D}_{k}$, whether or not all elements have been drawn at least once. If that is the case, the full set of determinants is moved to $\mathcal{D}_{\text{D}}$ and the corresponding contribution $\sum_{I_{k}}e_{I_{k}}$ is added to $E_{\text{D}}^{(2)}$. The statistical average and error bar are then updated accordingly. The expression of the $E^{(2)}$ estimator is now time-dependent and, at the $m$th MC step, the deterministic part is given by

[TABLE]

where

[TABLE]

On the other hand, the stochastic part is now given by

[TABLE]

where

[TABLE]

and $n_{I_{k}}^{(m)}$ denotes the number of times the determinant $I_{k}$ has been drawn at iteration $m$.

If desired, the calculation can be carried on until the stochastic part entirely vanishes. In that case, all the determinants are in $\mathcal{D}_{\text{D}}$ and the exact value of $E^{(2)}$ is obtained with zero statistical fluctuations.

Finally, to make sure that a given set $\mathcal{D}_{k}$ does not stay in the stochastic part because a very small number of its determinants has not been drawn, we have implemented an additional step as follows. At each MC iteration (where a new comb is created), the contribution $e_{I}$ of the first not-yet-sampled determinant (i.e. corresponding to the smallest $I$ value in the sorted determinant list) is calculated and stored. By doing this, the convergence of the hybrid stochastic-deterministic estimator is significantly improved. Moreover, after $N_{\text{det}}$ MC steps, it is now guaranteed that the exact deterministic value is reached.

III.4 Upper bound on the computational time

In the present method, the vast majority of the computational time is spent calculating $e_{I}$’s. A crucial point which makes the algorithm particularly efficient is the lazy evaluation of these quantities. This implies that, in practice, the stochastic calculation will never be longer than the time needed to compute all the individual $e_{I}$’s (i.e. the time necessary to complete the fully-deterministic calculation) due to the negligible time required by the MC sampling (drawing 100 million random numbers takes less than 3 seconds on a single CPU core).

Finally, it is noteworthy that the final expression of $E^{(2)}$ can be very easily decomposed into (strictly) independent calculations. The algorithm presented here is thus embarrassingly parallel (see Sec. IV.3 below).

IV

Numerical tests

The present algorithm has been implemented in our Quantum Package code.Scemama et al. (2015) The perturbatively-selected CI algorithm CIPSI,Huron, Rancurel, and Malrieu (1973); Evangelisti, Daudey, and Malrieu (1983) as described in Ref. Caffarel et al., 2016, is used to build the multideterminant reference space. In all the calculations performed in this section, we have chosen to use a comb with $M=100$. All the simulations were performed on the Curie supercomputer (TGCC/CEA/GENCI) where each node is a dual socket Xeon E5-2680 @ 2.70GHz with 64GiB of RAM, interconnected with an Infiniband QDR network.

IV.1 F2 molecule

As a first illustrative example, we consider the calculation of $E^{(2)}$ for the \ceF2 molecule in its ${}^{1}\Sigma_{g}^{+}$ electronic ground state at equilibrium geometry. The two $1s$ core electrons are kept frozen and the Dunning’s cc-pVQZ basis set is used. The Hilbert space is built by distributing the 14 active electrons within the 108 non-frozen molecular orbitals for a total of more than $10^{20}$ determinants.

Despite the huge size of the Hilbert space, the selected CI approach is able to reach the full CI (FCI) limit with a very good accuracy. The convergence of the variational energy $E^{(0)}$ and of the total energy (given by the sum of the variational and second-order contribution $E^{(0)}+E^{(2)}$) with respect to the number of selected determinants are presented in Fig. 3. The maximum number of determinants we have selected is $4\times 10^{6}$. For this value, $E^{(0)}$ is not converged but is already a reasonable approximation to the FCI energy with an error of about 18 m$E_{\text{h}}$. In sharp contrast, the total energy including the second-order correction converges very rapidly: millihartree accuracy is reached with about $2\times 10^{6}$ determinants. For $N_{\text{det}}=4\times 10^{6}$, the best value obtained is $-199.3594$ a.u., in quantitative agreement with the estimated FCI value of $-199.3598(2)$ a.u. obtained by Cleland et al. with FCIQMC.Cleland et al. (2012)

For this system and the maximum number of selected determinants considered, it is actually possible to calculate exactly $E^{(2)}$ by explicit evaluation of the entire sum (deterministic method). The corresponding wall-clock times (in minutes) using 50 nodes (800 cores) are reported directly in Fig. 3. For $N_{\text{det}}=10^{4}$, the calculation takes a few seconds while for the largest number of $N_{\text{det}}=4\times 10^{6}$ about 35 minutes are needed.

We now consider the hybrid stochastic-deterministic evaluation of $E^{(2)}$. The left graph of Fig. 4 shows the evolution of $E^{(2)}$ as a function of the wall-clock time (in minutes). Data are given for the cc-pVQZ basis and $N_{\text{det}}=4\times 10^{6}$. Similar curves are obtained with the two other basis sets. As one can see, the rate of convergence of the error is striking and, eventually, the exact value is obtained with very small fluctuations. If chemical accuracy is targeted (error of roughly 1 m$E_{\text{h}}$), three minutes are needed using 800 cores. This value has to be compared with the $\sim 35$ minutes needed to evaluate the exact value (see Fig. 3).

To have a better look at the fluctuations, the statistical error as a function of the wall-clock time is reported in Fig. 5. We have reported four curves to show the effects of the different strategies used in our algorithm. The first one (in green) is the curve one would typically obtain using a standard MC algorithm where the contributions are always recomputed (no lazy evaluation). Note that, for this particular curve, we have not performed the calculation, but we have plotted an arbitrary $\sigma_{1}/\sqrt{t}$ curve to illustrate its decay rate. The light blue curve is obtained using the MC estimator proposed in Sec. III.1. The slope is steeper than for the standard MC scheme thanks to the lazy evaluation strategy. The introduction of the comb (Sec. III.2) reduces the statistical error by an order of magnitude, and produces the dark blue curve. Finally, incorporating the hybrid deterministic/stochastic scheme (Sec. III.3) yields the red curve.

Quite remarkably, the overall convergence of the red curve is extremely rapid. Because of the irregular convergence, it is not easy to extract the exact mathematical form of the decay. However, it is clear that a typical polynomial decay is observed. Fitting the curve of the hybrid scheme gives a decrease of the error bar between $t^{-3.1}$ and $t^{-3.6}$, which is significantly faster than the $t^{-1/2}$ behavior of the standard MC algorithm. Note also that some discontinuities in the statistical error are regularly observed. Such sudden drops occur each time a subset $\mathcal{D}_{k}$ is entirely filled and its contribution is transferred to the deterministic part. Comparison with standard MC illustrates that obtaining an arbitrary accuracy with a standard MC sampling can rapidly become prohibitively expensive. Most importantly, the wall-clock time would rapidly become larger than the price to pay to compute exactly (i.e. deterministically) $E^{(2)}$, which is not the case with the here-proposed method.

IV.2 \ceCr2 molecule

We now consider the challenging example of the \ceCr2 molecule in its ${}^{1}\Sigma_{\text{g}}^{+}$ ground state. The internuclear distance is chosen to be close to its experimental equilibrium geometry, i.e. $R_{\text{Cr-Cr}}$ = 1.68 Å. Full-valence calculations including 28 active electrons (two frozen neon cores) are performed. The cc-pVDZ, TZ and QZ basis setsBalabanov and Peterson (2005) are employed and the associated active spaces corresponding to (28e,76o), (28e,126o) and (28e,176o) include more than $10^{29}$, $10^{36}$ and $10^{42}$ determinants, respectively. For all the basis sets, the molecular orbitals (MOs) were obtained with the GAMESS Schmidt et al. (1993) program using a CASSCF calculation with 12 electrons in 12 orbitals, and $2\times 10^{7}$ determinants were selected in the FCI space with the CIPSI algorithm implemented in Quantum Package. In the cc-pVQZ basis set, we had to remove the $h$ functions of the basis set since the version of GAMESS we used (prior to 2013) does not handle the corresponding two-electron integrals.

The right graph of Fig. 4 shows the convergence of $E^{(2)}$ as a function of the wall-clock time for the cc-pVTZ basis set and $N_{\text{det}}=2\times 10^{7}$. Again, similar curves are obtained with the two other basis sets. Similarly to \ceF2, the convergence is remarkably fast with a steep decrease of the statistical error with respect to the wall-clock time (for quantitative results, see Table 1). Note that the maximum energy range in the right graph of Fig. 4 is only 0.35 m$E_{\text{h}}$.

Table 2 reports the quantitative results obtained with the three basis sets. One can see that very accurate results for $E^{(2)}$ can be obtained even with the largest QZ basis set. For the three basis sets, the statistical error obtained is $10^{-6}$ $E_{\text{h}}$. However, it is clear that in practical applications we do not need such high level of accuracy as the finite-size basis effects as well as the high-order perturbative contributions are much larger. If, more reasonably, we target an accuracy of about 0.1 m$E_{\text{h}}$, we see in Table 1 that the wall-clock time needed is about 14 min, 1 hr and 2.5 hr with 800 CPU cores for the DZ, TZ, and QZ basis sets, respectively. Finally, we note that, in contrast with \ceF2, the absolute value of $E^{(2)}$ remains large even when relatively large MR wave functions are employed. This result clearly reflects the difficulty in treating accurately \ceCr2. We postpone to a forthcoming paper the detailed analysis of this system and the calculation of the entire potential energy curve.

IV.3 Parallel speedup

To measure the parallel speedup of the present implementation of our algorithm, we have measured the wall-clock time needed to reach a target statistical error of $10^{-6}$ a.u. with $800$, $4\,000$ and $8\,000$ cores ($50$, $250$ and $500$ nodes) using the \ceCr2/cc-pVQZ wave function with $N_{\text{det}}=2\times 10^{7}$. The speedup is calculated using the 800-core run as the reference, and the results are shown in Fig. 6. Going from $800$ to $4\,000$ cores gives a speedup of $4.95$, and the $8\,000$-core run exhibits a speedup of $9.82$. These values are extremely close to the ideal values of $5$ and $10$. Therefore, we believe that this method is a good candidate for running on exascale machines in a near future.

V

Conclusions

In this work, a hybrid stochastic-deterministic algorithm to compute the second-order energy contribution $E^{(2)}$ within Epstein-Nesbet MRPT has been introduced. Two main ideas are at the heart of the method. First, the reformulation of the standard expression of $E^{(2)}$, Eq. (8), into Eq. (16). Thanks to the unique property of the elementary contributions $e_{I}$, the latter expression turns out to be particularly well-suited for low-variance MC calculations. The second idea, which greatly enhances the convergence of the calculation, is to decompose $E^{(2)}$ as a sum of a deterministic and a stochastic part, the deterministic part being dynamically updated during the calculation.

We have observed that the size of the stochastic part (as well as the statistical error) decays in time with a polynomial behavior. If desired, the calculation can be carried on until the stochastic part entirely vanishes. In that case, the exact (deterministic) result is obtained with no error bar and no noticeable computational overhead compared to the fully-deterministic calculation. Such a remarkable result is in sharp contrast with standard MC calculations where the statistical error decreases indefinitely as the inverse square root of the simulation time.

The numerical applications presented for the \ceF2 and \ceCr2 molecules illustrate the great efficiency of the method. The largest calculation on \ceCr2 (cc-pVQZ basis set) has an active space of (28e,176o), corresponding to a Hilbert space consisting of approximately $10^{42}$ determinants and a multireference wave function containing $2\times 10^{7}$ determinants. Even in this extreme case, $E^{(2)}$ can easily be calculated with sub-millihartree accuracy using a fully- and massively-parallel version of the algorithm.

As a final comment, we would like to mention that, although we have only considered two illustrative examples in the present manuscript, our method has been shown to be highly successful in all the cases we have considered.

Acknowledgements.

We thank the referees for their valuable comments on the first version of our manuscript. This work was performed using HPC resources from CALMIP (Toulouse) under allocation 2016–0510 and from GENCI-TGCC (Grant 2016–08s015).