Zero-variance zero-bias quantum Monte Carlo estimators of the spherically and system averaged pair density

# Zero-variance zero-bias quantum Monte Carlo estimators of the spherically and system averaged pair density

Julien Toulouse Cornell Theory Center, Cornell University, Ithaca, New York 14853, USA.    Roland Assaraf Laboratoire de Chimie Théorique, Université Pierre et Marie Curie and Centre National de la Recherche Scientifique, 75005 Paris, France.    C. J. Umrigar Cornell Theory Center and Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New York 14853, USA.
July 9, 2019
###### Abstract

We construct improved quantum Monte Carlo estimators for the spherically- and system-averaged electron pair density (i.e. the probability density of finding two electrons separated by a relative distance ), also known as the spherically-averaged electron position intracule density , using the general zero-variance zero-bias principle for observables, introduced by Assaraf and Caffarel. The calculation of is made vastly more efficient by replacing the average of the local delta-function operator by the average of a smooth non-local operator that has several orders of magnitude smaller variance. These new estimators also reduce the systematic error (or bias) of the intracule density due to the approximate trial wave function. Used in combination with the optimization of an increasing number of parameters in trial Jastrow-Slater wave functions, they allow one to obtain well converged correlated intracule densities for atoms and molecules. These ideas can be applied to calculating any pair-correlation function in classical or quantum Monte Carlo calculations.

## I Introduction

Two-electron distribution functions occupy an important place in electronic structure theory between the simplicity of one-electron densities and the complexity of the many-electron wave function. In particular, the system-averaged electron pair density, the probability density of finding two electrons separated by the relative position vector , also known as the electron position intracule density , plays an important role in qualitative and quantitative descriptions of electronic systems. Position intracule densities have been extensively used to analyze shell structure, electron correlation, Hund’s rules and chemical bonding (see, e.g., Refs. CouNei-PPSL-61; LesKra-JCP-66; Koh-JCP-72; Kat-PRA-72; BoyCou-JPB-73; CooPou-TCA-77; ThaTriSmi-IJQC-84; RegTha-JPB-84; ShaTha-JPB-84; UgaBoy-IJQC-85; BoySarUga-JPB-88; SarUgaBoy-JPB-90; SarDomAguUga-JCP-92; WanTriSmi-JCP-92; CanBoyTha-JCP-93; WanSmi-IJQC-94; AriPorBueGal-JPB-95; MeyMulSch-JMS-96; CioLiu-JCP-96b; FraDurMes-JCP-97; MatKogRomDeh-PRA-98; SarGalBue-JCP-98; GalBueSar-JCP-99; CioLiu-JCP-99; GalBueSar-PRA-00; GalBueSar-JCP-02; GilOneBes-TCA-03; MerValUga-INC-03; GalBueSar-JCP-03; GalBueSar-JCP-05; GalBueSar-JCP-05a). Density functional theory-like approaches have been proposed based on the position intracule density Kog-JCP-90; GorSav-PRA-05; GorSav-PM-06; GorSav-JJJ-XX; Nag-JCP-06 or on the closely-related Wigner intracule density GilCriOneBes-PCCP-06; GilOne-JCP-05; Bes-JCP-06.

Position intracule densities have been extracted from experimental X-ray scattering intensities for small atoms and molecules ThaTriSmi-PRA-84; WanTriSmi-JCP-94; WatKamYamUdaMul-MP-04. They have been calculated in the Hartree-Fock (HF) approximation for systems ranging from small atoms to large molecules ThaTriSmi-IJQC-84; SarUgaBoy-JPB-90; WanSmi-IJQC-94; CioLiu-JCP-96a; CioLiu-JCP-96b; MatKogRomDeh-PRA-98; KogMatDehTha-JCP-99; LeeGil-CPL-99; GilLeeNarAda-JMS-00; GilOneBes-TCA-03. Calculations using common quantum chemistry correlated methods such as second-order Møller-Plesset perturbation theory, multi-configurational self-consistent-field (MCSCF) and configuration interaction approaches have been limited to atoms and small molecules CioLiu-JCP-98; CioLiu-JCP-99; Kog-JCP-02; CooPou-TCA-77; BoySarUga-JPB-88; WanTriSmi-JCP-92; MeyMulSch-JMS-96; MerValUga-INC-03. Very accurate calculations using Hylleraas-type explicitly-correlated wave functions have been done only for the helium and lithium isoelectronic series CouNei-PPSL-61; BoyCou-JPB-73; ThaSmi-JCP-77; RegTha-JPB-84; Tha-INC-87; CanBoyTha-JCP-93; DreKin-JCP-94; AriPorBueGal-JPB-95; CioSteTanUmr-JCP-95; Tha-CPL-03. Variational Monte Carlo (VMC) calculations using Jastrow-Slater wave functions have been used to compute correlated position intracule densities for atoms from helium to neon and some of their isoelectronic series CioSteTanUmr-JCP-95; SarGalBue-JCP-98; GalBueSar-JCP-99; GalBueSar-PRA-00; GalBueSar-JCP-02; GalBueSar-CPL-03; GalBueSar-CPL-03a; GalBueSar-JCP-03; GalBueSar-JCP-05; GalBueSar-JCP-05a. In this paper, we show that the calculation of position intracule densities using quantum Monte Carlo (QMC) methods can be made much more accurate and efficient, opening new possibilities of investigation.

