\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1} \usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amsmath,amssymb,amsfonts,physics,mhchem,xspace,subfigure} \usepackage[T1]{fontenc} \usepackage[utf8]{inputenc} \usepackage{txfonts} \usepackage[ colorlinks=true, citecolor=blue, breaklinks=true ]{hyperref} \urlstyle{same} \newcommand{\ie}{\textit{i.e.}} \newcommand{\eg}{\textit{e.g.}} \newcommand{\alert}[1]{\textcolor{red}{#1}} \definecolor{darkgreen}{HTML}{009900} \usepackage[normalem]{ulem} \newcommand{\toto}[1]{\textcolor{blue}{#1}} \newcommand{\trashAS}[1]{\textcolor{blue}{\sout{#1}}} \newcommand{\titou}[1]{\textcolor{red}{#1}} \newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}} \newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}} \newcommand{\SI}{\textcolor{blue}{supplementary material}} \newcommand{\mc}{\multicolumn} \newcommand{\fnm}{\footnotemark} \newcommand{\fnt}{\footnotetext} \newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}} \newcommand{\EPT}{E_{\text{PT2}}} \newcommand{\EDMC}{E_{\text{FN-DMC}}} \newcommand{\Ndet}{N_{\text{det}}} \newcommand{\hartree}{$E_h$} \newcommand{\LCT}{Laboratoire de Chimie Th\'eorique (UMR 7616), Sorbonne Universit\'e, CNRS, Paris, France} \newcommand{\ANL}{Argonne Leadership Computing Facility, Argonne National Laboratory, Argonne, IL 60439, United States} \newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France} \DeclareMathOperator{\erfc}{erfc} \begin{document} \title{Taming the fixed-node error in diffusion Monte Carlo via range separation} %\title{Enabling high accuracy diffusion Monte Carlo calculations with % range-separated density functional theory and selected configuration interaction} \author{Anthony Scemama} \email{scemama@irsamc.ups-tlse.fr} \affiliation{\LCPQ} \author{Emmanuel Giner} \email{emmanuel.giner@lct.jussieu.fr} \affiliation{\LCT} \author{Anouar Benali} \email{benali@anl.gov} \affiliation{\ANL} \author{Pierre-Fran\c{c}ois Loos} \email{loos@irsamc.ups-tlse.fr} \affiliation{\LCPQ} \begin{abstract} By combining density-functional theory (DFT) and wave function theory (WFT) via the range separation (RS) of the interelectronic Coulomb operator, we obtain accurate fixed-node diffusion Monte Carlo (FN-DMC) energies with compact multideterminant trial wave functions. These compact trial wave functions are generated via the diagonalization of the RS-DFT Hamiltonian. In particular, we combine here short-range correlation functionals with selected configuration interaction (SCI). As the WFT method is relieved from describing the short-range part of the correlation hole around the electron-electron coalescence points, the number of determinants in the trial wave function required to reach a given accuracy is significantly reduced as compared to a conventional SCI calculation. \titou{T2: work in progress.} \end{abstract} \maketitle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Introduction} \label{sec:intro} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Solving the Schr\"odinger equation for atoms and molecules is a complex task that has kept theoretical and computational chemists busy for almost hundred years now. \cite{Schrodinger_1926} In order to achieve this formidable endeavour, various strategies have been carefully designed and implemented in quantum chemistry software packages. One of this strategies consists in relying on wave function theory and, in particular, on the full configuration interaction (FCI) method. However, FCI delivers only an approximate solution of the Schr\"odinger equation within a finite one-electron basis. This solution is the eigenpair of an approximate Hamiltonian defined as the projection of the exact Hamiltonian onto the finite many-electron basis of all possible Slater determinants generated within this finite one-electron basis. The FCI wave function can be interpreted as a constrained solution of the true Hamiltonian forced to span the restricted space provided by the one-electron basis. In the complete basis set (CBS) limit, the constraint is lifted and the exact solution is recovered. Hence, the accuracy of a FCI calculation can be systematically improved by increasing the size of the one-electron basis set. Nevertheless, its exponential scaling with the number of electrons and with the size of the basis is prohibitive for most chemical systems. In recent years, the introduction of new algorithms \cite{Booth_2009} and the revival \cite{Abrams_2005,Bytautas_2009,Roth_2009,Giner_2013,Knowles_2015,Holmes_2016,Liu_2016,Garniron_2018} of selected configuration interaction (SCI) methods \cite{Bender_1969,Huron_1973,Buenker_1974} pushed the limits of the sizes of the systems that could be computed at the FCI level. \cite{Booth_2010,Cleland_2010,Daday_2012,Chien_2018,Loos_2018a,Loos_2019,Loos_2020b,Loos_2020c} However, the scaling remains exponential unless some bias is introduced leading to a loss of size consistency. \cite{Evangelisti_1983,Cleland_2010,Tenno_2017} Diffusion Monte Carlo (DMC) is another numerical scheme to obtain the exact solution of the Schr\"odinger equation with a different constraint. In DMC, the solution is imposed to have the same nodes (or zeroes) as a given trial (approximate) wave function. Within this so-called \emph{fixed-node} (FN) approximation, the FN-DMC energy associated with a given trial wave function is an upper bound to the exact energy, and the latter is recovered only when the nodes of the trial wave function coincide with the nodes of the exact wave function. The polynomial scaling with the number of electrons and with the size of the trial wave function makes the FN-DMC method particularly attractive. In addition, the total energies obtained are usually far below those obtained with the FCI method in computationally tractable basis sets because the constraints imposed by the FN approximation are less severe than the constraints imposed by the finite-basis approximation. %However, it is usually harder to control the FN error in DMC, and this %might affect energy differences such as atomization energies. %Moreover, improving systematically the nodal surface of the trial wave %function can be a tricky job as \trashAS{there is no variational %principle for the nodes}\toto{the derivatives of the FN-DMC energy %with respect to the variational parameters of the wave function can't %be computed}. The qualitative picture of the electronic structure of weakly correlated systems, such as organic molecules near their equilibrium geometry, is usually well represented with a single Slater determinant. This feature is in part responsible for the success of density-functional theory (DFT) and coupled cluster theory. DMC with a single-determinant trial wave function can be used as a single-reference post-Hatree-Fock method, with an accuracy comparable to coupled cluster.\cite{Dubecky_2014,Grossman_2002} The favorable scaling of QMC, its very low memory requirements and its adequacy with massively parallel architectures make it a serious alternative for high-accuracy simulations of large systems. As it is not possible to minimize directly the FN-DMC energy with respect to the variational parameters of the trial wave function, the fixed-node approximation is much more difficult to control than the finite-basis approximation. The conventional approach consists in multiplying the trial wave function by a positive function, the \emph{Jastrow factor}, taking account of the electron-electron cusp and the short-range correlation effects. The wave function is then re-optimized within variational Monte Carlo (VMC) in the presence of the Jastrow factor and the nodal surface is expected to be improved. Using this technique, it has been shown that the chemical accuracy could be reached within FN-DMC.\cite{Petruzielo_2012} Another approach consists in considering the FN-DMC method as a \emph{post-FCI method}. The trial wave function is obtained by approaching the FCI with a selected configuration interaction (SCI) method such as CIPSI for instance.\cite{Giner_2013,Caffarel_2016_2} \titou{When the basis set is increased, the trial wave function gets closer to the exact wave function, so the nodal surface can be systematically improved.\cite{Caffarel_2016} WRONG} This technique has the advantage that using FCI nodes in a given basis set is well defined, so the calculations are reproducible in a black-box way without needing any expertise in QMC. But this technique cannot be applied to large systems because of the exponential scaling of the size of the trial wave function. Extrapolation techniques have been used to estimate the FN-DMC energies obtained with FCI wave functions,\cite{Scemama_2018} and other authors have used a combination of the two approaches where highly truncated CIPSI trial wave functions are re-optimized in VMC under the presence of a Jastrow factor to keep the number of determinants small,\cite{Giner_2016} and where the consistency between the different wave functions is kept by imposing a constant energy difference between the estimated FCI energy and the variational energy of the CI wave function.\cite{Dash_2018,Dash_2019} Nevertheless, finding a robust protocol to obtain high accuracy calculations which can be reproduced systematically, and which is applicable for large systems with a multi-configurational character is still an active field of research. The present paper falls within this context. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Combining CIPSI with RS-DFT} \label{sec:rsdft-cipsi} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In single-determinant DMC calculations, the degrees of freedom used to reduce the fixed-node error are the molecular orbitals on which the Slater determinant is built. Different molecular orbitals can be chosen: Hartree-Fock (HF), Kohn-Sham (KS), natural (NO) orbitals of a correlated wave function, or orbitals optimized under the presence of a Jastrow factor. The nodal surfaces obtained with the KS determinant are in general better than those obtained with the HF determinant,\cite{Per_2012} and of comparable quality to those obtained with a Slater determinant built with NOs.\cite{Wang_2019} Orbitals obtained in the presence of a Jastrow factor are generally superior to KS orbitals.\cite{Filippi_2000,Scemama_2006,HaghighiMood_2017,Ludovicy_2019} The description of electron correlation within DFT is very different from correlated methods. In DFT, one solves a mean field problem with a modified potential incorporating the effects of electron correlation, whereas in correlated methods the real Hamiltonian is used and the electron-electron interactions are considered. Nevertheless, as the orbitals are one-electron functions, the procedure of orbital optimization in the presence of the Jastrow factor can be interpreted as a self-consistent field procedure with an effective Hamiltonian,\cite{Filippi_2000} similarly to DFT. So KS-DFT can be viewed as a very cheap way of introducing the effect of correlation in the orbital parameters determining the nodal surface of a single Slater determinant. Nevertheless, even when using the exact exchange correlation potential at the CBS limit, a fixed-node error necessarily remains because the single-determinant ansätz does not have enough flexibility to describe the nodal surface of the exact correlated wave function of a generic $N$-electron system. If one wants to have to exact CBS limit, a multi-determinant parameterization of the wave functions is required. %==================== \subsection{CIPSI} %==================== Beyond the single-determinant representation, the best multi-determinant wave function one can obtain is the FCI. FCI is a \emph{post-Hartree-Fock} method, and there exists several systematic improvements between the Hartree-Fock and FCI wave functions: increasing the maximum degree of excitation of CI methods (CISD, CISDT, CISDTQ, \emph{etc}), or increasing the complete active space (CAS) wave functions until all the orbitals are in the active space. Selected CI methods take a shorter path between the Hartree-Fock determinant and the FCI wave function by increasing iteratively the number of determinants on which the wave function is expanded, selecting the determinants which are expected to contribute the most to the FCI eigenvector. At every iteration, the lowest eigenpair is extracted from the CI matrix expressed in the determinant subspace, and the FCI energy can be estimated by computing a second-order perturbative correction (PT2) to the variational energy, $\EPT$. The magnitude of $\EPT$ is a measure of the distance to the exact eigenvalue, and is an adjustable parameter controlling the quality of the wave function. Within the \emph{Configuration interaction using a perturbative selection made iteratively} (CIPSI)\cite{Huron_1973} method, the PT2 correction is computed along with the determinant selection. So the magnitude of $\EPT$ can be made the only parameter of the algorithm, and we choose this parameter as the convergence criterion of the CIPSI algorithm. Considering that the perturbatively corrected energy is a reliable estimate of the FCI energy, using a fixed value of the PT2 correction as a stopping criterion enforces a constant distance of all the calculations to the FCI energy. In this work, we target the chemical accuracy so all the CIPSI selections were made such that $\abs{\EPT} < 1$ millihartree. %================================= \subsection{Range-separated DFT} \label{sec:rsdft} %================================= Following the seminal work of Savin,\cite{Savin_1996,Toulouse_2004} the Coulomb operator entering the interelectronic repulsion is split into two pieces: \begin{equation} \frac{1}{r} = w_{\text{ee}}^{\text{sr}, \mu}(r) + w_{\text{ee}}^{\text{lr}, \mu}(r) \end{equation} where \begin{align} w_{\text{ee}}^{\text{sr},\mu}(r) & = \frac{\erfc \qty( \mu\, r)}{r}, & w_{\text{ee}}^{\text{lr},\mu}(r) & = \frac{\erf \qty( \mu\, r)}{r} \end{align} are the singular short-range (sr) part and the non-singular long-range (lr) part, respectively, $\mu$ is the range-separation parameter which controls how rapidly the short-range part decays, $\erf(x)$ is the error function, and $\erfc(x) = 1 - \erf(x)$ is its complementary version. The main idea behind RS-DFT is to treat the short-range part of the interaction within KS-DFT, and the long-range part within a WFT method like FCI in the present case. The parameter $\mu$ controls the range of the separation, and allows to go continuously from the KS Hamiltonian ($\mu=0$) to the FCI Hamiltonian ($\mu \to \infty$). To rigorously connect WFT and DFT, the universal Levy-Lieb density functional \cite{Lev-PNAS-79,Lie-IJQC-83} is decomposed as \begin{equation} \mathcal{F}[n] = \mathcal{F}^{\text{lr},\mu}[n] + \bar{E}_{\text{Hxc}}^{\text{sr,}\mu}[n], \label{Fdecomp} \end{equation} where $n$ is a one-electron density, $\mathcal{F}^{\text{lr},\mu}$ is a long-range universal density functional and $\bar{E}_{\text{Hxc}}^{\text{sr,}\mu}$ is the complementary short-range Hartree-exchange-correlation (Hxc) density functional. \cite{Savin_1996,Toulouse_2004} One obtains the following expression for the ground-state electronic energy \begin{equation} \label{min_rsdft} E_0= \min_{\Psi} \qty{ \mel{\Psi}{\hat{T}+\hat{W}_\text{{ee}}^{\text{lr},\mu}+\hat{V}_{\text{ne}}}{\Psi} + \bar{E}^{\text{sr},\mu}_{\text{Hxc}}[n_\Psi] }, \end{equation} with $\hat{T}$ the kinetic energy operator, $\hat{W}_\text{ee}^{\text{lr},\mu}$ the long-range electron-electron interaction, $n_\Psi$ the one-electron density associated with $\Psi$, and $\hat{V}_{\text{ne}}$ the electron-nucleus potential. The minimizing multideterminant wave function $\Psi^\mu$ can be determined by the self-consistent eigenvalue equation \begin{equation} \label{rs-dft-eigen-equation} \hat{H}^\mu[n_{\Psi^{\mu}}] \ket{\Psi^{\mu}}= \mathcal{E}^{\mu} \ket{\Psi^{\mu}}, \end{equation} with the long-range interacting Hamiltonian \begin{equation} \label{H_mu} \hat{H}^\mu[n_{\Psi^{\mu}}] = \hat{T}+\hat{W}_{\text{ee}}^{\text{lr},\mu}+\hat{V}_{\text{ne}}+ \hat{\bar{V}}_{\text{Hxc}}^{\text{sr},\mu}[n_{\Psi^{\mu}}], \end{equation} where $\hat{\bar{V}}_{\text{Hxc}}^{\text{sr},\mu}$ is the complementary short-range Hartree-exchange-correlation potential operator. Once $\Psi^{\mu}$ has been calculated, the electronic ground-state energy is obtained as \begin{equation} \label{E-rsdft} E_0= \mel{\Psi^{\mu}}{\hat{T}+\hat{W}_\text{{ee}}^{\text{lr},\mu}+\hat{V}_{\text{ne}}}{\Psi^{\mu}}+\bar{E}^{\text{sr},\mu}_{\text{Hxc}}[n_{\Psi^\mu}]. \end{equation} Note that, for $\mu=0$, \titou{the long-range interaction vanishes}, \ie, $w_{\text{ee}}^{\text{lr},\mu=0}(r) = 0$, and thus RS-DFT reduces to standard KS-DFT and $\Psi^\mu$ is the KS determinant. For $\mu\to\infty$, the long-range interaction becomes the standard Coulomb interaction, \ie, $w_{\text{ee}}^{\text{lr},\mu\to\infty}(r) = r^{-1}$, and thus RS-DFT reduces to standard WFT and $\Psi^\mu$ is the FCI wave function. %%% FIG 1 %%% \begin{figure*} \centering \includegraphics[width=0.7\linewidth]{algorithm.pdf} \caption{Algorithm showing the generation of the RS-DFT wave function $\Psi^{\mu}$ starting from $\Psi^{(0)}$. The outer (macro-iteration) and inner (micro-iteration) loops are represented in red and blue, respectively. The steps common to both loops are represented in purple. DIIS extrapolations of the one-electron density are introduced in both the outer and inner loops in order to speed up convergence of the self-consistent process.} \label{fig:algo} \end{figure*} %%% %%% %%% %%% Hence, range separation creates a continuous path connecting smoothly the KS determinant to the FCI wave function. Because the KS nodes are of higher quality than the HF nodes, we expect that using wave functions built along this path will always provide reduced fixed-node errors compared to the path connecting HF to FCI which consists in increasing the number of determinants. We follow the KS-to-FCI path by performing FCI calculations using the RS-DFT Hamiltonian with different values of $\mu$. Our algorithm, depicted in Fig.~\ref{fig:algo}, starts with a single- or multi-determinant wave function $\Psi^{(0)}$ which can be obtained in many different ways depending on the system that one considers. One of the particularity of the present work is that we use the CIPSI algorithm to perform approximate FCI calculations with the RS-DFT Hamiltonian $\hat{H}^\mu$. \cite{GinPraFerAssSavTou-JCP-18} This provides a multi-determinant trial wave function $\Psi^{\mu}$ that one can ``feed'' to DMC. In the outer (macro-iteration) loop (red), at the $k$th iteration, a CIPSI selection is performed to obtain $\Psi^{\mu\,(k)}$ with the RS-DFT Hamiltonian $\hat{H}^{\mu\,(k)}$ parameterized using the current one-electron density $n^{(k)}$. At each iteration, the number of determinants in $\Psi^{\mu\,(k)}$ increases. One exits the outer loop when the absolute energy difference between two successive macro-iterations $\Delta E^{(k)}$ is below a threshold $\tau_1$ that has been set to $10^{-3}$ hartree in the present study which is consistent with the CIPSI threshold (see Sec.~\ref{sec:comp-details}). An inner (micro-iteration) loop (blue) is introduced to accelerate the convergence of the self-consistent calculation, in which the set of determinants in $\Psi^{\mu\,(k,l)}$ is kept fixed, and only the diagonalization of $\hat{H}^{\mu\,(k,l)}$ is performed iteratively with the updated density $n^{(k,l)}$. The inner loop is exited when the absolute energy difference between two successive micro-iterations $\Delta E^{(k,l)}$ is below a threshold $\tau_2$ that has been here set to $10^{-2} \times \tau_1$ hartree. The convergence of the algorithm was further improved by introducing a direct inversion in the iterative subspace (DIIS) step to extrapolate the density both in the outer and inner loops. \cite{Pulay_1980,Pulay_1982} Note that any range-separated post-HF method can be implemented using this scheme by just replacing the CIPSI step by the post-HF method of interest. Note that, thanks to the self-consistent nature of the algorithm, the final trial wave function $\Psi^{\mu}$ is independent of the starting wave function $\Psi^{(0)}$. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Computational details} \label{sec:comp-details} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% For all the systems considered here, experimental geometries have been used and they have been extracted from the NIST website \titou{[REF]}. Geometries for each system are reported in the {\SI}. All the calculations have been performed using Burkatzki-Filippi-Dolg (BFD) pseudopotentials \cite{Burkatzki_2007,Burkatzki_2008} with the associated double-, triple-, and quadruple-$\zeta$ basis sets (BFD-VXZ). The small-core BFD pseudopotentials include scalar relativistic effects. CCSD(T) and KS-DFT energies have been computed with \emph{Gaussian09},\cite{g16} using the unrestricted formalism for open-shell systems. All the CIPSI calculations have been performed with \emph{Quantum Package}.\cite{Garniron_2019,qp2_2020} We used the short-range version of \titou{the local-density approximation (LDA)} and Perdew-Burke-Ernzerhof (PBE) \cite{PerBurErn-PRL-96} exchange and correlation functionals defined in Ref.~\onlinecite{GolWerStoLeiGorSav-CP-06} (see also Refs.~\onlinecite{TouColSav-JCP-05,GolWerSto-PCCP-05}). The convergence criterion for stopping the CIPSI calculations has been set to $\EPT < 10^{-3}$ hartree or $ \Ndet > 10^7$. All the wave functions are eigenfunctions of the $\Hat{S}^2$ spin operator, as described in Ref.~\onlinecite{Applencourt_2018}. Quantum Monte Carlo calculations were made with QMC=Chem,\cite{scemama_2013} in the determinant localization approximation (DLA),\cite{Zen_2019} where only the determinantal component of the trial wave function is present in the expression of the wave function on which the pseudopotential is localized. Hence, in the DLA the fixed-node energy is independent of the Jastrow factor, as in all-electron calculations. Simple Jastrow factors were used to reduce the fluctuations of the local energy. The FN-DMC simulations are performed with the stochastic reconfiguration algorithm developed by Assaraf \textit{et al.}, \cite{Assaraf_2000} with a time step of $5 \times 10^{-4}$ a.u. \titou{All-electron move DMC?} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Influence of the range-separation parameter on the fixed-node error} \label{sec:mu-dmc} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% TABLE I %%% \begin{table} \caption{Fixed-node energies $\EDMC$ (in hartree) and number of determinants $\Ndet$ in \ce{H2O} with various trial wave functions $\Psi^{\mu}$.} \label{tab:h2o-dmc} \centering \begin{ruledtabular} \begin{tabular}{crlrl} & \multicolumn{2}{c}{BFD-VDZ} & \multicolumn{2}{c}{BFD-VTZ} \\ \cline{2-3} \cline{4-5} $\mu$ & $\Ndet$ & $\EDMC$ & $\Ndet$ & $\EDMC$ \\ \hline $0.00$ & $11$ & $-17.253\,59(6)$ & $23$ & $-17.256\,74(7)$ \\ $0.20$ & $23$ & $-17.253\,73(7)$ & $23$ & $-17.256\,73(8)$ \\ $0.30$ & $53$ & $-17.253\,4(2)$ & $219$ & $-17.253\,7(5)$ \\ $0.50$ & $1\,442$ & $-17.253\,9(2)$ & $16\,99$ & $-17.257\,7(2)$ \\ $0.75$ & $3\,213$ & $-17.255\,1(2)$ & $13\,362$ & $-17.258\,4(3)$ \\ $1.00$ & $6\,743$ & $-17.256\,6(2)$ & $256\,73$ & $-17.261\,0(2)$ \\ $1.75$ & $54\,540$ & $-17.259\,5(3)$ & $207\,475$ & $-17.263\,5(2)$ \\ $2.50$ & $51\,691$ & $-17.259\,4(3)$ & $858\,123$ & $-17.264\,3(3)$ \\ $3.80$ & $103\,059$ & $-17.258\,7(3)$ & $1\,621\,513$ & $-17.263\,7(3)$ \\ $5.70$ & $102\,599$ & $-17.257\,7(3)$ & $1\,629\,655$ & $-17.263\,2(3)$ \\ $8.50$ & $101\,803$ & $-17.257\,3(3)$ & $1\,643\,301$ & $-17.263\,3(4)$ \\ $\infty$ & $200\,521$ & $-17.256\,8(6)$ & $1\,631\,982$ & $-17.263\,9(3)$ \\ \end{tabular} \end{ruledtabular} \end{table} %%% %%% %%% %%% %%% FIG 2 %%% \begin{figure} \centering \includegraphics[width=\columnwidth]{h2o-dmc.pdf} \caption{Fixed-node energies of \ce{H2O} for different values of $\mu$, using the srLDA or srPBE density functionals to build the trial wave function.} \label{fig:h2o-dmc} \end{figure} %%% %%% %%% %%% The first question we would like to address is the quality of the nodes of the wave function $\Psi^{\mu}$ obtained with an intermediate range separation parameter (\ie, $0 < \mu < +\infty$). For this purpose, we consider a weakly correlated molecular system, namely the water molecule \titou{near its equilibrium geometry.} \cite{Caffarel_2016} We then generate trial wave functions $\Psi^\mu$ for multiple values of $\mu$, and compute the associated fixed-node energy keeping all the parameters having an impact on the nodal surface fixed (\titou{such as ??}). %====================================================== \subsection{Fixed-node energy of $\Psi^\mu$} \label{sec:fndmc_mu} %====================================================== From Table~\ref{tab:h2o-dmc} and Fig.~\ref{fig:h2o-dmc}, one can clearly observe that relying on FCI trial wave functions ($\mu \to \infty$) give FN-DMC energies lower than the energies obtained with a single KS determinant ($\mu=0$): a gain of $3.2 \pm 0.6$~m\hartree{} at the double-$\zeta$ level and $7.2 \pm 0.3$~m\hartree{} at the triple-zeta level are obtained. Coming now to the nodes of the trial wave functions $\Psi^{\mu}$ with intermediate values of $\mu$, Fig.~\ref{fig:h2o-dmc} shows that a smooth behaviour is obtained: starting from $\mu=0$ (\textit{i.e.} the KS determinant), the FN-DMC error is reduced continuously until it reaches a minimum for an optimal value of $\mu$, and then the FN-DMC error raises until it reaches the $\mu=\infty$ limit (\textit{i.e.} the FCI wave function). For instance, with respect to the FN-DMC energy of the FCI trial wave function in the double zeta basis set, with the optimal value of $\mu$, one can obtain a lowering of the FN-DMC energy of $2.6 \pm 0.7$~m\hartree{}, with an optimal value at $\mu=1.75$~bohr$^{-1}$. When the basis set is increased, the gain in FN-DMC energy with respect to the FCI trial wave function is reduced, and the optimal value of $\mu$ is slightly shifted towards large $\mu$. Eventually, the nodes of the wave functions $\Psi^\mu$ obtained with the short-range LDA exchange-correlation functional give very similar FN-DMC energies with respect to those obtained with the short-range PBE functional, even if the RS-DFT energies obtained with these two functionals differ by several tens of m\hartree{}. %====================================================== \subsection{Link between RS-DFT and Jastrow factors } \label{sec:rsdft-j} %====================================================== The data obtained in \ref{sec:fndmc_mu} show that RS-DFT can provide CI coefficients giving trial wave functions with better nodes than FCI wave functions. Such behaviour can be compared to the common practice of re-optimizing the Slater part of a trial wave function in the presence of the Jastrow factor. In the present paragraph, we would like therefore to elaborate some more on the link between RS-DFT and wave function optimization within the presence of a Jastrow factor. Let us assume a fixed Jastrow factor $J({\bf{r}_1}, \hdots , {\bf{r}_N})$, and a corresponding Slater-Jastrow wave function $\Phi = e^J \Psi$, where $\Psi = \sum_I c_I D_I$ is a general linear combination of Slater determinants $D_I$. The only remaining variational parameters in $\Phi$ are therefore the Slater part $\Psi$. Let us call $\Psi^J$ the linear combination of Slater determinant minimizing the variational energy \begin{equation} \Psi^J = \text{argmin}_{\Psi}\frac{ \mel{ \Psi }{ e^{J} H e^{J} }{ \Psi } }{\mel{ \Psi }{ e^{2J} }{ \Psi } }. \end{equation} Such a wave function $\Psi^J$ satisfies the generalized hermitian eigenvalue equation \begin{equation} e^{J} H e^{J} \Psi^J = E e^{2J} \Psi^J, \end{equation} but also the non-hermitian transcorrelated eigenvalue problem\cite{many_things} \begin{equation} \label{eq:transcor} e^{-J} H e^{J} \Psi^J = E \Psi^J, \end{equation} which is much easier to handle despite its non-hermicity. Of course, the FN-DMC energy of $\Phi$ depends therefore only on the nodes of $\Psi^J$. In a finite basis set and with a quite accurate Jastrow factor, it is known that the nodes of $\Psi^J$ may be better than that of the FCI wave function, and therefore, we would like to compare $\Psi^J$ and $\Psi^\mu$. To do so, we have made the following numerical experiment. First, we extract the 200 determinants with the largest weights in the FCI wave function out of a large CIPSI calculation. Within this set of determinants, we solve the self-consistent equations of RS-DFT [see Eq.~\eqref{rs-dft-eigen-equation}] with different values of $\mu$. This gives the CI expansions $\Psi^\mu$. Then, within the same set of determinants we optimize the CI coefficients in the presence of a simple one- and two-body Jastrow factor. This gives the CI expansion $\Psi^J$. Then, we can easily compare $\Psi^\mu$ and $\Psi^J$ as they are developed on the same Slater determinant basis. In figure~\ref{fig:overlap}, we plot the overlaps $\braket{\Psi^J}{\Psi^\mu}$ obtained for water and the fluorine dimer, and in figure~\ref{dmc_small} the FN-DMC energy of the wave functions $\Psi^\mu$ together with that of $\Psi^J$. %%% FIG 3 %%% \begin{figure} \centering \includegraphics[width=\columnwidth]{overlap.pdf} \caption{\ce{H2O}, double-zeta basis set, 200 most important determinants of the FCI expansion (see Sec.~\ref{sec:rsdft-j}). Overlap of the RS-DFT CI expansions $\Psi^\mu$ with the CI expansion optimized in the presence of a Jastrow factor $\Psi^J$.} \label{fig:overlap} \end{figure} %%% %%% %%% %%% %%% FIG 4 %%% \begin{figure} \centering \includegraphics[width=\columnwidth]{h2o-200-dmc.pdf} \caption{\ce{H2O}, double-zeta basis set, 200 most important determinants of the FCI expansion (see \ref{sec:rsdft-j}). FN-DMC energies of $\Psi^\mu$ (red curve), together with the FN-DMC energy of $\Psi^J$ (blue line). The width of the lines represent the statistical error bars.} \label{dmc_small} \end{figure} %%% %%% %%% %%% There is a clear maximum of overlap at $\mu=1$~bohr$^{-1}$, which coincides with the minimum of the FN-DMC energy of $\Psi^\mu$. Also, it is interesting to notice that the FN-DMC energy of $\Psi^J$ is compatible with that of $\Psi^\mu$ with $0.5 < \mu < 1$~bohr$^{-1}$. This confirms that introducing short-range correlation with DFT has an impact on the CI coefficients similar to the Jastrow factor. %%% TABLE II %%% \begin{table} \caption{\ce{H2O}, double-zeta basis set. Integrated on-top pair density $\expval{ n_2({\bf r},{\bf r}) }$ for $\Psi^J$ and $\Psi^\mu$ with different values of $\mu$. } \label{table_on_top} \begin{ruledtabular} \begin{tabular}{cc} $\mu$ & $\expval{ n_2({\bf r},{\bf r}) }$ \\ \hline 0.00 & 1.443 \\ 0.25 & 1.438 \\ 0.50 & 1.423 \\ 1.00 & 1.378 \\ 2.00 & 1.325 \\ 5.00 & 1.288 \\ $\infty$ & 1.277 \\ \hline $\Psi^J$ & 1.404 \\ \end{tabular} \end{ruledtabular} \end{table} %%% %%% %%% %%% %%% FIG 5 %%% \begin{figure} \centering \includegraphics[width=\columnwidth]{on-top-mu.pdf} \caption{\ce{H2O}, double-zeta basis set. On-top pair density $n_2({\bf r},{\bf r})$ along the O---H axis, for $\Psi^J$ (dashed curve) and $\Psi^\mu$ with different values of $\mu$. } \label{fig:n2} \end{figure} %%% %%% %%% %%% %%% FIG 6 %%% \begin{figure} \centering \includegraphics[width=\columnwidth]{density-mu.pdf} \caption{\ce{H2O}, double-zeta basis set. Density $n({\bf r})$ along the O---H axis, for $\Psi^J$ (dashed curve) and $\Psi^\mu$ with different values of $\mu$. } \label{fig:n1} \end{figure} %%% %%% %%% %%% In order to refine the comparison between $\Psi^\mu$ and $\Psi^J$, we report several quantities related to the one- and two-body density of $\Psi^J$ and $\Psi^\mu$ with different values of $\mu$. First, we report in table~\ref{table_on_top} the integrated on-top pair density $\expval{ n_2({\bf r},{\bf r}) }$ \begin{equation} \expval{ n_2({\bf r},{\bf r}) } = \int \text{d}{\bf r} \,\,n_2({\bf r},{\bf r}) \end{equation} where $n_2({\bf r}_1,{\bf r}_2)$ is the two-body density (normalized to $N(N-1)$ where $N$ is the number of electrons) obtained for both $\Psi^\mu$ and $\Psi^J$. Then, in order to have a pictorial representation of both the on-top pair density and the density, we report in Fig.~\ref{fig:n1} and Fig.~\ref{fig:n2} the plots of the total density $n({\bf r})$ and on-top pair density $n_2({\bf r},{\bf r})$ along one O---H axis of the water molecule. From these data, one can clearly notice several trends. First, from Tab.~\ref{table_on_top}, we can observe that the overall on-top pair density decreases when $\mu$ increases. This is expected as the two-electron interaction increases in $H^\mu[n]$. Second, the relative variations of the on-top pair density with $\mu$ are much more important than that of the one-body density, the latter being essentially unchanged between $\mu=0$ and $\mu=\infty$ while the former can vary by about 10$\%$ in some regions. %TODO TOTO In the high-density region of the O---H bond, the value of the on-top pair density obtained from $\Psi^J$ is superimposed with $\Psi^{\mu=0.5}$, and at a large distance the on-top pair density is the closest to $\mu=\infty$. The integrated on-top pair density obtained with $\Psi^J$ lies between the values obtained with $\mu=0.5$ and $\mu=1$~bohr$^{-1}$, constently with the FN-DMC energies and the overlap curve. These data suggest that the wave functions $\Psi^{0.5 \le \mu \le 1}$ and $\Psi^J$ are close, and therefore that the operators that produced these wave functions (\textit{i.e.} $H^\mu[n]$ and $e^{-J}He^J$) contain similar physics. Considering the form of $\hat{H}^\mu[n]$ (see Eq.~\eqref{H_mu}), one can notice that the differences with respect to the usual Hamiltonian come from the non-divergent two-body interaction $\hat{W}_{\text{ee}}^{\text{lr},\mu}$ and the effective one-body potential $\hat{\bar{V}}_{\text{Hxc}}^{\text{sr},\mu}[n]$ which is the functional derivative of the Hartree-exchange-correlation functional. The roles of these two terms are therefore very different: with respect to the exact ground state wave function $\Psi$, the non divergent two body interaction increases the probability to find electrons at short distances in $\Psi^\mu$, while the effective one-body potential $\hat{\bar{V}}_{\text{Hxc}}^{\text{sr},\mu}[n_{\Psi^{\mu}}]$, provided that it is exact, maintains the exact one-body density. This is clearly what has been observed from the plots in Fig.~\ref{fig:n1} and Fig.~\ref{fig:n2}. Regarding now the transcorrelated Hamiltonian $e^{-J}He^J$, as pointed out by Ten-No,\cite{Ten-no2000Nov} the effective two-body interaction induced by the presence of a Jastrow factor can be non-divergent when a proper Jastrow factor is chosen. Therefore, one can understand the similarity between the eigenfunctions of $H^\mu$ and the Jastrow-Slater optimization: they both deal with an effective non-divergent interaction but still produce a reasonable one-body density. As a conclusion of the first part of this study, we can notice that: \begin{itemize} \item with respect to the nodes of a KS determinant or a FCI wave function, one can obtain a multi determinant trial wave function $\Psi^\mu$ with a smaller fixed node error by properly choosing an optimal value of $\mu$ in RS-DFT calculations, \item the optimal value of $\mu$ depends on the system and the basis set, and the larger the basis set, the larger the optimal value of $\mu$, \item numerical experiments (overlap $\braket{\Psi^\mu}{\Psi^J}$, one-body density, on-top pair density, and FN-DMC energy) indicate that the RS-DFT scheme essentially plays the role of a simple Jastrow factor, \textit{i.e.} mimicking short-range correlation effects. The latter statement can be qualitatively understood by noticing that both RS-DFT and transcorrelated approaches deal with an effective non-divergent electron-electron interaction, while keeping the density constant. \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Energy differences in FN-DMC: atomization energies} \label{sec:atomization} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Atomization energies are challenging for post-Hartree-Fock methods because their calculation requires a perfect balance in the description of atoms and molecules. Basis sets used in molecular calculations are atom-centered, so they are always better adapted to atoms than molecules and atomization energies usually tend to be underestimated with variational methods. In the context of FN-DMC calculations, the nodal surface is imposed by the trial wavefunction which is expanded on an atom-centered basis set, so we expect the fixed-node error to be also tightly related to the basis set incompleteness error. Increasing the size of the basis set improves the description of the density and of electron correlation, but also reduces the imbalance in the quality of the description of the atoms and the molecule, leading to more accurate atomization energies. %============================ \subsection{Size consistency} %============================ An extremely important feature required to get accurate atomization energies is size-consistency (or strict separability), since the numbers of correlated electron pairs in the isolated atoms are different from those of the molecules. The energy computed within density functional theory is size-consistent, and as it is a mean-field method the convergence to the complete basis set (CBS) limit is relatively fast. Hence, DFT methods are very well adapted to the calculation of atomization energies, especially with small basis sets. But going to the CBS limit will converge to biased atomization energies because of the use of approximate density functionals. On the other hand, FCI is also size-consistent, but the convergence of the FCI energies to the CBS limit is much slower because of the description of short-range electron correlation using atom-centered functions. But ultimately the exact energy will be reached. In the context of selected CI calculations, when the variational energy is extrapolated to the FCI energy\cite{Holmes_2017} there is no size-consistency error. But when the truncated SCI wave function is used as a reference for post-Hartree-Fock methods such as SCI+PT2 or for QMC calculations, there is a residual size-consistency error originating from the truncation of the wave function. QMC energies can be made size-consistent by extrapolating the FN-DMC energy to estimate the energy obtained with the FCI as a trial wave function.\cite{Scemama_2018,Scemama_2018b} Alternatively, the size-consistency error can be reduced by choosing the number of selected determinants such that the sum of the PT2 corrections on the fragments is equal to the PT2 correction of the molecule, enforcing that the variational potential energy surface (PES) is parallel to the perturbatively corrected PES, which is a relatively accurate estimate of the FCI PES.\cite{Giner_2015} Another source of size-consistency error in QMC calculations originates from the Jastrow factor. Usually, the Jastrow factor contains one-electron, two-electron and one-nucleus-two-electron terms. The problematic part is the two-electron term, whose simplest form can be expressed as \begin{equation} J_\text{ee} = \sum_i \sum_{j