All-electron quantum Monte Carlo with Jastrow single determinant Ansatz: application to the sodium dimer

TL;DR

This study demonstrates that combining variational and lattice regularized diffusion Monte Carlo methods with a Jastrow single determinant ansatz yields highly accurate potential energy surfaces for the sodium dimer, closely matching experimental data.

Contribution

The paper introduces a novel approach using a Jastrow Antisymmetrized Geminal Product ansatz with optimized basis functions to achieve chemical accuracy in quantum Monte Carlo calculations of diatomic molecules.

Findings

LRDMC reaches ~1 kcal/mol accuracy for dissociation energy

Optimized basis reduces computational effort and statistical fluctuations

Method achieves very good agreement with experimental vibrational frequencies

Abstract

In this work, we report potential energy surfaces (PESs) of the sodium dimer calculated by variational (VMC) and lattice regularized diffusion Monte Carlo (LRDMC). The VMC calculation is accurate for determining the equilibrium distance and the qualitative shape of the experimental PES. Remarkably, after the application of the LRDMC projection to this single determinant ansatz, namely the Jastrow Antisymmetrized Geminal Product (JAGP), chemical accuracy (~ 1kcal/mol) is reached, and the obtained dissociation energy, equilibrium internuclear distance, and harmonic vibrational frequency are in very good agreement with the experimental ones. This outcome crucially depends on the quality of the optimization used to determine the best possible trial function within the chosen ansatz. The strategy adopted in this work is to minimize the variational energy by initializing the trial function with the DFT single determinant ansatz expanded exactly in the same atomic basis used for the corresponding VMC and LRDMC calculations. This atomic basis is ad-hoc reshaped for QMC calculations. Indeed, we multiply the standard Gaussian type atomic orbitals by a one-body Jastrow factor, satisfying in this way the electron-ion cusp conditions. This allows us to use a very small basis almost converged in the complete basis set limit, by reducing the computational effort as well as the statistical fluctuations on the total energy. In order to achieve these important advantages, we have defined a very efficient DFT algorithm in the mentioned basis, by estimating the corresponding matrix elements on a mesh, and by using a much finer mesh grid in the vicinity of nuclei.

Tables

Table 1. Table 1: Ground state energies of the sodium atom obtained by HF, VMC, and LRDMC, and the estimated exact energy. GF denotes the guiding function.

Grid	Method	Energy (Ha)	Correlation (%)
Reference	HF	-161.8589 333Reference Schäfer et al., 1992.	0
	VMC-JSD	-162.20717(33) 444Reference Buendía et al., 2006.	88.0(1)
	VMC-JAGP	-162.1434(7) 555Reference Casula and Sorella, 2003.	71.9(2)
	DMC (GF=JAGP)	-162.2370(1) 666Reference Casula and Sorella, 2003.	95.6(3)
	DMC (GF=STO-HF)	-162.23966(22) 777Reference Nemec et al., 2010.	96.2(2)
	Exact	-162.2546 888Reference Chakravorty et al., 1993.	100
	DFT-LDA	-161.4161	-
	VMC-JDFT	-162.20442(21)	87.32(5)
	VMC-JSD	-162.21423(16)	89.80(4)
Fine grid	VMC-JAGP	-162.22045(16)	91.37(4)
(0.02 Bohr)3	LRDMC (GF=JDFT)	-162.23924(23)	96.12(6)
	LRDMC (GF=JSD)	-162.24017(57)	96.35(14)
	LRDMC (GF=JAGP)	-162.24249(16)	96.94(4)
Fine grid	DFT-LDA	-165.5715	-
(0.05 Bohr)3	VMC-JDFT	-154.554(11)	-
Fine grid	DFT-LDA	Unstable	-
(0.10 Bohr)3	VMC-JDFT	Unstable	-
	DFT-LDA	-162.9714	-
	VMC-JDFT	-162.20151(21)	86.58(5)
	VMC-JSD	-162.21474(17)	89.93(4)
Cubic interpolation	VMC-JAGP	-162.22079(16)	91.46(4)
(0.10 Bohr)3 $+$ (0.01 Bohr)3	LRDMC (GF=JDFT)	-162.23764(24)	95.71(6)
	LRDMC (GF=JSD)	-162.24078(20)	96.51(5)
	LRDMC (GF=JAGP)	-162.24250(16)	96.94(4)
Cubic interpolation	DFT-LDA	-161.5461	-
(0.05 Bohr)3 $+$ (0.01 Bohr)3	VMC-JDFT	-162.20368(18)	87.13(5)
	LRDMC (GF=JDFT)	-162.23779(24)	95.75(6)
Cubic interpolation	DFT-LDA	-173.13140	-
(0.20 Bohr)3 $+$ (0.01 Bohr)3	VMC-JDFT	-162.18028(28)	81.22(7)
Linear interpolation	DFT-LDA	-174.67577	-
(0.20 Bohr)3 $+$ (0.01 Bohr)3	VMC-JDFT	-162.17665(29)	80.30(7)

Table 2. Table 2: Summary of the VMC-JAGP and LRDMC calculations of the sodium dimer.

$d_{\text{Na-Na}}$ (Å)	$E_{\text{VMC-AGP}}$ (Ha/Na2)	$E_{\text{LRDMC}}$ (Ha/Na2)
1.8	-324.37708(34)	-324.43673(43)
1.9	-324.39731(25)	-324.45269(45)
2.0	-324.40784(33)	-324.46654(44)
2.2	-324.42894(27)	-324.48354(46)
2.4	-324.44191(22)	-324.49596(44)
2.5	-324.44593(30)	-324.50076(44)
2.7	-324.45197(25)	-324.50514(45)
2.9	-324.45554(25)	-324.50944(44)
3.0789	-324.45664(23)	-324.51033(44)
3.1	-324.45668(21)	-324.50982(43)
3.3	-324,45505(23)	-324.50911(44)
3.5	-324.45303(29)	-324.50616(44)
5.0	-324.44120(28)	-324.48999(43)
6.0	-324.44028(28)	-324.48502(42)
7.0	-324.44054(23)	-324.48538(42)
8.0	-324.44156(20)	-324.48511(43)
9.0	-324.44202(30)	-324.48443(42)
10.0	-324.44219(30)	-324.48459(43)
2 $\times$ Na atom	-324.44158(32)	-324.48500(31)

Table 3. Table 3: Summary of the obtained dissociation energies, equilibrium distances and harmonic vibrational frequencies in this work. Those obtained by previous ab initio calculations and the experimental values are also listed.

