Self-Consistent Electron-Nucleus Cusp Correction for Molecular Orbitals

TL;DR

This paper introduces a method to impose the correct electron-nucleus cusp in molecular orbitals by augmenting Gaussian basis sets with Slater functions, significantly reducing energy variance in quantum Monte Carlo calculations.

Contribution

The method uniquely combines Gaussian and Slater basis functions with self-consistent optimization to improve electron-nucleus cusp representation in molecular orbitals.

Findings

Variance reduced by an order of magnitude for helium atom

Method effective for atoms and small molecules

Self-consistent optimization enhances cusp correction

Abstract

We describe a method for imposing the correct electron-nucleus (e-n) cusp in molecular orbitals expanded as a linear combination of (cuspless) Gaussian basis functions. Enforcing the e-n cusp in trial wave functions is an important asset in quantum Monte Carlo calculations as it significantly reduces the variance of the local energy during the Monte Carlo sampling. In the method presented here, the Gaussian basis set is augmented with a small number of Slater basis functions. Note that, unlike other e-n cusp correction schemes, the presence of the Slater function is not limited to the vicinity of the nuclei. Both the coefficients of these cuspless Gaussian and cusp-correcting Slater basis functions may be self-consistently optimized by diagonalization of an orbital-dependent effective Fock operator. Illustrative examples are reported for atoms (\ce{H}, \ce{He} and \ce{Ne}) as well as for a small molecular system (\ce{BeH2}). For the simple case of the \ce{He} atom, we observe that, with respect to the cuspless version, the variance is reduced by one order of magnitude by applying our cusp-corrected scheme.

Figures

Tables

Table 1. Table 1: Variational energy and its variance for various wave functions of the \ce H atom. The energy is obtained with the decontraced STO-3G Gaussian basis set. The OS and SCD cusp-corrected energies are obtained by adding a Slater basis function of unit exponent ( α ~ H = 1 subscript ~ 𝛼 H 1 \tilde{\alpha}_{\text{H}}=1 ) to the Gaussian basis set. The energy and variance at each iteration of the SCD process is also reported. max ⁡ | 𝒄 | 𝒄 \max{\absolutevalue{\bm{c}}} is the maximum absolute value of the Gaussian basis coefficients and c ~ ~ 𝑐 \tilde{c} is the value of the coefficient of the Slater function.

Basis	Cusp correction	Iteration	Energy	Variance	$\max \left|𝒄\right|$	$\tilde{c}$
Gaussian	—		$-0.495741$	$2.23\times {10}^{-1}$	0.66158	0
Mixed	OS		$-0.499270$	$4.49\times {10}^{-2}$	0.07486	1.95629
Mixed	SCD	#1	$-0.499270$	$4.49\times {10}^{-2}$	0.07486	1.95629
		#2	$-0.499970$	$3.07\times {10}^{-6}$	0.00254	2.00225
		#3	$-0.500000$	$4.88\times {10}^{-9}$	0.00006	1.99691

Table 2. Table 2: Energy and corresponding variance of the \ce He atom computed with various methods. The trial wave function is the HF wave function computed in the 6-31G basis. The error bar corresponding to one standard deviation is reported in parenthesis.

Cusp	Energy (a.u.)			Variance (a.u.)
correction	Deterministic	VMC	DMC	VMC	DMC
—	$-\mathrm{2.855\hspace{0.17em}160}$	$-\mathrm{2.855\hspace{0.17em}12}\left(6\right)$	$-\mathrm{2.903\hspace{0.17em}9}\left(1\right)$	$3.99\left(3\right)$	$4.47\left(18\right)$
OS		$-\mathrm{2.857\hspace{0.17em}89}\left(6\right)$	$-\mathrm{2.903\hspace{0.17em}4}\left(3\right)$	$0.605\left(6\right)$	$0.498\left(2\right)$
SCD		$-\mathrm{2.858\hspace{0.17em}17}\left(9\right)$	$-\mathrm{2.903\hspace{0.17em}2}\left(2\right)$	$0.610\left(3\right)$	$0.498\left(1\right)$

Equations

$${\Psi }_{\text{T}}(\mathbf{R})=e^{J(\mathbf{R})}\sum _I{c}_I{D}_I^{\uparrow }({\mathbf{R}}^{\uparrow })D_I^{\downarrow }({\mathbf{R}}^{\downarrow }),$$

$${\varphi }_i(\mathbf{r})=\sum _{\mu }^N{c}_{\mu i}{\chi }_{\mu }(\mathbf{r})$$

$$\text{evaluated}{\text{partialderivative}\text{expectationvalue}\Psi (\mathbf{R})r_i}_{r_i=r_A}=-Z_A{\text{expectationvalue}\Psi (\mathbf{R})|}_{r_i=r_A},$$

$$\text{evaluated}{\text{partialderivative}\text{expectationvalue}{\varphi }_i(\mathbf{r})r}_{r=r_A}=0.$$

$$\evaluated{\partialderivative{\expectationvalue*{\tilde{\phi_{i}}(\bm{r})}}{r}}_{r=r_{A}}=-Z_{A}\expectationvalue*{\tilde{\phi_{i}}(\bm{r})}\big{|}_{r=r_{A}}.$$

$$\widetilde{{\varphi }_i}(\mathbf{r})={\varphi }_i(\mathbf{r})+\hat{P}{\phi }_i(\mathbf{r}),$$

$${\phi }_i(\mathbf{r})=\sum _A^M{\tilde{c}}_{Ai}{{\tilde{\chi }}_A}^i(\mathbf{r}),$$

$${{\tilde{\chi }}_A}^i(\mathbf{r})=\sqrt{\frac{{\tilde{\alpha }}_i^3}{\pi }}\exp [-{\tilde{\alpha }}_i\text{absolutevalue}\mathbf{r}-{\mathbf{r}}_A]$$

$$\hat{P}=\hat{I}-\sum _{\mu }\text{outerproduct}{\chi }_{\mu }{\chi }_{\mu }$$

$$\sum_{B}\Bigg{[}-\frac{\delta_{AB}}{Z_{A}}\partial_{r}{\tilde{\chi}_{A}}^{i}(\bm{r}_{A})-{\tilde{\chi}_{B}}^{i}(\bm{r}_{A})+\sum_{\mu}\tilde{S}_{B\mu}^{i}\chi_{\mu}(\bm{r}_{A})\Bigg{]}\tilde{c}_{Bi}=\phi_{i}(\bm{r}_{A}),$$

$${\delta }_{AB}=\{\begin{array}{ll}1,& A=B,\\ 0,& A\ne B\end{array}$$

$${\partial }_r{{\tilde{\chi }}_A}^i({\mathbf{r}}_A)\equiv \text{evaluated}{\text{partialderivative}{{\tilde{\chi }}_A}^i(\mathbf{r})r}_{r=r_A}.$$

$$\hat{f}|\ast ⟩\widetilde{{\varphi }_i}={\tilde{\epsilon }}_i|\ast ⟩\widetilde{{\varphi }_i}.$$

$$\sum _{\nu }{F}_{\mu \nu }{c}_{\nu i}+\sum _A{\tilde{c}}_{Ai}\text{quantity}({\tilde{F}}_{\mu A}^i-\sum _{\lambda }{F}_{\mu \lambda }{\tilde{S}}_{\lambda A}^i)={\tilde{\epsilon }}_i{c}_{\mu i},$$

$$F_{\mu \nu }$$

$${\tilde{F}}_{\mu A}^i-\sum _{\lambda }{F}_{\mu \lambda }{\tilde{S}}_{\lambda A}^i\approx {\tilde{h}}_{\mu A}^i-\sum _{\lambda }{H}_{\mu \lambda }^{\text{c}}{\tilde{S}}_{\lambda A}^i$$

$$H_{\mu \nu }^{\text{c}}$$

$$\sum _{\nu }{\tilde{F}}_{\mu \nu }^i{c}_{\nu i}={\tilde{\epsilon }}_i{c}_{\mu i},$$

$${\tilde{F}}_{\mu \nu }^i=\{\begin{array}{ll}{F}_{\mu \mu }+{\tilde{D}}_{\mu }^i,& \text{if μ=ν},\\ F_{\mu \nu },& \text{otherwise},\end{array}$$

$${\tilde{D}}_{\mu }^i=c_{\mu i}^{-1}\sum _A{\tilde{c}}_{Ai}\text{quantity}({\tilde{h}}_{\mu A}^i-\sum _{\lambda }{H}_{\mu \lambda }^{\text{c}}{\tilde{S}}_{\lambda A}^i).$$

$$P_{\mu \nu }=\sum _i^{\text{occ}}{c}_{\mu i}{c}_{\nu i}.$$

$${\tilde{\alpha }}_i=\frac{{\varphi }_i({\mathbf{r}}_A)}{\mathring{{\varphi }_i}({\mathbf{r}}_A)}{Z}_A$$

$${\tilde{S}}_{\mathbf{a}B}=\text{innerproduct}{\chi }_{\mathbf{a}}{\tilde{\chi }}_B=\int {\chi }_{\mathbf{a}}(\mathbf{r}){\tilde{\chi }}_B(\mathbf{r})d\mathbf{r},$$

$${\chi }_{\mathbf{a}}(\mathbf{r})=(x-A_x{)}^{a_x}(y-A_y{)}^{a_y}(z-A_z{)}^{a_z}\exp [-\alpha \text{absolutevalue}{\mathbf{r}-\mathbf{A}}^2]$$

$${\tilde{S}}_{(1,0,0)B}=\frac{1}{2\alpha }\text{partialderivative}{\tilde{S}}_{0B}{A}_x.$$

$$\exp (-\zeta r)=\frac{\zeta }{\sqrt{\pi }}{\int }_0^{\infty }\exp [-u^{-2}{r}^2-\frac{{\zeta }^2{u}^2}{4}]du,$$

$$\footnotesize\begin{split}\tilde{S}_{\mathbf{0}B}&=\frac{\tilde{\beta}^{5/2}}{\pi}\int_{0}^{\infty}\exp[-\frac{\tilde{\beta}^{2}u^{2}}{4}]\quantity(\int\exp[-\alpha\absolutevalue{\bm{r}-\mathbf{A}}^{2}]\exp[-u^{-2}\absolutevalue{\bm{r}-\mathbf{B}}^{2}]d\bm{r})du\\ &=\frac{\tilde{\beta}^{5/2}}{\pi}\int_{0}^{\infty}\exp[-\frac{\tilde{\beta}^{2}u^{2}}{4}]\quantity(\frac{\pi}{\alpha+u^{-2}})^{3/2}\exp[-\frac{\alpha\,u^{-2}}{\alpha+u^{-2}}\mathbf{A}\mathbf{B}^{2}]du\\ &=\frac{\pi\tilde{\beta}^{3/2}}{2\alpha^{2}\absolutevalue{\mathbf{A}\mathbf{B}}}\Bigg{\{}\quantity(\sqrt{\frac{\tilde{\beta}^{2}}{4\alpha}}+\sqrt{\alpha\absolutevalue{\mathbf{A}\mathbf{B}}})\operatorname{erfc}\quantity(\sqrt{\frac{\tilde{\beta}^{2}}{4\alpha}}+\sqrt{\alpha\absolutevalue{\mathbf{A}\mathbf{B}}})\exp[\frac{\tilde{\beta}^{2}}{4\alpha}+\tilde{\beta}\absolutevalue{\mathbf{A}\mathbf{B}}]\\ &\qquad\qquad\quad-\quantity(\sqrt{\frac{\tilde{\beta}^{2}}{4\alpha}}-\sqrt{\alpha\absolutevalue{\mathbf{A}\mathbf{B}}})\operatorname{erfc}\quantity(\sqrt{\frac{\tilde{\beta}^{2}}{4\alpha}}-\sqrt{\alpha\absolutevalue{\mathbf{A}\mathbf{B}}})\exp[\frac{\tilde{\beta}^{2}}{4\alpha}-\tilde{\beta}\absolutevalue{\mathbf{A}\mathbf{B}}]\Bigg{\}},\end{split}$$

$${\tilde{T}}_{\mathbf{a}B}=-\frac{1}{2}\text{innerproduct}{\nabla }^2{\chi }_{\mathbf{a}}{\tilde{\chi }}_B,$$

$$\footnotesize\begin{split}\tilde{T}_{\mathbf{0}B}&=\frac{\pi\tilde{\beta}^{5/2}}{2\alpha^{3/2}\absolutevalue{\mathbf{A}\mathbf{B}}}\Bigg{\{}\quantity(1+\frac{\tilde{\beta}^{2}}{4\alpha}-\sqrt{\frac{\tilde{\beta}^{2}}{4}}\absolutevalue{\mathbf{A}\mathbf{B}})\operatorname{erfc}\quantity(\sqrt{\frac{\tilde{\beta}^{2}}{4\alpha}}-\sqrt{\alpha}\absolutevalue{\mathbf{A}\mathbf{B}})\exp[\frac{\tilde{\beta}^{2}}{4\alpha}-\tilde{\beta}\absolutevalue{\mathbf{A}\mathbf{B}}]\\ &\qquad\qquad\quad-\quantity(1+\frac{\tilde{\beta}^{2}}{4\alpha}+\sqrt{\frac{\tilde{\beta}^{2}}{4}}\absolutevalue{\mathbf{A}\mathbf{B}})\operatorname{erfc}\quantity(\sqrt{\frac{\tilde{\beta}^{2}}{4\alpha}}+\sqrt{\alpha}\absolutevalue{\mathbf{A}\mathbf{B}})\exp[\frac{\tilde{\beta}^{2}}{4\alpha}+\tilde{\beta}\absolutevalue{\mathbf{A}\mathbf{B}}]\Bigg{\}}.\end{split}$$

$${\tilde{V}}_{\mathbf{a}B}=\text{matrixelement}{\chi }_{\mathbf{a}}\text{absolutevalue}{\mathbf{r}-\mathbf{C}}^{-1}{\tilde{\chi }}_B,$$

Full text

††journal: Advances in Quantum Chemistry

1 Introduction

In the last decade, the advent of massively parallel computational platforms and their ever-growing number of computing nodes has unveiled new horizons for studying quantum systems. It is now widely recognized that there is an imperative need to develop methods that take full advantages of these new supercomputer architectures and scale up to an arbitrary number of cores. A class of methods known to scale up nicely are stochastic approaches, and especially quantum Monte Carlo (QMC) methods which are steadily becoming the go-to computational tool for reaching high accuracy in large-scale problems (see, for example [1, 2, 3, 4, 5]). In practice, to make QMC feasible for large systems, it is essential to resort to accurate trial wave functions leading both to an efficient sampling of the configuration space and to low energy fluctuations. A precious guide to build up such functions is to take into account the universal features known about the exact many-electron wave function [6, 7, 8, 9, 10, 11, 12, 13].

In standard QMC implementations, the trial wave functions are usually defined as [14, 15, 16]

[TABLE]

where $D_{I}^{\sigma}$ and $\bm{R}^{\sigma}$ are determinants and coordinates of the spin-$\sigma$ electrons, respectively. The fermionic nature of the wave function is imposed using a single- or multi-determinant expansion of Slater determinants [17, 18, 19, 20, 21, 22] made of Hartree-Fock (HF) or Kohn-Sham (KS) molecular orbitals (MOs)

[TABLE]

built as a linear combination of $N$ Gaussian basis functions $\chi_{\mu}(\bm{r})$. $J(\bm{R})$ is called the Jastrow factor and its main purpose is to catch the bulk of the dynamic electron correlation.

At short interparticle distances, the Coulombic singularity dominates all other terms and, near the two-particle coalescence point, the behaviour of the exact wave function $\Psi$ becomes independent of other details of the system [13]. In particular, early work by Kato [23, 24], and elaborations by Pack and Byers Brown [25], showed that, as one electron at $\bm{r}_{i}$ approaches a nucleus of charge $Z_{A}$ at $\bm{r}_{A}$, we have

[TABLE]

where $\left.\expectationvalue{\Psi(\bm{R})}\right|_{r_{i}=r_{A}}$ is the spherical average of $\Psi(\bm{R})$ about $\bm{r}_{i}=\bm{r}_{A}$.

To remove divergences in the local energy at the electron-nucleus (e-n) coalescence points, cusp conditions such as (3) must be satisfied. (Note that the use of pseudopotentials would also remove the divergence at the e-n coalescence points as routinely done in QMC calculations, but at the price of introducing systematic errors such as the pseudopotential approximate representation and the localization error.) These divergences are especially harmful in DMC calculations, where they can lead to a large increase of the statistical variance, population-control problems and significant biases [15].

There are two possible ways to enforce the correct e-n cusp. One way to do it is to enforce the e-n cusp within the Jastrow factor in Eq. (1). This has the disadvantage of increasing the number of parameters in $J(\bm{R})$, and their interdependence can be tricky as one must optimise the large number of linear and non-linear parameters contained in $J(\bm{R})$ via a stochastic (noisy) optimization of the energy and/or its variance. However, it is frequently done in the literature thanks to some recent progress [26, 27, 28]. Another way is to enforce the cusp within the multideterminant expansion of Eq. (1).

However, because one usually employs Gaussian basis functions [29] (as in standard quantum chemistry packages), the MOs $\phi_{i}(\bm{r})$ are cuspless, i.e.

[TABLE]

One solution would be to use a different set of basis functions [30] as, for instance, Slater basis functions [31, 32, 33]. However, they are known to be troublesome, mainly due to the difficulty of calculating multicentric two-electron integrals which require expensive numerical expansions or quadratures. Nevertheless, some authors [34] have explored using wave functions built with Slater basis functions [35] while imposing the right e-n cusp afterwards. (Note that it is also possible to enforce the correct e-n cusp during the SCF process although it is rarely done [36].) These types of calculations can be performed with an electronic structure package such as ADF [37]. However, as far as we know, it is hard to perform large-scale calculations with Slater basis functions and the virtual space is usually poorly described. Moreover, Gaussian bases are usually of better quality than Slater-based ones due to the extensive knowledge and experience gathered by quantum chemists over the last fifty years while building robust, compact and effective Gaussian basis sets [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51].

Conventional cusp correction methods usually replace the part of $\chi_{\mu}(\bm{r})$ or $\phi_{i}(\bm{r})$ close to the nuclei within a cusp-correction radius by a polynomial or a spline function which fulfils Kato cusp conditions and ensures a well-behaved local energy [52, 53, 54, 55]. For atoms, one can also substitute Gaussian core orbitals by tabulated Slater-based ones [56, 57, 58]. In the same vein, Toulouse and Umrigar have fitted Slater basis functions with a large number of Gaussian functions and replaced them within the QMC calculation [28]. However, it is hardly scalable for large systems due to its lack of transferability and the ever-growing number of primitive two-electron integrals to compute.

Here, we propose to follow a different, alternative route by augmenting conventional Gaussian basis sets with cusp-correcting Slater basis functions. Mixed Gaussian-Slater basis sets have been already considered in the past with limited success due to the difficultly of computing efficiently mixed electron repulsion integrals [59, 60, 61, 62, 63, 64, 65, 66]. However, we will show that, because of the way we introduce the cusp correction, the integrals required here are not that scary. For the sake of simplicity, we will focus on the HF formalism in the present study, although our scheme can also be applied in the KS framework.

2 Cusp-corrected orbitals

A sufficient condition to ensure that $\Phi$ fulfills the e-n cusp (3) is that each (occupied and virtual) MO $\tilde{\phi_{i}}(\bm{r})$ satisfies the e-n cusp at each nuclear position $\bm{r}_{A}$:

[TABLE]

Note that this is true only if no linear term in $r$ is introduced within the Jastrow factor. Without loss of generality, we also assume that the basis functions have been already orthogonalized via the standard procedure [67], i.e. $\innerproduct{\chi_{\mu}}{\chi_{\nu}}=\delta_{\mu\nu}$, where $\delta_{\mu\nu}$ is the Kronecker delta [68].

Here, we enforce the correct e-n cusp by adding a cusp-correcting orbital to each MO:

[TABLE]

with

[TABLE]

where $M$ is the number of nuclear centers and

[TABLE]

is a $s$-type Slater function centered on nucleus $A$ with an orbital-dependent exponent $\tilde{\alpha}_{i}$. In Eq. (6), the projector

[TABLE]

(where $\hat{I}$ is the identity operator) ensures orthogonality between $\phi_{i}(\bm{r})$ and the cusp-correcting orbital $\varphi_{i}(\bm{r})$.

It is easy to show that ensuring the right e-n cusp yields the following linear system of equations for the coefficients $c_{Ai}$:

[TABLE]

where

[TABLE]

is the Kronecker delta [68] and the explicit expression of the matrix elements $\tilde{S}_{\mu A}^{i}=\innerproduct{\chi_{\mu}}{{\tilde{\chi}_{A}}^{i}}$ is given in Appendix and

[TABLE]

Equation (10) can be easily solved using standard linear algebra packages, and provides a way to obtain a cusp-corrected orbital $\tilde{\phi_{i}}(\bm{r})$ from a given MO $\phi_{i}(\bm{r})$. For reasons that will later become apparent, we will refer to this procedure as a one-step (OS) calculation. In the next section, we are going to explain how one can optimize self-consistently the coefficients $\tilde{c}_{Ai}$.

3 Self-consistent dressing of the Fock matrix

So far, the coefficient $\tilde{c}_{Ai}$ have been set via Eq. (10). Therefore, they have not been obtained via a variational procedure as their only purpose is to enforce the e-n cusp. However, they do depend on $\phi_{i}(\bm{r}_{A})$, hence on the MO coefficients $c_{\mu i}$. We will show below that one can optimize simultaneously the coefficients $\tilde{c}_{Ai}$ and $c_{\mu i}$ by constructing an orbital-dependent effective Fock matrix.

As it is ultimately what we wish for, the key point is to assume that $\tilde{\phi_{i}}(\bm{r})$ is an eigenfunction of the Fock operator $\hat{f}$, i.e.

[TABLE]

Note that, even at convergence of a conventional HF or KS calculation, the equality (13) is never fulfilled (unless the basis happens to span the exact orbital). This under-appreciated fact has been used by Deng et al. to design a measure of the quality of a MO [69].

Next, we project out Eq. (13) over $\bra*{\chi_{\mu}}$ yielding

[TABLE]

where

[TABLE]

In the general case, because we must use basis functions with non-zero derivatives at the nucleus, finding the matrix elements $\tilde{F}_{\mu A}^{i}$ is challenging and costly. However, because we are interested in the e-n cusp, we have found that a satisfactory approximation is

[TABLE]

where

[TABLE]

and $\hat{h}$ is the core Hamiltonian. (The expression of the matrix elements $\tilde{h}_{\mu A}^{i}$ are given in Appendix.) Note that, in Eq. (16), it is important to use the same approximation for both terms ($\tilde{F}_{\mu A}^{i}\approx\tilde{h}_{\mu A}^{i}$ and $F_{\mu\nu}\approx H^{\text{c}}_{\mu\nu}$) in order to preserve the subtle balance between the two terms.

The eigenvalue problem given by Eq. (14) can be recast as

[TABLE]

where we have “dressed” the diagonal of the Fock matrix

[TABLE]

with

[TABLE]

The process is repeated until our convergence criterion is met, i.e. the largest absolute value of the elements of the commutator $\tilde{\bm{F}}^{i}\bm{P}-\bm{P}\tilde{\bm{F}}^{i}$ is lower than a given threshold, where $\tilde{\bm{F}}^{i}$ is the dressed Fock matrix [Eq. (19)] and $\bm{P}$ is the density matrix with

[TABLE]

In the remainder of this paper, we will refer to this procedure as self-consistent dressing (SCD).

Similar to the Perdew-Zunger self-interaction correction [70], the orbitals $\tilde{\phi_{i}}(\bm{r})$ are eigenfunctions of different Fock operators and therefore no longer necessarily orthogonal. Practically, we have found that the e-n cusp correction makes the cusp-corrected MOs $\tilde{\phi_{i}}$ slightly non-orthogonal. However, this is not an issue as, within QMC, one evaluates the energy via MC sampling which only requires the evaluation of the MOs and their first and second derivatives.

Obviously, as evidenced by Eq. (19), when $c_{\mu i}$ is small, the dressing of the Fock matrix is numerically unstable. Therefore, we have chosen not to dress the Fock matrix if $c_{\mu i}$ is smaller than a user-defined threshold $\tau$. We have found that a value of $10^{-5}$ is suitable for our purposes, and we use the same value for the convergence threshold. Moreover, we have found that setting [53]

[TABLE]

(where $\mathring{\phi_{i}}(\bm{r})$ corresponds to the $s$-type components of $\phi_{i}(\bm{r})$ centered at $\bm{r}_{A}$) yields satisfactory results. In the case where $\mathring{\phi_{i}}(\bm{r}_{A})=0$, the MO is effectively zero at $\bm{r}=\bm{r}_{A}$ and, therefore, does not need to be cusp corrected. As in conventional self-consistent calculations, it is sometimes useful to switch on the convergence accelerator DIIS [71, 72], and we have done so in some cases. The general skeleton of the algorithm is given in Algorithm 1.

4 Illustrative examples

4.1 Atoms

Let us illustrate the present method with a simple example. For pedagogical purposes, we have computed the wave function of the hydrogen atom within a small Gaussian basis (decontracted STO-3G basis). In Fig. 1, we have plotted the local energy associated with this wave function as well as its OS and SCD cusp-corrected versions. The numerical results are reported in Table 1. As expected, the “cuspless” local energy (red curve) diverges for small $r$ with a variational energy off by $4.3$ millihartree compared to the exact value of $-1/2$. The OS cusp-correcting procedure which introduces a Slater basis function of unit exponent (but does not re-optimise any coefficients) cures the divergence at $r=0$ and significantly improves (by roughly one order of magnitude) both the variational energy and the variance. Moreover, we observe that the long-range part of the wave function is also improved compared to the Gaussian basis set due to the presence of the Slater basis function which has the correct asymptotic decay. The SCD cusp-correcting procedure further improve upon the OS scheme, and we reach a variance lower than $10^{-8}$ after only 3 iterations. The values of the coefficients of the Gaussian and Slater functions reported in Table 1 clearly show that, as expected, the Gaussian functions are getting quickly washed away and replaced by the Slater function.

In Fig. 2, we have plotted the cuspless HF 1s core orbitals of the helium (left) and neon (right) atoms, and their cusp-corrected versions obtained with various schemes. For the \ceHe atom, we compare the cusp-corrected orbitals produced by the OS and SCD procedures. One can clearly see that the qualitative difference between the cusp-corrected orbitals is small (at least graphically). For the \ceNe atom, one cannot graphically distinguished between the two cusp-correcting schemes. Figure 3 reports the local energy of the cuspless and cusp-corrected HF wave functions as an electron is moved through the nucleus of a \ceNe atom located at the origin. (The other electrons have been positioned randomly.) The right panel of Fig. 3 corresponds to a zoom around the origin where the local energy associated with the cuspless wave function is strongly oscillatory and ultimately diverges towards $-\infty$ as $r\to 0$. We observe that the cusp-correcting algorithm removes both the divergence of the local energy at the origin but also smooths out its erratic oscillations in the neighborhood of the origin, while remaining identical to the local energy obtained with the cuspless wave function for large $r$. In particular, one can see that the node (i.e. zero) of the wave function around $r=-2$ is not significantly altered by the addition of the Slater basis function (albeit not strictly identical).

Table 2 reports the energy and the corresponding variance of the \ceHe atom computed at the VMC and DMC level. The trial wave function is the HF wave function computed in the 6-31G basis with or without cusp correction. Similar to the case of the \ceH atom discussed above, the OS and SCD schemes reduce significantly both the variational energy and the variance. The energy decreases by roughly $2.8$ and $3.1$ millihartree (compared to the uncorrected scheme) using OS and SCD, respectively. Likewise, we observe that the variance is reduced by one order of magnitude by applying our cusp-corrected schemes, the difference between OS and SCD being negligible.

4.2 Molecules

As a molecular example, we consider the beryllium hydride molecule \ceBeH2 at experimental geometry. The Gaussian basis set is Pople’s 6-31G basis and the Slater exponents have been obtained via Eq. (22). The two HF valence orbitals (a${}_{\text{g}}$ and b${}_{\text{1u}}$) of this linear molecule are depicted in Fig. 4. The black line corresponds to the difference between the cuspless and cusp-corrected orbitals magnified by one order of magnitude. Note that the cusp-correcting scheme does not correct the second MO on the \ceBe nucleus because the value of this MO is effectively zero on this center.

The corresponding local energies as an electron is moved through the nuclei (marked with thin black lines) of the \ceBeH2 molecule are represented in Fig. 5. (The other electrons have been positioned randomly.) Similarly to the results of the previous section, the cusp-correcting scheme removes the divergences of the local energy at the nuclei. Note, however, that a discontinuity appears in the local energy at the nuclear centers (see inset graphs of Fig. 5). It is well-known that these discontinuities do not lead to any problems within QMC calculations. Note also that the node of the wave function around $z=-3/4$ is significantly shifted due to the introduction of the Slater basis function.

5 Conclusion

We have introduced a procedure to enforce the electron-nucleus (e-n) cusp by augmenting conventional (cuspless) Gaussian basis sets with cusp-correcting Slater basis functions. Two types of procedure has been presented. In the one-step (OS) procedure, the coefficients of the Slater functions are obtained by ensuring the correct e-n cusp at each nucleus. We have also designed a self-consistent procedure to optimize simultaneously the coefficients of the Gaussian and Slater basis functions by diagonalization of an orbital-dependent effective Fock operator.

The same procedure could potentially be employed to correct the long-range part of the electronic density with obvious application within DFT. We are currently working on a similar methodology to enforce the electron-electron cusp in explicitly correlated wave functions. We hope to be able to report on this in the near future.

Appendix A Dressing integrals

Thanks to the approximation (16), we eschew the calculations of two-electron integrals and we only need to consider three types of one-electron integrals in order to dress the Fock matrix in Eq. (14): overlap, kinetic energy and nuclear attraction. (For a OS calculation, only the former is mandatory.) Their particularity is that they are “mixed” integrals as they involve one Gaussian function (with arbitrary angular momentum) and one (momentumless) Slater function [73, 74]. Note that one can easily generalized the present procedure and consider a (contracted) linear combination of $s$-type Slater functions.

A.1

Overlap integrals

We must find mixed Gaussian-Slater overlap integrals of the form

[TABLE]

where ${\tilde{\chi}_{B}}(\bm{r})$ is given by Eq. (8) (where we have removed the superscript $i$ for the sake of clarity) and

[TABLE]

is a primitive Gaussian function of exponent $\alpha$ and angular momentum $\bm{a}=(a_{x},a_{y},a_{z})$ centered in $\mathbf{A}=(A_{x},A_{y},A_{z})$, its total angular momentum being given by $a=a_{x}+a_{y}+a_{z}$. Contracted integrals can be obtained by a straightforward summation of the primitive integrals weighted by their contraction coefficients [73, 75, 74].

We shall start by reporting the expression of the fundamental integral $\tilde{S}_{\mathbf{0}B}$ (where $\mathbf{0}=(0,0,0)$). Higher angular momentum integrals can be obtained by differentiation with respect to the center coordinates. For example, we have

[TABLE]

Using the Gaussian representation of a Slater function

[TABLE]

we obtain

[TABLE]

where $\mathbf{A}\mathbf{B}=\mathbf{A}-\mathbf{B}$ and $\operatorname{erfc}(x)$ is the complementary error function [68]. It is easy to show that the use of the Gaussian representation (26) allows us to reduce the integral to the conventional Gaussian-type function case, which has been extensively discussed in the literature [67, 76]. One only needs to perform the last integration that can be easily performed using a computer algebra system such as Mathematica [77].

A.2

Kinetic energy integrals

The same technique is applied to the kinetic energy integral

[TABLE]

which yields

[TABLE]

A.3

Nuclear attraction integrals

For the nuclear attraction integrals of the form

[TABLE]

in addition to the Gaussian representation of the Slater function (see Eq. (31)), we also use the well-known Gaussian representation of the Coulomb operator:

[TABLE]

We obtain

[TABLE]

where $F_{0}(t)$ is the Boys function [78, 79, 80]. To the best of our knowledge, this expression cannot be integrated further (expect in some particular cases) but it can be efficiently evaluated by numerical quadrature using the Gauss-Legendre rule. In the case of the nuclear attraction integrals, due to the form of the integrand in Eq. (32), we use the conventional Gaussian recurrence relations to evaluate higher angular momentum integrals. These involve the evaluation of the generalized Boys functions $F_{m}(t)$ which can be computed efficiently using well-established algorithms [78, 79, 80]. We refer the reader to Ref. [76] for more details.