The position intracule density associated with an -electron (real) wave function , where is the dimensional vector of electron coordinates (ignoring spin for now), is defined as the quantum-mechanical average of the delta-function operator

 I(u)=12∑i≠j∫dRΨ(R)2δ(rij−u), (1)

where and the sum is over all electron pairs, and its spherical average is

 I(u)=∫dΩu4πI(u). (2)

Among the most important properties of are the normalization sum rule (giving the total number of electron pairs)

 ∫∞0du4πu2I(u)=N(N−1)2, (3)

the electron-electron cusp condition ThaSmi-CPL-76

 dI(u)du∣∣∣u=0=I(0), (4)

and, for finite systems, the exponential decay at large determined by the spherically averaged one-electron density evaluated at the distance from the chosen origin

 I(u)∼u→∞(N−1)2n(u)∝u→∞e−2√2Iu, (5)

where is the vertical ionization energy (see Refs. LevPerSah-PRA-84; AlmBar-PRB-85; ErnBurPer-JCP-96). The moments of the radial intracule density are related to physical observables KogMat-JCP-01, in particular is just the electron-electron Coulomb interaction energy

 Wee=M−1=∫∞0du4πuI(u). (6)

In standard correlated methods based on an expansion of the wave function in Slater determinants, the important short-range part of the position intracule density converges very slowly with the one-electron and many-electron basis. The advantage of employing QMC methods FouMitNeeRaj-RMP-01 lies in the possibility of using compact, explicitly-correlated wave functions which are able to describe properly the short-range part of . However, the problem in this approach is that one calculates the average of a delta-function operator which has an infinite variance. As for other probability densities, the standard procedure in QMC approaches simply consists of counting the number of electron pairs separated by the distance within encountered in the Monte Carlo run. More precisely, e.g. in VMC calculations, is estimated as the statistical average

 Ihisto(u)=⟨IhistoL(u,R)⟩Ψ2, (7)