Method	$D_{\text{e}}$ (mHa)	$d_{\text{eq}}$ (Å)	${\omega }_{\text{e}}$ (cm-1)
UHF	3.20	3.60	78.93
UCCSD(T)	26.49	3.179	154.79
QCISD 999Reference Ben-Hai et al., 2007. $D_{\text{e}}$ = 0.72193 (eV), $d_{\text{eq}}$ = 0.31813 (nm), and ${\omega }_{\text{e}}$ = 151.63 (cm-1).	26.53	3.181	151.63
VMC	13.58(10)	3.050(6)	155.5(1.7)
LRDMC	25.28(43)	3.083(11)	163.4(3.4)
Full valence CI 101010Reference Magnier et al., 1993. $D_{\text{e}}$ = 5892 (cm-1), $d_{\text{eq}}$ = 5.83 (Bohr), and ${\omega }_{\text{e}}$ = 159.1 (cm-1).	26.85	3.09	159.1
Experiment 111111Reference Liu et al., 1995. $D_{\text{e}}$ = 6022.6 (cm-1), $d_{\text{eq}}$ = 5.82 (Bohr), and ${\omega }_{\text{e}}$ = 159.1 (cm-1).	27.44	3.08	159.1
Experiment 121212Reference Huber, 2013. $D_{\text{e}}$ = 0.7298 (eV), $d_{\text{eq}}$ = 0.3079 (nm), and ${\omega }_{\text{e}}$ = 159.12 (cm-1).	26.82	3.079	159.12

Table 4. Table 4: Dissociation energies of the sodium dimer at the experimental equilibrium distance ( d Na-Na subscript 𝑑 Na-Na d_{\text{Na-Na}} = 3.0789 Å). GF denotes the guiding function.

Method	GF	$E_{\text{dimer-2atoms}}$ (mHa)
DMC	STO-HF	23.87(57)13131314.981$\pm$0.357 kcal/mol. See the supplementary material of ref. Nemec et al.,2010
LRDMC	JDFT	23.23(60)
LRDMC	JSD	23.47(50)
LRDMC	JAGP	25.34(54)

Table 5. Table 5: The absolute energies of the sodium atom and dimer obtained by LRDMC (GF=JDFT, JAGP) and an experiment. The nodal surface errors (NS errors) of the absolute energy and the binding energies due to the error cancellations are also shown. The exact energy of the sodium atom and the exact binding energy are cited from Ref.  Chakravorty  et al. , 1993 and   Huber, 2013 , respectively. The exact energy of the sodium dimer was calculated these values.

	LRDMC (JDFT)		LRDMC (JAGP)		Experiment
	Energy (Ha)	NS error (mHa)	Energy (Ha)	NS error (mHa)	Exact energy (Ha)
2 atoms	-324.4753(48)	33.92	-324.4850(31)	24.21	-324.5092141414Reference Chakravorty et al., 1993.
Dimer	-324.4985(36)	37.51	-324.5103(44)	25.69	-324.5360
$D_0$ (mHa)	23.23(60)		25.34(54)		26.82151515Reference Huber, 2013.
NS error (mHa)	3.59		1.48		-

Equations

$$J_1\left({\vec{r}}_1,\dots {\vec{r}}_N\right)=\exp \left(\sum _{i,I,l}{g}_{I,l}{\chi }_{I,l}^J\left({\vec{r}}_i\right)\right)\cdot \prod _i{\tilde{J}}_1\left({\vec{r}}_i\right),$$

$${\tilde{J}}_1\left(\vec{r}\right)=\exp \left(\sum _I-{\left(2Z_I\right)}^{3/4}u\left({2Z_I}^{1/4}\left|\vec{r}-{\vec{R}}_I\right|\right)\right),$$

$$u\left(r\right)=\frac{b}{2}\left(1-e^{-r/b}\right).$$

$$J_2\left({\vec{r}}_1,\dots {\vec{r}}_N\right)=\exp \left(\sum _{i<j}v\left(r_{i,j}\right)\right),$$

$$v\left(r\right)=\frac{r}{2}{\left(1-Fr\right)}^{-1}.$$

$$J_{3/4}\left({\vec{r}}_1,\dots {\vec{r}}_N\right)=\exp \left(\sum _{i<j}{\Phi }_J\left({\vec{r}}_i,{\vec{r}}_j\right)\right)$$

$${\Phi }_J\left({\vec{r}}_i,{\vec{r}}_j\right)=\sum _{l,m,a,b}{g}_{l,m}^{a,b}{\chi }_{a,l}^J\left({\vec{r}}_i\right){\chi }_{b,m}^J\left({\vec{r}}_j\right),$$

$${\Psi }_{\text{AGP}}\left({\vec{r}}_1,\dots ,{\vec{r}}_N\right)=\det \left({\Phi }_{\text{AGP}}\left({\vec{r}}_i,{\vec{r}}_j\right)\right)$$

$${\Phi }_{\text{AGP}}\left({\vec{r}}_i,{\vec{r}}_j\right)=\sum _{l,m,a,b}{\lambda }_{a,b}^{l,m}{\varphi }_{a,l}\left({\vec{r}}_i\right){\varphi }_{b,m}\left({\vec{r}}_j\right)$$

$${\Phi }_{\text{AGP}}\left({\vec{r}}_i,{\vec{r}}_j\right)=\sum _k^L{\lambda }_k{\tilde{\psi }}_k\left({\vec{r}}_i\right){\tilde{\psi }}_k\left({\vec{r}}_j\right),$$

$${\tilde{\psi }}_k\left(\vec{r}\right)=\sum _a^M\sum _l^{L_a}{c}_{a,l}^k{\varphi }_{a,l}\left(\vec{r}\right),$$

$$H=-\frac{1}{2}\sum _i^N{\Delta }_i+V\left(\vec{x}\right)+{\sum }_{I<J}\frac{Z_I{Z}_J}{\left|{\vec{R}}_I-{\vec{R}}_J\right|},$$

$$V_{\text{ee}}\left(\vec{x}\right)={\sum }_{i<j}\frac{1}{\left|{\vec{r}}_i-{\vec{r}}_j\right|},$$

$$V_{\text{ei}}\left(\vec{x}\right)=-\sum _i{\nu }_{\text{ei}}\left({\vec{r}}_i\right)=-\sum _i\sum _I\frac{Z_I}{\left|{\vec{r}}_i-{\vec{R}}_I\right|},$$

$${\Delta }_i\approx {\Delta }_i^a={\Delta }_i^{a,p}+{\Delta }_i^{a^{\prime },1-p},$$

$$\begin{array}{rl}{\Delta }_i^{a,p}f\left(x_i,y_i,z_i\right)& =\frac{1}{a^2}[p\left(x_i+a/2\right)\left(f\left(x_i+a\right)-f\left(x_i\right)\right)\\ & +p\left(x_i-a/2\right)\left(f\left(x_i-a\right)-f\left(x_i\right)\right)]+x_i\leftrightarrow y_i\leftrightarrow z_i.\end{array}$$

$$p\left(\vec{r}\right)=\exp \left(-4{\left|\vec{r}-{\vec{R}}_{\text{c}}\right|}^2\right),$$

$${\nu }_{\text{ei}}\left({\vec{r}}_i\right)\to {\nu }_{\max ,i}^a\left(\vec{x}\right)=\mathrm{Max}\left[{\nu }_{\text{zv},i}^a\left(\vec{x}\right),{\nu }_{\text{ei}}\left({\vec{r}}_i\right)\right],$$

$${\nu }_{\text{zv},i}^a\left(\vec{x}\right)={\nu }_{\text{ei}}\left({\vec{r}}_i\right)+\frac{\left({\nabla }_{i,a}^2-{\nabla }_i^2\right){\Psi }_{\text{G}}\left(\vec{x}\right)}{2{\Psi }_{\text{G}}\left(\vec{x}\right)},$$

$${\nu }_{\text{ei}}\left({\vec{r}}_i\right)=\sum _I\frac{Z_I}{\text{Max}\left(\left|{\vec{r}}_i-{\vec{R}}_I\right|,a\right)}$$

$$⟨x|H^a|{\Psi }_{\text{G}}⟩=\sum _{x^{\prime }}{H}_{x,x^{\prime }}^a⟨x^{\prime }|{\Psi }_{\text{G}}⟩$$

$$H_{x,x^{\prime}}^{\text{FN},\gamma}=\left\{\begin{array}[]{l}{H_{x,x}}+\left({1+\gamma}\right){\mathcal{V}_{{\rm{SF}}}}\left({\vec{x}}\right)\,\,\,\,{\rm{for}}\,\,\,\,x^{\prime}=x,\\ {H_{x,x^{\prime}}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\rm{for}}\,\,\,\,x^{\prime}\neq x,{s_{x,x^{\prime}}}<0,\\ -\gamma{H_{x,x^{\prime}}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\rm{for}}\,\,\,\,x^{\prime}\neq x,{s_{x,x^{\prime}}}>0,\end{array}\right.$$

$$E_{\mathrm{MA}}=\frac{⟨{\Psi }_{\text{G}}|H_{\text{FN}}^a|{{\rm Y}}_0⟩}{⟨{\Psi }_{\text{G}}|{{\rm Y}}_0⟩}$$

$$H_{i,j}^{a,b}=\int d\vec{r}{\tilde{\varphi }}_j^b\left(\vec{r}-{\vec{R}}_b\right)\hat{H}\left(\vec{r}\right){\tilde{\varphi }}_i^a\left(\vec{r}-{\vec{R}}_a\right),$$

$$S_{i,j}^{a,b}=\int d\vec{r}{\tilde{\varphi }}_j^b\left(\vec{r}-{\vec{R}}_b\right){\tilde{\varphi }}_i^a\left(\vec{r}-{\vec{R}}_a\right),$$

$${\tilde{\varphi }}_j^b\left(\vec{r}-{\vec{R}}_b\right)={\varphi }_j^b\left(\vec{r}-{\vec{R}}_b\right){\tilde{J}}_1\left(\vec{r}\right),$$

$$\underset{\vec{r}\to {\vec{R}}_b}{\lim }\frac{\nabla {\tilde{\varphi }}_j^a}{{\tilde{\varphi }}_j^a}=-Z_b\frac{\vec{r}-{\vec{R}}_b}{|\vec{r}-{\vec{R}}_b|}$$

$$\hat{H}\left(\vec{r}\right)=\hat{T}+V_{\text{H}}\left(\vec{r}\right)+V_{\text{ext}}\left(\vec{r}\right)+V_{\text{XC}}\left(\vec{r}\right),$$

$$\nabla V_{\text{H}}\left(\vec{r}\right)=-4\pi \rho \left(\vec{r}\right),$$

$$V_{\text{H}}\left(x_i,y_j\right)$$

Full text

I Introduction

First principle quantum Monte Carlo (QMC) techniques, such as variational quantum Monte Carlo (VMC) and diffusion quantum Monte Carlo (DMC) are among the state-of-art numerical methods used to obtain highly accurate many-body wavefunctions Foulkes et al. (2001). Recent developments in QMC enable us to calculate not only the ground state energy but also vibrational frequencies Luo et al. (2014); Zen et al. (2013) and excited states Mussard et al. (2018); Hunt et al. (2018), as well as to study phase diagrams of materials Mazzola et al. (2018) and determine quantitative properties of a metal-insulator transition Stella et al. (2011) or an excitonic behavior Varsano et al. (2017). Although QMC has been mainly applied to model compounds such as atoms or small molecules due to the large computational cost, it should be much more feasible for ”real materials” (e.g. protein, surface, glass, etc.) in the near future, because the QMC algorithm scales very well with the number $N$ of electrons -at most $N^{4}$- and sustains almost ideal scaling in massively parallel architectures Foulkes et al. (2001).

In order to apply QMC for ”real materials”, it is convenient to replace core electrons with pseudopotentials because they have a little effect on chemical properties, and their replacement can reduce the QMC computational cost Ceperley (1986); Hammond et al. (1987); Ma et al. (2005) by a factor proportional to Z5.5∼6.5, where Z is the atomic number. Nevertheless, all-electron calculations are important as they represent useful benchmarks for highly accurate methods, removing the problem to work with a very accurate pseudopotential, though several progress have been made recently Trail and Needs (2015); Krogel et al. (2016); Trail and Needs (2017); Bennett et al. (2017, 2018); Annaberdiyev et al. (2018). Unfortunately, within QMC, all-electron calculations are rarely applied for atoms of large atomic number, mainly because they are too much computationally demanding at least in the simplest formulation of the VMC and DMC algorithms. Indeed, some progress has been obtained in VMC by considering more sophisticated trial moves in the Metropolis algorithm. Umrigar et al.Umrigar (1993); Stedman et al. (1998) have proposed an accelerated Metropolis method to reduce fluctuations in VMC of full-core atoms, wherein electrons close to the atomic cores are displaced much more slowly than those in the valence region. Analogously in DMC  Foulkes et al. (2001) the time step is decreased only around nuclei, by improving the efficiency of the algorithm as compared with a standard all-electron calculation with a single very small time step. Moreover, a very accurate trial wave function is necessary for a reasonably efficient QMC all-electron simulation because otherwise large absolute values of kinetic and potential energies around nuclei usually induce large statistical fluctuations. For the latter problem, in this work, we determine the trial function by using an appropriate density-functional theory (DFT) algorithm working in the atomic basis set that exactly satisfies the electron-ion cusp condition. Indeed, we have experienced that the large fluctuations in the corresponding VMC and DMC calculations can be substantially reduced by adopting a single determinant trial wave function that is particularly accurate in the vicinity of the nuclei.

In order to obtain a good trial wave function, the Kohn-Sham Hamiltonian matrix elements involving a rapidly varying electron density (wavefunction) in the vicinity of nuclei should be accurately evaluated in DFT. The matrix elements are composed of integrals of the kinetic, the Coulomb (Hartree and electron-ion), and the exchange-correlation (XC) parts Martin (2004). When the wave function is expanded using Gaussian-type orbitals (GTOs), integrals of the kinetic and the Coulomb parts can be done analytically Szabo and Ostlund (1982). However, the integral should be done numerically if the wave function is expanded in Slater-type orbitals (STOs), as well as in our modified GTO basis. In this case, Poisson’s equation is solved to determine the Hartree potential by integrating the Coulomb kernel over the electron density calculated at each point of the mesh. This scheme is employed in our TurboRVB code Sorella . If Poisson’s equation is solved in real space (e.g. the finite element method Beck (2000)), the integral can be evaluated without no other approximation than the finite mesh. Moreover, an arbitrary fine mesh grid can be used in the vicinity of nuclei, within the so-called multigrid approach Beck (2000). On the other hand, the multigrid approach is not easily implemented when the very efficient fast Fourier transform (FFT) is used to solve Poisson’s equation. In this case, the grid should be uniform in all the space, which increases a computational cost. This drawback can be solved by the so-called pseudo-charge method within the LAPW technique Singh (1994); Martin (2004), but the corresponding implementation is rather involved. In this work, we determine the Hartree potential with standard FFT convolution on a coarse mesh, then interpolate these values on a much finer mesh in the vicinity of nuclei. We show that this is enough to determine a good trial function that can be used as a suitable starting guess for QMC energy optimization. Although the DFT energy obtained with the above approximation is not exactly consistent with the one corresponding to a very dense uniform mesh, the VMC energies and the variances of the obtained initial trial wave functions are almost indistinguishable each other. We emphasize that this is just due to the simple and efficient interpolation scheme of the Hartree potential that we have introduced in this work.

This method is applied to the sodium dimer, which has been extensively studied both experimentally Wood and Hackett (1909); Fredrickson (1929); Demtröder et al. (1969); Demtröder and Stock (1975); Kusch and Hessel (1978); Verma et al. (1983); Kaminsky (1977); Qi et al. (2007) and theoretically Magnier et al. (1993); Ho et al. (2000); Matsunaga and Zavitsas (2004); Maroulis (2004); Harrison and Lawson (2005); Ben-Hai et al. (2007); Musiał and Bartlett (2011); Musia (2012); Morales et al. (2012) in the past decades. Several all-electron VMC and DMC studies have been reported for various atoms and molecules, Reynolds et al. (1982, 1986); Vrbik et al. (1988); Lester Jr and Hammond (1990); Umrigar (1993); Kenny et al. (1995); Lüchow and Anderson (1996); Yoshida and Miyako (1997); Huang et al. (1997); Shlyakhter et al. (1999); Sarsa et al. (2002); Casula and Sorella (2003); Casula et al. (2004); Ma et al. (2005); Caffarel et al. (2005); Buendía et al. (2006); Nemec et al. (2010); Scemama et al. (2014); Powell and Dawes (2016) but, to our knowledge, only one paper has reported the sodium dimer Nemec et al. (2010), wherein the dissociation energy at the experimental equilibrium distance has been calculated. Moreover, the full potential energy surface (PES) and other spectroscopic properties such as harmonic vibration frequency have not been calculated using all-electron VMC and DMC yet. All-electron calculations for the sodium dimer is informative as a reference because it is known that the use of a pseudopotential sometimes induces discrepancy due to the presence of the semi-core electrons Maezono et al. (2003). We successfully calculated PESs of the sodium dimer with small statistical errors, and the obtained dissociation energy, equilibrium internuclear distance, and harmonic vibrational frequency are in very good agreement with the experimental values. The main outcome of this work is that, after the optimization of the energy, a single determinant ansatz, the so-called JAGP described in the next section, can accurately describe this very weak and challenging chemical bond, within QMC technique. This is very important because a single determinant ansatz can be extended to much larger systems, even within the computationally demanding QMC methods. On the contrary, the multireference approach would be certainly impossible in this case, because it requires a number of determinants exponentially large in the number of electrons, and a corresponding computational burden.

II Methodology

II.1 Variational and Lattice regularized Diffusion Monte Carlo

The Jastrow single determinant Ansatz, a Jastrow-Slater determinant (JSD) and Jastrow antisymmetrized geminal power (JAGP) Casula and Sorella (2003) variational wave functions, are defined by the product of two terms, namely a Jastrow $J$ and an antisymmetric part ($\Psi=J{\Psi_{\text{AGP/SD}}}$). The Jastrow term is composed of one-body, two-body and three/four-body factors ($J={J_{1}}{J_{2}}{J_{3/4}}$). The one-body and two-body factors are used to fulfill the electron-ion and electron-electron cusp conditions, respectively. The one-body Jastrow factor reads:

[TABLE]

[TABLE]

where ${\vec{r}_{i}}$ are the electron positions, $R_{I}$ are the atomic positions with corresponding atomic number $Z_{I}$, $l$ runs over single-particle orbitals ${\chi_{I,l}^{J}}$ centered on the atom $I$, and ${u\left(r\right)}$ contains a variational parameter $b$:

[TABLE]

The two-body Jastrow factor reads:

[TABLE]

where ${r_{i,j}}=\left|{{{\vec{r}}_{i}}-{{\vec{r}}_{j}}}\right|$ is the distance between two electrons, and $v\left(r\right)$ contains a variational parameter $F$:

[TABLE]

The three-body Jastrow factor reads:

[TABLE]

[TABLE]

where the indices $l$ and $m$ again indicate different orbitals centered on corresponding atoms $a$ and $b$. In the present study, the coefficients of the three/four-body Jastrow factor were set to zero in case of $a\neq b$. The antisymmetric part reads:

[TABLE]

and the geminal function is expanded over an atomic basis:

[TABLE]

where indices $l$ and $m$ indicate different orbitals centered on atoms $a$ and $b$, and $i$ and $j$ are coordinates of spin up and down electrons, respectively. The antisymmetric part can also be represented by molecular orbitals Casula et al. (2004):

[TABLE]

[TABLE]

where $M$ is the number of atoms, ${{L_{a}}}$ is the number of atomic orbitals belonging to atom $a$, ${c_{a,l}^{k}}$ are the coefficients of the atomic orbitals, and $L$ is the number of molecular orbitals. If $L$ is equal to the half of the total number of electrons ($N/2$), the antisymmetric part coincides with the Slater determinant Casula and Sorella (2003); Casula et al. (2004). In this study, the cc-pDVZ basis set taken from EMSL Basis Set Library Feller (1996); Schuchardt et al. (2007) was used for the atomic orbitals ${\phi_{a,l}}\left({{{\vec{r}}_{i}}}\right)$ of sodium both in JSD and JAGP ansatz. According to the scheme recently introduced by Mazzola et. al. Mazzola et al. (2018), the orbitals whose exponents ($Z$) are larger than 300 were cut to avoid numerical instabilities. The modified cc-pDVZ basis for the sodium was finally composed of 8s7p1d for the determinant part. The basis set was treated as an uncontracted one and the exponents are relaxed in the optimization procedure. We note that the large exponent elements removed from the basis set are taken into account implicitly by means of the one-body Jastrow term Mazzola et al. (2018) that indeed allows us to fulfill exactly the electron-ion cusp conditions. The variational JSD and JAGP wave functions were optimized using the stochastic configuration in combination with the linear method Sorella et al. (2007); Umrigar et al. (2007) that enable us to optimize thousands of parameters simultaneously even within a stochastic optimization technique.

Lattice regularized diffusion Monte Carlo (LRDMC) is a projection technique that allows a systematic improvement of the variational ansatz, yielding the corresponding one with the lowest energy and the same signs in configuration space. This energy is the so-called ”fixed-node” DMC energy and can be obtained with the standard short time discretizationUmrigar (1993)- i.e. the conventional approach- or by the so-called lattice regularization, namely by discretizing on a lattice the continuous Hamiltonian Ten Haaf et al. (1995); Calandra Buonaura and Sorella (1998); Sorella and Capriotti (2000). We summarize the method here by emphasizing some important improvements for the all-electron case studied here. The interested readers should refer to Refs. Casula et al., 2005, 2006, 2010; Becca and Sorella, 2017 for details. In LRDMC, the original continuous Hamiltonian is regularized by an approximate one ${H^{a}}$, such that ${H^{a}}\to H$ for $a\to 0$, where $a$ is the parameter used to discretize the continuous space. We consider the Hamiltonian in atomic units:

[TABLE]

where $N$ is the number of electrons, $\vec{x}$ is 3$N$ dimension electron coordination, $\vec{x}=\left({{{\vec{r}}_{1}},{{\vec{r}}_{2}},\ldots,{{\vec{r}}_{N}}}\right)$, and $V\left({\vec{x}}\right)={V_{\text{ee}}}\left({\vec{x}}\right)+{V_{\text{ei}}}\left({\vec{x}}\right)$ is the standard many-body potential, which includes the electron-electron interaction:

[TABLE]

and the electron-ion interaction:

[TABLE]

where ${{{\vec{R}}}}$ and ${{{\vec{r}}}}$ are the ionic and electron positions, respectively. The kinetic part is approximated by a finite difference form:

[TABLE]

where $\Delta_{i}^{a,p}$ is defined by a mesh size $a$ and a function $p\left({\vec{r}}\right)$:

[TABLE]

In this work, we have adopted a more convenient and simpler form for the the function $p\left({\vec{r}}\right)$ that is chosen as:

[TABLE]

where ${{\vec{R}}_{\text{c}}}$ is the position of the nucleus closest to the electron in $\vec{r}$. The function $p$ decays much faster than the original form ( $p(\vec{r})=1/\left({1+{Z^{2}}/4{{\left|{\vec{r}-{{\vec{R}}_{\rm{c}}}}\right|}^{2}}}\right)$ ) and enables us to use the larger lattice space $a^{\prime}$ in a wider valence region, because the small lattice space $a$ is used only if the electron is very close to the nucleus. The constant $a^{\prime}/a$ is set to an irrational number in order to sample all the continuous space of the original HamiltonianCasula et al. (2005). The potential term $V\left(x\right)={V_{\text{ee}}}\left(x\right)+{V_{\text{ei}}}\left(x\right)$ is also discretized by the parameter $a$ to realize a smooth convergence for $a\to 0$. The electron-electron potential is not necessarily discretized, but the electron-ion one is modified as:

[TABLE]

[TABLE]

where ${\Psi_{\text{G}}}\left({\vec{x}}\right)$ is a guiding function, and $\nabla_{i,a}^{2}=\Delta_{i}^{a}$. Although the electron-ion potential (${\nu_{\text{ei}}}\left({{{\vec{r}}_{i}}}\right)$) in the right-hand side of Eq. (18) was regularized by:

[TABLE]

to cut the Coulomb singularity at small distances, we have noticed that this regularization is not necessary within the so-called fixed-node approximation during this work. This is because $\nu_{\text{zv},i}^{a}\left({\vec{x}}\right)$ does not diverge even when the electron-ion distance is small. Therefore, the Eq. (18) ensures the removal of the singularity unless in the vicinity of the nodal surface. The fixed node approximation removes also this singularity and therefore the algorithm remains always stable even without the use of Eq.(20). Now that the Hamiltonian is discretized, the efficient lattice Green function Monte Carlo algorithm Ten Haaf et al. (1995); Calandra Buonaura and Sorella (1998); Sorella and Capriotti (2000), which is valid on a lattice model, can be applied straightforwardly:

[TABLE]

The resulting algorithm is called LRDMC. The corresponding Green function matrix elements with the important sampling are ${G_{x,x^{\prime}}}={\Psi_{\text{G}}}\left({\vec{x}^{\prime}}\right)\left({{\Lambda_{x,x^{\prime}}}-H_{x^{\prime},x}^{a}}\right)/{\Psi_{\text{G}}}\left({\vec{x}}\right)$, where $\Lambda$ is a diagonal matrix with ${{\Lambda_{x,x^{\prime}}}}=\lambda$ and $\lambda$ should be a sufficiently large to obtain the ground state. The LRDMC algorithm is as follows Casula et al. (2005): given a walker with configuration ${\vec{x}}$ and weight $w$, a new configuration $x^{\prime}\neq x$ is ${p_{x,x^{\prime}}}={G_{x,x^{\prime}}}{\rm{/}}{b_{x}}$, where ${b_{x}}=\sum\nolimits_{x^{\prime}\neq x}{{G_{x^{\prime},x}}}$ is a normalization factor to make the Green function a transition probability. The walker weight is updated by a factor $w\to w\exp\left({-{\tau_{x}}{e_{L}}\left({\vec{x},\left[{{\Psi_{\text{G}}}}\right]}\right)}\right)$, where ${\tau_{x}}=-\log\left({rand}\right)/{b_{x}}$ is a diffusion time determined by a random number $0<rand\leq 1$, and ${e_{L}}\left({\vec{x},\left[{{\Psi_{\text{G}}}}\right]}\right)=\sum\limits_{x^{\prime}}{H_{x^{\prime},x}^{a}{\Psi_{\text{G}}}\left({\vec{x}^{\prime}}\right)}/{\Psi_{\text{G}}}\left({\vec{x}}\right)$ is a local energy with the guiding function. Of course, the usual branching scheme and many walker technique can also be used. Unfortunately, the Green function cannot be made strictly positive for fermions, therefore, the fixed-node approximation should be introduced Becca and Sorella (2017). To avoid the sign problem, the Hamiltonian is modified using the spin-flip term ${\mathcal{V}_{{\rm{sf}}}}\left({\vec{x}}\right)=\sum\limits_{x^{\prime}:{s_{x,x^{\prime}}}>0}{{H_{x,x^{\prime}}}{\Psi_{\text{G}}}\left({\vec{x}^{\prime}}\right)}/{\Psi_{\text{G}}}\left({\vec{x}}\right)$:

[TABLE]

where ${s_{x,x^{\prime}}}={\Psi_{\text{G}}}\left({\vec{x}}\right){H_{x,x^{\prime}}}{\Psi_{\text{G}}}\left({\vec{x}^{\prime}}\right)$ and $\gamma\geq 0$ is a real parameter. Finally, a mixed average of the fixed-node Hamiltonian:

[TABLE]

can be calculated by the weights and local energies after sufficient number of projections, where ${\left|{{\Upsilon_{0}}}\right\rangle}$ is the ground-state wave function. It has been confirmed that the mixed average energy is consistent with the fixed-node energy of the standard DMC in the limit $a\to 0$ Casula et al. (2005).

II.2 DFT algorithm in the same basis used for QMC

The trial functions for the Jastrow single determinant ansatz were determined from DFT calculations by using a single determinant ansatz expanded exactly in the same atomic basis used for the corresponding VMC and LRDMC calculations. Within DFT, the Hamiltonian and the overlap matrix elements required for the solution of the Kohn-Sham equations are represented as:

[TABLE]

[TABLE]

where $\tilde{\phi}_{i}^{a}\left({\vec{r}-{{\vec{R}}_{a}}}\right)$ and $\tilde{\phi}_{j}^{b}\left({\vec{r}-{{\vec{R}}_{b}}}\right)$ are $i$-th, $j$-th GTO for atom $a$ and $b$ multiplied by one-body Jastrow factor, namely:

[TABLE]

where ${{\tilde{J}}_{1}}\left({\vec{r}}\right)$ is the same as in Eq. (2). This formulation allows us to use a very small basis almost converged in the complete basis set limit, by reducing in this way the computational effort as well as the statistical fluctuations on the total energy Mazzola et al. (2018). Indeed, as it is simple to show, each element of the basis set satisfies the so-called electron-ion cusp condition, namely that when $\vec{r}$ is close to any atomic position ${\vec{R}}_{b}$:

[TABLE]

for all $a,b$.

In order to construct the trial wave function efficiently, we have defined an efficient DFT algorithm in the mentioned basis, by estimating the corresponding matrix elements on a mesh, and by using a much finer mesh grid in the vicinity of nuclei. The Hamiltonian operator is composed of:

[TABLE]

where $\hat{T}$ is a kinetic operator, ${V_{\text{H}}}\left({\vec{r}}\right)$ is the Hartree (electron-electron) potential, ${V_{\text{ext}}}\left({\vec{r}}\right)$ is the electron-ion potential (that may or may not include a true external potential)  and ${V_{\text{XC}}}\left({\vec{r}}\right)$ is the exchange-correlation potential. Given that the wavefunction is expanded in atomic orbitals such as GTO, the kinetic, electron-ion, and exchange-correlation terms are readily calculated at any point in real space. On the other hand, the Hartree potential is not determined in this manner, since it can be evaluated more conveniently on a uniform grid by solving Poisson’s equation with FFT:

[TABLE]

where $\rho\left({\vec{r}}\right)$ is the charge density. Therefore, the use of FFT with a fine grid in the vicinity of nuclei necessarily involves the same fine grid in the interstitial regions where the electron density smoothly changes, which gratuitously increases the computational cost. In our implementation, we have found a good compromise between accuracy and efficiency in the following way. The Hartree potential is calculated first on a *coarse uniform * grid by solving Poisson’s equation with the FFT algorithm. In the second step, the Hartree potential is interpolated on a fine grid only in the vicinity of nuclei using standard interpolation methods such as linear or cubic. A schematic figure of the linear interpolation in the two-dimensional case is shown in Fig. 1. The values on the fine grid are interpolated using the nearest four points, namely:

[TABLE]

where ${X_{I}}$ and ${Y_{J}}$ represent coarse grid points, $x_{i=nI+{s_{x}}}$ and $y_{j=nJ+{s_{y}}}$ represent interpolated points in the vicinity of nuclei, $n$ is the ratio of the interpolated grid to the coarse one (i.e $n-1$ coincides with the number of interpolated fine points between the coarse-grid ones), and $0\leq{s_{x}},{s_{y}}<n$. The values on the fine grid is interpolated by the nearest eight points in three dimensional case. The cubic interpolation is performed using the nearest twenty-four points. As a result, the matrix elements of Hamiltonian can be evaluated by combining a coarse grid and an interpolated fine grid in the vicinity of nuclei. Notice that a similar interpolation was done in pseudopotential calculation Ono and Hirose (1999, 2005), wherein the interpolation scheme was used to evaluate inner products between wave functions and nonlocal parts of pseudopotentials. The total DFT energy corresponding to the chosen interpolation for the Hartree potential is sizably different from the one obtained with a very fine grid (namely converged). However, VMC and LRDMC energies obtained with the Kohn-Sham Slater determinants with or without the interpolation scheme proposed, are very close, indicating that our DFT scheme provides very good Kohn-Sham orbitals, despite the observed error in the DFT energy.

III Validation of the interpolation scheme

To investigate the quality of the trial wave functions obtained by the interpolation technique, ground state energies of the sodium atom were calculated using DFT, VMC, and LRDMC. DFT calculations were performed with a single fine grid or using the interpolation scheme, wherein the LDA functional developed by Perdew and Zunger Perdew and Zunger (1981) was employed. Three types of single-grid DFT calculations were performed with (0.02 Bohr)3, (0.05 Bohr)3, and (0.10 Bohr)3 grids. For comparison, three types of DFT calculations using the cubic interpolation method were performed, namely, (0.01 Bohr)3 grid was used for the core electron region, while (0.05 Bohr)3, (0.10 Bohr)3 or (0.20 Bohr)3 grids were used for the valence electron region, wherein the core-electron region, centered on the sodium atom, has a volume of (2.00 Bohr)3. The calculation using the linear interpolation method was also performed using (0.01 Bohr)3 and (0.20 Bohr)3 grids. Then, three types of VMC and LRDMC calculations were performed starting from the resultant wave functions. VMC-JDFT denotes that only the Jastrow factor was optimized using the Jastrow-Slater ansatz, namely, the nodal surface was determined by the DFT. VMC-JSD and VMC-JAGP denote that both Jastrow and determinant parts were optimized using Jastrow-Slater and Jastrow antisymmetrized geminal power ansatz, respectively. LRDMC (GF=JDFT, JSD, JAGP) denotes that the wave functions optimized using each ansatz were used for the guiding functions (GF). All results are summarized in Table  1.

In the fine grid calculations, without using the interpolation technique discussed above, a well-converged result was obtained only by using (0.02 Bohr)3 single grid. Indeed, DFT calculation with (0.05 Bohr)3 grid resulted in a much worse DFT-LDA and corresponding VMC-JDFT energies and a very coarse (0.10 Bohr)3 grid implies numerical instabilities. On the other hand, a very coarse (0.10 Bohr)3 grid is already sufficient to obtain a reasonable trial wave function when the cubic interpolation method is used, and (0.05 Bohr)3 is essentially converged. As expected, the DFT energy obtained by the interpolation method (-162.9714 Ha) is not consistent with that obtained by the very fine mesh (-161.4161 Ha) due to the approximation in the Hartree potential. However, the wave function obtained in this way can be used as a trial wave function for accurate VMC and DMC calculations because the nodal surface is almost the same as the fine-grid one: The VMC-JDFT calculations (i.e. only the amplitude is optimized) show that the VMC energy obtained by the interpolation grid (-162.20151(21) Ha) is almost the same as the fine-grid one (-162.20442(21) Ha), where the deviation is only a few mHartree. The LRDMC (GF=JDFT) calculations also show a very good agreement (-162.23764(24) Ha and -162.23924(23) Ha for the cubic interpolation and fine grid, respectively). Furthermore, our LRDMC (GF=JDFT) also reproduces the reference energy (-162.23966(22) Ha) that was obtained by an all-electron DMC (GF=STO-HF) calculation using a very large basis set (quadruple-$\xi$-fourfold-polarized: QZ4P). These results indicate that the nodal surface determined by the interpolated DFT is as good as the fine-grid and the large-basis one. Thus, the interpolation method enables us to obtain a reasonable trial wave function with a low computational cost. Notice that, this interpolation method was applied also with (0.20 Bohr)3 and (0.01 Bohr)3 double mesh grids with much worse results as far as the quality of the nodal surface and corresponding VMC energies are concerned. Nevertheless , with such a sizable error, it can be clearly appreciated that the cubic interpolation performs better than the linear one.

The wave function can be further improved by optimizing the determinant part in presence of the Jastrow factor. As shown in Table  1, VMC-JSD and VMC-AGP show lower variational energies than VMC-DFT, and LRDMC(GF=JSD, JAGP) also show lower variational energies than LRDMC(GF=JDFT) thanks to the improvement of the nodal surfaces. Remarkably our LRDMC energy corresponding to our best VMC-JAGP is very close to the estimated exact total energy, namely -162.2546 (Ha) Chakravorty et al. (1993).

IV Application to the sodium dimer

PESs of the sodium dimer were calculated by VMC-JAGP and LRDMC, by using the developed interpolation scheme, and were compared with previous experiments and calculations. First, a PES was calculated using JAGP ansatz by changing the internuclear distance from 1.8 $\AA$ to 10.0 $\AA$, then, a PES was again calculated by LRDMC starting from the optimized wave functions. The energies obtained by LRDMC for each $a$ were extrapolated by quartic polynomial fits $E\left(a\right)={E_{0}}+b{a^{2}}+c{a^{4}}$ to obtain the $a\to 0$ limit (${E_{0}}$), wherein $a$ = 0.03, 0.04, 0.05, 0.06, 0.07, and 0.08 are employed (Fig. 2) in all these calculations. The VMC-JAGP and LRDMC energies are summarized in Table 2, and the obtained PESs are shown in Fig. 3. The PESs obtained by HF, MP2, CCSD(T) 111 The aug-cc-VQZ basis sets were used for these calculations.

calculated using Gaussian 09, Revision E.01 Frisch et al. (2013), and the experimental values reported by Verma et. al  Verma et al. (1983) are also shown in Fig. 3 for comparison. Notice that the energy of the molecule at large distance coincides with twice the energy of a single atom, namely the size consistency is perfectly fulfilled within VMC-JAGP and LRDMC calculations (see the bottom of Table 2).

We have analyzed the PESs, by using the simple analytic Murrell-Sorbie (MS) function Sorbie and Murrell (1975) that has been widely used for describing PES of neutral dimers,

[TABLE]

where $D_{\text{e}}$ is the dissociation energy without the zero point vibration energy (ZPVE), $\rho=d-d_{\text{eq}}$, $d$ is the internuclear distance between the sodium atoms, $d_{\text{eq}}$ is the equilibrium internuclear distance, and ${a_{1}}$, ${a_{2}}$, and ${a_{3}}$ are fitting parameters. The $D_{\text{e}}$, $d_{\text{eq}}$, ${a_{1}}$, ${a_{2}}$, and ${a_{3}}$ were determined using scipy.optimize.curve_fit module implemented in the Python SciPy library Jones et al. . Then, a harmonic vibration frequency (${\omega_{\text{e}}}$ cm*-1*) was calculated according to the following relation Ben-Hai et al. (2007):

[TABLE]

where $c$ is the light velocity, and $\mu$ is the reduced mass. The obtained values are summarized in Table 3.

The VMC calculation reproduces the qualitative shape of the experimental PES (Fig. 3), and is accurate for determining the equilibrium distance and the harmonic frequency (Table 3). However, it is not enough for obtaining the accurate dissociation energy and for reproducing the binding character in the range of 5.0 $\AA$ and 8.0 $\AA$ (i.e. showing higher energies than the dissociation limit). Notice that the VMC-JAGP is usually enough to reproduce almost correct binding energies for the second-row dimers Marchi et al. (2009). This suggests that DMC is extremely important for molecules of large atomic number and that our Jastrow is not enough accurate for describing this weak chemical bond.

The experimental PES is accurately described after the application of the LRDMC projection to the JAGP (Fig. 3). The equilibrium distance ($d_{\text{eq}}$ = 3.083(11) Å) and the harmonic frequency (${\omega_{\text{e}}}$ = 163.3(3.4) cm*-1*) obtained by our LRDMC calculations are well consistent with the experimental values ($d_{\text{eq}}$ = 3.08 Å, ${\omega_{\text{e}}}$ = 159.1 cm*-1*, and $d_{\text{eq}}$ = 3.079 Å, ${\omega_{\text{e}}}$ = 159.12 cm*-1* cited from Ref. Liu et al., 1995 and  Huber, 2013, respectively). The dissociation energy ($D_{\text{e}}$ = 25.28(43) mHa) is also consistent with the experimental ones ($D_{\text{e}}$ = 26.82 mHa Huber (2013) and $D_{\text{e}}$ = 27.44 mHa Liu et al. (1995)), and several theoretical works, such as coupled cluster calculations ($D_{\text{e}}$ = 26.49 mHa for CCSD(T) and $D_{\text{e}}$ = 26.53 mHa for QCISD Ben-Hai et al. (2007)) and full valence configuration interaction calculation ($D_{\text{e}}$ = 26.85 mHa Magnier et al. (1993)), converged within the chemical accuracy ($\sim$ 1kcal/mol $\approx$ 1.6 mHa).

Although the deviation of the dissociation energy is small enough, it is worth discussing how to obtain a more accurate result. Nemec et. al argued in their work Nemec et al. (2010) that the deviation was due to insufficient nodal error cancellation between the atoms and the dimer. They reported that the dissociation energy of the sodium dimer was underestimated ($D_{\text{e}}$ = 23.87(57) mHa at $d$ = 3.0789 Å222The original value is 14.918$\pm$0.357 kcal/mol at $d$ = 3.0789 Å. See their supplementary material.) by an all-electron DMC calculation starting from STO Nemec et al. (2010). The reduction of the error cancellation is important to obtain a more accurate result. In order to compare directly our results with the previous DMC one Nemec et al. (2010), LRDMC (GF=JDFT, JSD, JAGP) for the sodium dimer at the same distance ($d$ = 3.0789 Å) were performed. The results are shown in Fig. 4 and summarized in Table 4. LRDMC (GF=JDFT) and LRDMC (GF=JSD) give $D_{\text{e}}$ = 23.23(60) mHa and $D_{\text{e}}$ =23.47(50) mHa, respectively, which are statistically consistent with the previous DMC calculation ($D_{\text{e}}$ = 23.87(57) mHa). On the other hand, our LRDMC (GF=JAGP) greatly improves the error cancellation and provides, to our knowledge, the best binding energy ($D_{\text{e}}$ = 25.34(60) mHa) so far available within QMC techniques. Fig. 4 shows that the best variational energies are obtained when LRDMC is applied to the JAGP guiding functions both for the atom and the dimer, meaning that the corresponding nodal surfaces are better than previous calculations. Table. 5 summarizes the absolute energies of the sodium atom and dimer, and the residual errors in the absolute energies and corresponding binding energies within the fixed-node approximation. The nodal surface errors in LRDMC (GF=JDFT) are 33.92 mHa and 37.51 mHa for two sodium atoms and for the dimer, respectively. This implies 3.59 mHa smaller binding energy ($D_{\text{e}}$ = 23.23(60) mHa) than the experimental value ($D_{\text{e}}$ = 26.82 mHa) due to insufficient error cancellation. On the other hand, the nodal surface errors in LRDMC (GF=JAGP) become smaller, 24.21 mHa and 25.69 mHa for two sodium atoms and the dimer, respectively, thanks to the multiconfigurational nature of JAGP Zen et al. (2014); Genovese et al. (2019) (i.e. static correlation). This leads to a much better binding energy ($D_{\text{e}}$ = 25.34(60) mHa) because of the improvement in the error cancellation. Figure 5 shows the energy diagram and the results of the error cancellations. Compared to LRDMC (GF=JDFT), LRDMC (GF=JAGP) reduces by 9.71 mHa the nodal error for the two sodium atoms, while by 11.82 mHa for the dimer, resulting in a better error cancellation and a corresponding more accurate binding energy. While the value of Ref. Huber, 2013 is used for the exact binding energy in this discussion, the conclusion does not change when the other experimental value (e.g. 27.44 mHa Liu et al. (1995)) is employed. The fact that LRDMC (GF=JAGP) lowers the total energy more effectively in the dimer rather than in the atom indicates that the DFT nodes have some error also in the valence region because one can assume an almost exact nodal error cancellation in the core region Nemec et al. (2010). Thus, we expect that the use of more flexible wave functions such as backflow López Ríos et al. (2006) or pfaffian Bajdich et al. (2006) should further reduce the error and should lead to an almost exact error cancellation (i.e. better binding energy) and essentially exact binding energies of dimers as well as PESs.

V Summary

In this work, we report potential energy surfaces (PES) of the sodium dimer calculated by variational (VMC) and lattice regularized diffusion Monte Carlo (LRDMC). Remarkably, after the application of the LRDMC projection to the Jastrow Antisymmetrized Geminal Product (JAGP) ansatz, chemical accuracy is reached, and the obtained dissociation energy, equilibrium internuclear distance, and harmonic vibration frequency are in good agreement with the experimental ones. The trial wave functions for the VMC and LRDMC calculations were prepared by the DFT single determinant ansatz expanded exactly in the same atomic basis used for the QMC calculation, which we have conveniently devised to satisfy exactly the electron-ion cusp conditions. This allows us to use a very small basis, which is however almost converged in the complete basis set limit. In this way, it is possible to reduce the computational effort as well as the statistical fluctuations on the total energy. We have found that the improvement in the description of the electron correlation and the weak chemical bond of the sodium dimer, is mainly achieved thanks to the energy optimization strategy that we have developed in this work. For the all-electron calculation, the DFT step is computationally very demanding, at least in the convenient basis we have chosen. Therefore, we have developed an efficient DFT algorithm in the mentioned basis, by estimating the corresponding matrix elements on a mesh, and by using a much finer mesh grid only in the vicinity of the nuclei. In this way, we can have a very good description of this chemical bond and evaluate the corresponding PES with a high degree of accuracy. We believe that our work represents an important step to define a quantum Monte Carlo method that will have the same reliability and accuracy of modern quantum chemistry packages in the future, with the considerable advantage that QMC with the single determinant ansatz used in this work scales very well with the number of electrons and has an almost ideal scaling for massively parallel computations.

VI Acknowledgements

The computations in this work have been mainly performed using the facilities of Research Center for Advanced Computing Infrastructure at Japan Advanced Institute of Science and Technology (JAIST). The structure of the sodium dimer was depicted by VESTA 3 package Momma and Izumi (2011). K. Nakano is grateful for a financial support from Simons Foundation and for technical support by Prof. K. Hongo (JAIST). K. Nakano and S. Sorella are grateful for useful discussion with C. Genovese (SISSA). R. Maezono is grateful for financial supports from MEXT-KAKENHI (17H05478 and 16KK0097), from FLAGSHIP2020 (project nos. hp180206 and hp180175 at K-computer), from Toyota Motor Corporation, from I-O DATA Foundation, and from the Air Force Office of Scientific Research (AFOSR-AOARD/FA2386-17-1-4049).