over the trial wave function density , of the local “histogram” estimator

 IhistoL(u,R)=12∑i≠j1[u−Δu/2,u+Δu/2[(rij)4πu2Δu, (8)

where if and otherwise. Here and in the following, designates the statistical average of the local quantity over configurations sampled from . For , this histogram estimator has a finite but large variance for small or small . Although the use of importance sampling can help decrease the variance LanRotVrb-JCP-97; SarGalBue-JCP-98; SarGalBue-CPC-99, the calculation of position intracule densities in QMC remains very inefficient: very long runs have to be performed to reach an acceptably small statistical uncertainty. Moreover, the histogram estimator has a discretization error: even in the limit of an infinite sample , remains only an approximation of first order in to the position intracule density of the wave function . Note however that it is possible to greatly reduce the discretization error by choosing a flexible analytic form for that obeys known conditions (such as the cusp condition of Eq. (4)) and fitting in each interval the integral of (instead of ) to the computed data. This approach has been used to calculate radial electron densities for atoms FilGonUmr-INC-96 and spherically-averaged intracule densities CioSteTanUmr-JCP-95 but the method is not as effective for multidimensional densities.

In this work, we develop improved QMC estimators for the spherically averaged position intracule density based on the zero-variance zero-bias (ZVZB) principle for observables introduced by Assaraf and Caffarel for calculations of electronic forces AssCaf-PRL-99; AssCaf-JCP-00; AssCaf-JCP-03, which has been recently applied to calculations of one-electron densities AssCafSce-PRE-07. These new ZVZB estimators of have variances several orders of magnitude smaller than the ones obtained with the standard estimator of Eq. (8) and thus dramatically increase the efficiency of these calculations. Moreover, these estimators do not suffer from any discretization error and can even reduce the systematic error due to the approximate trial wave function. Like related techniques proposed for calculations of one-electron or two-electron densities using deterministic ab initio methods HilSucFei-PRA-78; SucDra-PRA-79; Kat-PRA-80; Tri-JPB-80; Har-IJQC-80; Dra-JPB-81; Har-IJQC-85; Har-IJQC-86; MomShi-JCP-87; Ish-CPL-89; ChaCio-JCP-94; CioSteTanUmr-JCP-95; LiuParNag-PRA-95; RasChi-JCP-96a; RasChi-JCP-96b; WanSchSmi-PRA-00, these estimators replace the average of the local delta-function operator by the average of an operator which is nonlocal in real space. They can be viewed as a generalization of the improved QMC estimators previously proposed for computations of averages of probability densities at particle coalescences VrbDepRot-JCP-88; LanRotVrb-JCP-97; BreMelMor-PRA-98; AleCol-JMS-99.

Moreover, we make use in this work of the recently-developed linear energy minimization method TouUmr-JCP-07; UmrTouFilSorHen-PRL-07 to finely optimize the Jastrow parameters, the configuration state functions (CSFs) coefficients and the orbital coefficients of our trial Jastrow-Slater wave functions. Optimization of the determinantal part of the wave function with an increasing number of CSFs allows us to obtain a systematic improvement of wave functions. This provides a practical route for calculating intracule densities in VMC and fixed-node (FN) diffusion Monte Carlo (DMC) with progressively smaller systematic errors.

The paper is organized as follows. In Sec. II, we review the principle of zero-variance zero-bias improved estimators for an arbitrary observable in QMC. In Sec. III, we give improved estimators for the case of the position intracule density. Sec. IV contains computational details of the calculations, and Sec. V discusses results for the He and C atoms and for the C and N molecules to illustrate the technique. Sec. VI contains our conclusions.

Hartree atomic units are used throughout this work.

## Ii Zero-variance zero-bias improved estimators

We review here the principle of the zero-variance zero-bias improved QMC estimators for an arbitrary observable that does not commute with the Hamiltonian, developed in Refs. AssCaf-PRL-99; AssCaf-JCP-00; AssCaf-JCP-03. Throughout this section, all the averages are considered in the limit of an infinite Monte Carlo (MC) sample .

### ii.1 Estimators in variational Monte Carlo

In VMC, the exact energy of some exact eigenfunction of the electronic Hamiltonian is estimated by the statistical average of the local energy

 E=⟨EL(R)⟩Ψ2, (9)

using an approximate trial wave function . The systematic error (or bias) of this estimator and its variance , whose square root is proportional to the statistical uncertainty, both vanish quadratically as a function of the error in the trial wave function (where designates the Hilbert space norm)

 δE=O(|δΨ|2), (10a) σ2(EL)=O(|δΨ|2), (10b)

which is referred to as the quadratic zero-variance zero-bias property of the local energy. This is easily shown by writing in the expressions of the average and the variance, and expanding to second-order in .

Similarly, the exact expectation value of an arbitrary observable can be estimated by the statistical average of the local observable

 O=⟨OL(R)⟩Ψ2, (11)

but, as is generally not an eigenstate of , the systematic error of this estimator vanishes only linearly with , while its variance generally does not even vanish with

 δO=O(|δΨ|), (12a) σ2(OL)=O(1), (12b)

which often makes the calculations of observables inaccurate and inefficient.

However, Assaraf and Caffarel AssCaf-JCP-03 have pointed out that the quadratic zero-variance zero-bias property of the energy can be extended to an arbitrary observable by expressing it as an energy derivative. This is based on the Hellmann-Feynman theorem which states that, if one defines the -dependent Hamiltonian

 ^Hλ=^H+λ^O, (13)

with some associated exact -dependent eigenfunction

 Ψλ0=Ψ0+λΨ′0+⋯, (14)

where and corresponding exact -dependent energy , then the exact value of the observable is given by the derivative of the energy with respect to at : . Introducing an approximate -dependent trial wave function

 Ψλ=Ψ+λΨ′+⋯, (15)

where , the -dependent energy can be estimated as the statistical average of the -dependent local energy

 Eλ=⟨EλL(R)⟩Ψλ2, (16)

over the probability density , and the Hellmann-Feynman theorem suggests now to define a zero-variance zero-bias (ZVZB) estimator for the observable as the derivative of with respect to at

 OZVZB = ⟨OZVZBL(R)⟩Ψ2 (17) = (dEλdλ)λ=0 = +⟨ΔOZBL(R)⟩Ψ2,

where the zero-variance (ZV) contribution

 ΔOZVL(R) = [⟨R|^H|Ψ′⟩Ψ′(R)−EL(R)]Ψ′(R)Ψ(R), (18)

does not contribute to the average, (in the limit of an infinite MC sample , from the Hermiticity of the Hamiltonian) but can decrease the variance, and the zero-bias (ZB) contribution

 ΔOZBL(R)=2[EL(R)−E]Ψ′(R)Ψ(R), (19)

usually has little effect on the variance but can decrease the systematic error. The average of this ZB term vanishes if the trial wave function is an exact eigenstate or, more generally, if the energy is stationary with respect to an infinitesimal variation of the trial wave function along the derivative : . For the special case of the calculation of the derivatives of the energy with respect to the nuclear coordinates, the average of the ZB term is known in the electronic structure literature as the Pulay contribution Pul-MP-69.

If one wants simply to compute an observable for a given trial wave function without changing the average, the ZB correction has to be omitted. In this case, the ZV estimator has a systematic error that still vanishes linearly with but a variance that now vanishes quadratically with and the error in its derivative

 δOZV=O(|δΨ|), (20a) σ2(OZVL)=O(|δΨ|2+|δΨ′|2+|δΨ||δΨ′|). (20b)

If one also includes the ZB correction, then both the systematic error and the variance of the ZVZB estimator vanish quadratically with and

 δOZVZB=O(|δΨ|2+|δΨ||δΨ′|), (21a) σ2(OZVZBL)=O(|δΨ|2+|δΨ′|2+|δΨ||δΨ′|). (21b)

This is easily shown by writing and in the expressions of the average and the variance, and expanding with respect to and . One can verify that the exact cancellation of the dominant contributions to the systematic error and the variance of Eqs. (12a) and (12b) by the ZV and ZB corrections is a direct consequence of the relationship between the exact wave function and its derivative in first-order perturbation theory: .

Thus, the ZVZB estimators permit in principle accurate and efficient calculations for observables, provided that a good approximation for the derivative is available. We note that optimizing the parameters in the trial wave function by energy minimization presents the advantage of decreasing the sensitivity of the systematic error to the quality of the approximate derivative , since in this case only needs to be a good approximation of in the orthogonal complement of the space spanned by the wave function derivatives with respect to the parameters at the optimal parameter values , i.e. . Indeed, the contribution to the average of the ZB term of any component of along is proportional to the energy gradient with respect to parameter and thus vanishes at the stationary point: . This point is at the origin of the improved accuracy in the calculations of forces in QMC in Refs. CasMelRap-JCP-03; LeeMelRap-JCP-05.

### ii.2 Estimators in diffusion Monte Carlo

DMC calculations within the FN approximation go beyond VMC calculations by replacing the VMC distribution by the more accurate mixed FN-DMC distribution in statistical averages, i.e. . All the DMC variances have the same lowest order behaviors with respect to the error in the trial wave function as in VMC. The systematic error of the energy now vanishes as the product of the error in the trial wave function and the error in the FN wave function

 δE=O(|δΨ||δΨFN|). (22)

For an arbitrary observable, the systematic error still vanishes only linearly with or

 δO=O(|δΨ|+|δΨFN|), (23)

but the use of the hybrid estimator “”, , allows one to remove the dominant linear term

 δOhybrid=O(|δΨFN|+|δΨ|2+|δΨFN|2+|δΨ||δΨ% FN|).

The same improved estimators defined in the previous section can also be used straightforwardly in FN-DMC calculations. Note that, in this case, the average of the ZV term of Eq. (18) no longer vanishes on an infinite MC sample. Nevertheless, this term gives the same order of reduction of the variance as in VMC. The systematic error and the variance are still given by Eqs. (20a) and (20b), though the prefactors can be different in VMC and DMC. We note that it is possible to define a new ZV correction of vanishing average in FN-DMC AssCaf-JCP-00 but this does not change the leading order of either the systematic error or the variance.

Adding the ZB term of Eq. (19), the systematic error of the ZVZB estimator vanishes quadratically

 δOZVZB = O(|δΨ|2+|δΨ||δΨ′|+|δΨ||δΨFN| (25) +|δΨ′||δΨFN|),

with no linear term in contrast to the hybrid estimator of Eq. (LABEL:deltaOhyb). The behavior of the variance of the ZVZB estimator in FN-DMC is still given by Eq. (21b). One can also define a hybrid ZVZB estimator, , whose systematic error also vanishes quadratically as

 δOZVZBhybrid = O(|δΨ|2+|δΨ||δΨFN|+|δΨ′||δΨFN|),

with no term in in contrast to Eq. (25).

### ii.3 Expressions of the estimators for local Hamiltonians

We give now more explicit expressions of the estimators in terms of the convenient logarithmic derivative .

In the case of local Hamiltonians where contains kinetic and potential energy contributions, the ZV term takes the form

 ΔOZVL(R) =

where is the drift velocity of the trial wave function. The ZB term is simply

 ΔOZBL(R)=2ΔEL(R)Q(R), (28)

where . For each observable, a specific form for has to be chosen. In principle, one could optimize for each system by minimizing the variance of the ZVZB estimator. Besides , the ZV and ZB terms require only the evaluation of the drift velocity and the local energy of the trial wave function that are already available in QMC programs.

Clearly, a nonlocal operator in the Hamiltonian, e.g. a nonlocal pseudopotential, would yield an additional term in the simplified expression of the ZV correction of Eq. (LABEL:DOLZVHloc). Although the inclusion of this additional term is required for the variance to rigorously vanish quadratically as in Eq. (20b), in practice the estimator using the simple form of Eq. (LABEL:DOLZVHloc) can already achieve a considerable variance reduction.

## Iii Zero-variance zero-bias improved estimators for the spherically averaged position intracule density

In the case of the spherically averaged position intracule density, the local observable is

 IL(u,R)=12∑i≠j∫dΩu4πδ(rij−u), (29)

which has an infinite variance. We will give two possible choices for which define good improved estimators.

### iii.1 First improved estimator

The minimal choice for that cancels the infinite variance of is

 Q1(u,R)=−18π∑i≠j∫dΩu4π1|rij−u|, (30)

which gives the following ZV contribution

 ΔIZV1L(u,R) = 18π∑i≠j∫dΩu4π[∇2rij1|rij−u| (31) +2vi(R)⋅∇rij1|rij−u|].

Using the identity , it is seen that the first term in exactly cancels the delta function in , and the ZV estimator is simply

 IZV1L(u,R) = 14π∑i≠jvi(R)⋅∫dΩu4π∇rij1|rij−u| (32) = −14π∑i≠jvi(R)⋅rijr3ij1[u,+∞[(rij).

The form of of Eq. (30) could also have been guessed from the behavior at small of the first-order solution of the perturbation theory with respect to (see Refs. Sch-AP-59; Tri-JPB-80). In contrast to the histogram estimator of Eq. (8) where only electron pairs with relative distance in the small interval contribute to , for the ZV estimator of Eq. (32) all the electron pairs with relative distance contribute to . For , the variance of this ZV estimator is finite and much smaller than the variance of the histogram estimator as shown in Sec. V. For , this ZV estimator reduces to the estimator proposed in Ref. LanRotVrb-JCP-97 to compute one-electron densities at the nucleus and applied in Ref. SarGalBue-JCP-98 for calculations of intracule densities at the coalescence point.

The ZB correction corresponding to of Eq. (30) is

 ΔIZB1L(u,R)=−ΔEL(R)4π∑i≠j∫dΩu4π1|rij−u| =−ΔEL(R)4π∑i≠j[1[0,u[(rij)u+1[u,+∞[(rij)rij]. (33)

The resulting ZVZB estimator allows one to reduce the systematic error due to the trial wave function as shown in Sec. V. Note that, in contrast to the histogram estimator, these ZV and ZVZB improved estimators do not suffer from any discretization error: they are estimators of for an infinitely precise value of . Both the ZV estimator of Eq. (32) and its ZB correction of Eq. (33) are very simple to program and very fast to compute on a grid over .

Like the histogram estimator, the ZV estimator of Eq. (32) gives simply zero for distances beyond the largest electron-electron distance encountered in the Monte Carlo run. Also, because the form of of Eq. (30) has been chosen only to cure the deficiency of the histogram estimator at small , the ZB correction of Eq. (33) actually makes the histogram estimator worse at large . The ZB correction decays only as as instead of the correct exponential decay and consequently gives large variances at large . Thus, it is inaccurate to extract the long-range behavior of the intracule density using the ZVZB estimator . The second choice for presented next remedies this limitation.

### iii.2 Second improved estimator

A more general choice for that cancels the infinite variance of is

 Q2(u,R)=−18π∑i≠j∫dΩu4πg(rij,u)|rij−u|, (34)

where is some function satisfying for . For example, using

 g(rij,u)=e−ζ|rij−u|, (35)

with , has the advantage of ensuring a correct exponential decay of the estimator at infinity for a finite system.

The corresponding ZV estimator is

 IZV2L(u,R) = 18π∑i≠j∫dΩu4π[2vi(R)⋅∇rije−ζ|rij−u||rij−u| (36) +2∇rij1|rij−u|⋅∇rije−ζ|rij−u| +1|rij−u|∇2rije−ζ|rij−u|] = −14π∑i≠j[vi(R)⋅rijr3ij(rijA(rij,u) +rijuB(rij,u))−ζ22A(rij,u)],

and its ZB correction is

 ΔIZB2L(u,R) = −ΔEL(R)4π∑i≠j∫dΩu4πe−ζ|rij−u||rij−u| (37) = −ΔEL(R)4π∑i≠jA(rij,u),

where

 A(rij,u)=e−ζ|rij−u|−e−ζ(rij+u)2ζriju, (38)
 B(rij,u)=12[sgn(rij−u)e−ζ|rij−u|−e−ζ(rij+u)], (39)

and where sgn is the sign function. These refined ZV estimator and ZVZB estimator now both decay correctly exponentially as . The constant is chosen according to Eq. (5), i.e. where is an estimate of the vertical ionization energy. On the other hand, the computational cost is larger per Monte Carlo step than that of the simpler estimators of the previous subsection.

## Iv Computational details

We have chosen to illustrate the efficiency and accuracy of the new improved QMC estimators by calculating the intracule density of the He and C atoms and the C and N molecules in their ground states. The molecules are considered at their experimental geometries ( Bohr CadWah-ADNDT-74 and Bohr HubHer-BOOK-79).

We start by generating a standard ab initio wave function using the quantum chemistry program GAMESS SchBalBoaElbGorJenKosMatNguSuWinDupMon-JCC-93, typically a HF wave function or a MCSCF wave function with a complete active space (CAS) generated by distributing valence electrons in valence orbitals [CAS(,)]. For each system considered, we choose the one-electron Slater basis so as to ensure that the HF intracule density is reasonably converged with respect to the basis. For the He atom, we use the basis of Clementi and Roetti CleRoe-ADNDT-74, with exponents reoptimized at the HF level by Koga et al. KogTatTha-PRA-93. For the C atom, we use the CVB1 basis of Ema et al. EmaGarRamLopFerMeiPal-JCC-03. For the C and N molecules, we use the CVB2 basis of the same authors. In GAMESS, each Slater function is actually approximated by a fit to 14 Gaussian functions HehStePop-JCP-69; Ste-JCP-70; KolReiAss-JJJ-XX

This standard ab initio wave function is then multiplied by a Jastrow factor, imposing the electron-electron cusp condition, but with essentially all other free parameters chosen to be zero to form our starting trial Jastrow-Slater wave function, and QMC calculations are performed with the program CHAMP Cha-PROG-XX or QMCMOL Qmc-PROG-XX using the true Slater basis set rather than its Gaussian expansion. The Jastrow parameters, the orbital coefficients and the configuration state functions (CSFs) coefficients of this trial wave function are optimized in VMC using the very efficient, recently-developed linear energy minimization method TouUmr-JCP-07; UmrTouFilSorHen-PRL-07 and an accelerated Metropolis algorithm Umr-PRL-93; Umr-INC-99. Once the trial wave function has been optimized, we compute the intracule density in VMC, and in DMC using the fixed-node and the short-time approximations (see, e.g., Refs. And-JCP-75; And-JCP-76; ReyCepAldLes-JCP-82; MosSchLeeKal-JCP-82). We use an imaginary time-step of hartree in an efficient DMC algorithm featuring very small time-step errors UmrNigRun-JCP-93.

## V Results and discussion

### v.1 He atom

The improvement due to the ZV and ZVZB estimators is illustrated for the simple case of He atom.

Figure 1 shows the spherically averaged intracule density as a function of the electron-electron distance , calculated in VMC using the standard histogram estimator of Eq. (8), the ZV1 improved estimator of Eq. (32) and the ZV1ZB1 improved estimator of Eq. (33), employing only 100 000 MC configurations sampled from a HF trial wave function (without a Jastrow factor). In this subsection, we intentionally use a poor trial wave function to demonstrate the large reduction in the bias resulting from the ZVZB estimator. The calculations take a few seconds on a present-day single-processor personal computer. The intracule density obtained with an accurate, explicitly-correlated Hylleraas-type wave function with 491 terms FreHuxMor-PRA-84; UmrGon-PRA-94; CioSteTanUmr-JCP-95 is also shown as a reference. The statistical uncertainty obtained with the histogram estimator is very large, especially at small , making it impossible to extract the on-top intracule density . The ZV1 estimator spectacularly reduces the statistical uncertainty which becomes invisible at the scale of the plot, its variance being smaller by nearly 4 orders of magnitude for small and large and by about orders of magnitude for intermediate . The ZV1ZB1 estimator in turn spectacularly reduces the systematic error due to the use of a HF trial wave function, the obtained intracule density agreeing well with the accurate reference. In particular, the maximum at small caused by the cusp at is accurately reproduced. On the other hand, the use of the ZV1ZB1 estimator increases slightly the variance at small and more importantly at large where the ZB1 correction of Eq. (33) decays too slowly as .

The behavior of the estimators in the long-range tail of is explored in detail in Fig. 2 which compares the ZV1ZB1 and ZV2ZB2 estimators of Eqs. (36) and (37) for . As mentioned in Sec. III.2, the relative statistical uncertainty obtained with the ZV1ZB1 is very large (about at ), whereas the ZV2ZB2 estimator, which has the correct exponential decay at large , spectacularly reduces the statistical uncertainty, its variance being smaller by about 3 orders of magnitude at . The improvement is even more dramatic for larger .

We now discuss the scaling of the intracule density and its statistical uncertainty with respect to the nuclear charge . For this purpose, we have calculated the intracule densities of the elements of the He isoelectronic series from (He) to (Ar) using Jastrow-Slater wave functions. We found that, in the relevant range of (), the intracule density roughly obeys the scaling law with , in agreement with the value predicted by a simple hydrogenic model KogMat-JCP-97. The statistical uncertainty of computed with the ZV1 estimator roughly obeys the same law, so that the relative statistical accuracy on does not deteriorate with in the series (at least up to ).

Finally, we note in passing that Fig. 2 also reveals another advantage of our improved estimators: because the estimator of is strongly correlated with the estimator of , as obvious for instance in Eq. (32), the obtained intracule density curves are always very smooth, even at a scale much smaller than the statistical uncertainty. This makes the calculated intracule densities directly suitable for subsequent manipulations.

### v.2 C atom

We continue the discussion of the improved estimators for the more interesting case of the C atom.

The repercussion of the shell structure of the system on the intracule is apparent in Fig. 3 which shows the integrand of the electron-electron Coulomb interaction energy over the electron-electron distance , i.e. [see Eq. (6)]. The ZVZB2 estimator is used with an accurate Jastrow-Slater CAS(4,4) wave function in which the Jastrow, CSF and orbital parameters have been simultaneously optimized. Note however that, at the scale of the plot, the HF intracule would look the same. The curve displays two maxima: the maximum at short distance () corresponds essentially to the K-K electron pair and the maximum at longer distance () corresponds essentially to the K-L and L-L electron pairs. Fig. 3 also shows the same spin (SS) contribution

 ISS(u)=I↑↑(u)+I↓↓(u), (40)

and opposite spin (OS) contribution

 IOS(u)=I↑↓(u)+I↓↑(u), (41)

where is the spherical average of the spin-resolved intracule densities

 Iσσ′(u)=12∑i≠j∫dRΨ(R)2δ(rij−u)δsi,σδsj,σ′, (42)

where the spin coordinates of the electrons in the spin-assigned wave function are fixed at for and for . Not surprisingly, the intracule density at short distance is dominated by the OS contribution, while at long distance the SS and OS contributions become nearly identical.

Figure 4 shows the correlation part of the radial intracule density , also called correlation hole, using the ZVZB2 improved estimator and the same Jastrow-Slater CAS(4,4) trial wave function. Correlation decreases the radial probability density at the two previously-mentioned maxima and increases it at longer distances. The same spin contribution and the opposite spin contribution are also shown in Fig. 4. The OS component constitutes the main contribution to the correlation hole.

The effect of the chosen trial wave function on the accuracy on the correlation hole is examined in Fig. 5. Four trial wave functions of increasing accuracy have been tested: HF, Jastrow HF (with only Jastrow parameters optimized), Jastrow single-determinant (SD) (with Jastrow and orbital parameters optimized) and Jastrow CAS(4,4) (with Jastrow, CSF and orbital parameters optimized). Remarkably, the ZVZB2 estimator gives a correlation hole with a correct overall structure even with the uncorrelated HF wave function. For more quantitative predictions, the use of a Jastrow factor is however necessary. Optimization of the orbitals in a Jastrow-Slater single-determinant wave function brings a further significant improvement in the interesting short-range region (). The Jastrow-Slater multi-determinant CAS(4,4) and the Jastrow-Slater single-determinant intracules agree closely.

### v.3 C2 molecule

We now discuss the more difficult case of the C molecule. The ground-state wave function of this system has a strong multi-configurational character due to energetic near-degeneracies among the valence orbitals (“nondynamical correlation” in chemists’ jargon, or “strong correlation” in physicists’ jargon), making it a challenging system despite its small size.

Figure 6 plots versus using the ZVZB2 estimator with an accurate Jastrow-Slater CAS(8,8) wave function (Jastrow, CSF and orbital parameters optimized). The curve displays three maxima: the maximum at short distance () corresponds essentially to the intra-atomic K-K electron pairs; the maximum at long distance () is located at about the interatomic distance and corresponds mainly to the interatomic K-K electron pairs; the maximum at intermediate distance () must correspond mainly to intra-atomic and interatomic K-L and L-L electron pairs. As for the C atom, the intracule density at short distance is dominated by the OS contribution, while at long distance the SS and OS contributions become nearly identical.

Figure 7 shows the correlation hole of the C molecule using the ZV2 improved estimator calculated in VMC with a MCSCF CAS(8,8) wave function and a series of Jastrow CAS(8,8) wave functions for different levels of optimization. The MCSCF CAS(8,8) wave function does not correlate the core electrons and thus gives essentially no correlation hole at distances corresponding to intra-atomic core electron pairs (). At longer valence distances, the MCSCF CAS(8,8) wave function gives only a gross overall shape of the correlation hole. Introduction of a Jastrow factor yields a correct correlation hole at core distances, which as expected is about twice as deep as the core correlation hole of the C atom, and reduces the correlation hole at valence distances. The action of the Jastrow factor is thus not limited to very short electron-electron distances, but in fact importantly modifies the correlation hole up to distances 4. The Jastrow factor has also an indirect action via reoptimization of the CSF and orbital coefficients in its presence. The reoptimization of the CSF coefficients changes slightly the correlation hole at valence distances. The reoptimization of the CSF and orbital coefficients has a more important impact of the correlation hole at valence distances and also at core distances.

We next scrutinize in detail the convergence of the intracule density with respect to the determinantal part of the trial wave function. Figure 8a shows the correlation hole using the ZVZB2 improved estimator calculated in VMC with five trial Jastrow-Slater wave functions of increasing accuracy: Jastrow HF (with only the Jastrow parameters optimized), Jastrow SD (with the Jastrow and orbital parameters optimized), and multi-determinant Jastrow CAS(8,5), Jastrow CAS(8,7) and Jastrow CAS(8,8) (with the Jastrow, CSF and orbital parameters optimized). At core distances (), the correlation hole is essentially converged with only a single-determinant Jastrow-Slater wave function, provided that the orbitals are reoptimized together with the Jastrow factor. In contrast, at longer valence distances, the correlation hole depends very strongly on the determinantal part of the wave function. In particular, the multi-determinant wave functions yield a depletion of the intracule density at about the bond length, a feature not present at the single-determinant level. We have verified that this minimum disappears when using a pseudopotential removing the electrons, and we thus interpret it as a decrease of probability of finding two interatomic core electrons separated by the bond distance. Overall, it appears necessary to use at least a multi-determinant CAS(8,7) wave function which includes configurations constructed from the antibonding orbitals to reach reasonable convergence. We have checked that further excitations beyond the CAS(8,8) wave function barely change the correlation hole. The corresponding correlation holes calculated in FN-DMC with the same ZVZB2 estimator and trial wave functions are reported in Fig. 8b. The DMC intracules do not differ much from the corresponding VMC intracules, the accuracy of the intracule densities in the valence region being still essentially controlled by the trial wave function.

### v.4 N2 molecule

We finish our illustration of the calculation of intracule densities with the N molecule, which has recently been the object of experimental and theoretical investigations WatKamYamUdaMul-MP-04.

The correlation holes of this molecule calculated in VMC with the ZVZB2 improved estimator for a Jastrow HF, and fully optimized Jastrow SD and Jastrow CAS(10,8) wave functions are plotted in Fig. 9. In sharp contrast to the C molecule, no depletion of intracule density is observed at the bond length. Also, the correlation hole is here essentially converged with only a single-determinant Jastrow-Slater wave function provided that the orbitals are optimized with the Jastrow factor. Unlike the C molecule, going to a multi-determinant wave function does not have a large effect on the correlation hole at valence distances. The correlation hole calculated here agrees reasonably well with recent multi-reference configuration interaction and coupled-cluster calculations, and with experiment WatKamYamUdaMul-MP-04.

## Vi Conclusions

We have presented improved QMC estimators for the spherically averaged position intracule density , constructed using the general zero-variance zero-bias principle for observables that do not commute with the Hamiltonian. By replacing the average of the local delta-function operator by the average of a smooth non-local operator, these estimators decrease the variance of the standard histogram estimator by several orders of magnitude, and thus make the calculation of this quantity in QMC vastly more efficient. Interestingly, they permit calculations of for very short and very large interelectronic distances that are never realized in the Monte Carlo run. These new estimators can also decrease the systematic error of the intracule density due to the approximate trial wave function. Other advantages of these estimators are the absence of any discretization error with respect to and the possibility to obtain very smooth curves for . These improved estimators, together with the achievement of systematically reducing the systematic error in both VMC and DMC calculations by optimization of trial wave functions with an increasing number of parameters, have allowed us to obtain accurate correlated intracule densities for atoms and molecules.

The estimators presented here can be used with trivial adaptations for QMC calculations of the companion entity of , namely the spherically averaged extracule density , representing the probability density of finding two electrons with center of mass at a radial distance with respect to the chosen origin. Similar improved estimators can be constructed for the three-dimensional intracule density , extracule density and the full pair density . More generally, the variance reduction technique presented here can be applied to calculations of any pair-correlation function in classical and quantum Monte Carlo calculations.

###### Acknowledgements.
We would like to thank Andreas Savin and Paola Gori-Giorgi for stimulating discussions. We also thank Alexander Kollias for providing us the Gaussian fits of Slater basis functions of Ref. KolReiAss-JJJ-XX. J. T. acknowledges financial support from a Marie Curie Outgoing International Fellowship (039750-QMC-DFT). This work was also supported by the Centre National de la Recherche Scientifique and by the National Science Foundation (DMR-0205328, EAR-0530301). Most of the calculations were been performed at the Cornell Theory Center and on the Intel cluster at the Cornell Nanoscale Facility (a member of the National Nanotechnology Infrastructure Network supported by the National Science Foundation).

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters