\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{\alert}[1]{\textcolor{red}{#1}} \definecolor{darkgreen}{HTML}{009900} \usepackage[normalem]{ulem} \newcommand{\LCT}{Laboratoire de Chimie Théorique (UMR 7616), Sorbonne Université, 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é de Toulouse, CNRS, UPS, France} \begin{document} \title{Enabling high accuracy diffusion Monte Carlo calculations with range-separated density functional theory and selected configuration interaction} \author{Anthony Scemama} \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} \end{abstract} \maketitle \section{Introduction} \label{sec:intro} The full configuration interaction (FCI) method leads to the exact solution of the Schrödinger equation with an approximate Hamiltonian expressed in a finite basis of Slater determinants. The FCI wave function can be interpreted as the exact solution of the true Hamiltonian obtained with the additional constraint that it can only span the space provided by the basis. At the complete basis set (CBS) limit, the constraint vanishes and the exact solution is obtained. Hence the FCI method enables a systematic improvement of the calculations by increasing the size of the basis set. Nevertheless, its exponential scaling with the number of electrons and with the size of the basis is prohibitive to treat large 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} 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, but the scaling remains exponential unless some bias is introduced leading to a loss of size extensivity. The Diffusion Monte Carlo (DMC) method is a numerical scheme to obtain the exact solution of the Schrödinger equation with an additional constraint, imposing the solution to have the same nodal hypersurface as a given trial wave function. This approximation is known as the \emph{fixed-node} approximation. When the nodes of the trial wave function coincide with the nodes of the exact wave function, the exact energy and wave function are obtained. The DMC method is attractive because its scaling is polynomial with the number of electrons and with the size of the trial wave function. Moreover, the total energies obtained are usually below those obtained with FCI because the fixed-node approximation imposes less constraints on the solution than the finite-basis approximation. In many cases, the systems under study are well described by a single Slater determinant. Single-determinant DMC can be used as a post-Hatree-Fock single-reference method with an accuracy comparable to coupled cluster. The favorable scaling of QMC and its adequation with massively parallel architectures makes it an attractive alternative for large systems. It has been shown that the nodal surfaces obtained with Kohn-Sham determinants are in general better than those obtained with the Hartree-Fock determinant,\cite{Per_2012} and of comparable quality to those obtained with natural orbitals of single-determinant correlated calculations.\cite{Wang_2019} However, the fixed-node approximation is much more difficult to control than the finite-basis approximation, as it is not possible to minimize directly the DMC energy with respect to the variational parameters of the trial wave function. 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 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 DMC.\cite{Petruzielo_2012} Another approach consists in considering the DMC method as a \emph{post-FCI method}. The trial wave function is obtained by approaching the FCI with a selected configuration interaction method such as CIPSI for instance.\cite{Giner_2013,Caffarel_2016_2} When the basis set is increased, the trial wave function tends to the exact wave function, so the nodal surface can be systematically improved.\cite{Caffarel_2016} This technique has the advantage that using FCI nodes in a given basis set is well defined and has a unique solution. The optimization of the wave function is deterministic, so the calculations are reproducible and don't require the expertise of a QMC expert. However, this technique can't be applied to large systems because of the exponential scaling of the size of the wave function. Extrapolation techniques have been applied to estimate the DMC energy of a FCI wave function in a large basis set,\cite{Scemama_2018} and other authors have used a combination of the two approaches where CIPSI trial wave functions are re-optimized under the presence of a Jastrow factor.\cite{Giner_2016,Dash_2018,Dash_2019} \section{Combining range-separated DFT with CIPSI} \label{sec:rsdft-cipsi} Starting from a Hartree-Fock determinant in a small basis set, we have seen that we can systematically improve the trial wave function in two directions. The first one is by increasing the size of the atomic basis set, and the second one is by increasing the determinant expansion towards the FCI limit. A third direction of improvement exists. It is the path which connects the Hartree-Fock determinant to the Kohn-Sham determinant, which usually has a better nodal surface. This third path is obtained by the so-called \emph{range-separated} DFT framework which splits the electron-electron interaction in a long-range part and a short-range part. HELP MANU! \subsection{CIPSI} \emph{Configuration interaction using a perturbative selection made iteratively} (CIPSI)\cite{Huron_1973} belongs to the class of selected CI methods, whose goal is to build compressed configuration interaction (CI) wave functions with a controlled accuracy, together with an estimate of the associated eigenvalue obtained with perturbation theory. Using a conventional iterative diagonalization method such as Davidson's,\cite{Davidson_1975} the size of the CI space is pre-defined, and at each iteration the length of the vectors is equal to the size of the CI space. Within selected CI methods, at each iteration the size of the determinant space is increased such that one keeps only the most important determinants that will lead to an accurate description of the state of interest. The lowest eigenpairs are extracted from the CI matrix expressed in this determinant subspace, and the CI energy can be estimated by computing a second-order perturbative correction (PT2) to the variational energy computed in the complete CI space. The magnitude of the PT2 correction is a measure of the distance to the exact eigenvalue, and is an adjustable parameter to control the quality of the wave function. \subsection{Range-separated DFT in a nutshell} \label{sec:rsdft} \subsubsection{Exact equations} \label{sec:exact} The exact ground-state energy of a $N$-electron system with nuclei-electron potential $v_\mathrm{ne}(\textbf{r})$ can be expressed by the following minimization over $N$-representable densities $n$~\cite{Lev-PNAS-79,Lie-IJQC-83} \begin{equation} E_0 = \min_n \left\{ \mathcal{F}[n] + \int v_\mathrm{ne}(\textbf{r}) n(\textbf{r}) \mathrm{d} \textbf{r} \right\}, \label{E0} \end{equation} with the standard constrained-search universal density functional \begin{equation} \mathcal{F}[n] = \min_{\Psi\rightarrow n} \left \langle\Psi|\hat{T}+\hat{W}_\mathrm{ee}|\Psi \right \rangle, \label{Fn} \end{equation} where $\hat{T}$ and $\hat{W}_\mathrm{ee}$ are the kinetic-energy and Coulomb electron-electron interaction operators, respectively. The minimizing multideterminant wave function in Eq.~\eqref{Fn} will be denoted by $\Psi[n]$. In RS-DFT, the universal density functional is decomposed as~\cite{Sav-INC-96,TouColSav-PRA-04} \begin{equation} \mathcal{F}[n] = \mathcal{F}^{\mathrm{lr},\mu}[n] + \bar{E}_{\mathrm{Hxc}}^{\mathrm{sr,}\mu}[n], \label{Fdecomp} \end{equation} where $\mathcal{F}^{\mathrm{lr},\mu}[n]$ is a long-range (lr) universal density functional \begin{equation} \label{lr_univ_fonc} \mathcal{F}^{\mathrm{lr},\mu}[n]= \min_{\Psi\rightarrow n} \left \langle\Psi|\hat{T}+\hat{W}_\mathrm{ee}^{\mathrm{lr},\mu}|\Psi \right \rangle, \end{equation} and $\bar{E}_{\mathrm{Hxc}}^{\,\mathrm{sr,}\mu}[n]$ is the complementary short-range (sr) Hartree-exchange-correlation (Hxc) density functional. In Eq.~(\ref{lr_univ_fonc}), $\hat{W}_\mathrm{ee}^{\mathrm{lr}}$ is the long-range electron-electron interaction defined as \begin{equation} \hat{W}_\mathrm{ee}^{\mathrm{lr},\mu} = \frac{1}{2} \iint w_{\mathrm{ee}}^{\mathrm{lr},\mu}(r_{12}) \hat{n}_2(\textbf{r}_1,\textbf{r}_2) \mathrm{d} \textbf{r}_1 \mathrm{d} \textbf{r}_2, \end{equation} with the error-function potential $w_{\mathrm{ee}}^{\mathrm{lr},\mu}(r_{12})=\mathrm{erf}(\mu\, r_{12} )/r_{12}$ (expressed with the interelectronic distance $r_{12} = ||\textbf{r}_1-\textbf{r}_2||$) and the pair-density operator $\hat{n}_2(\textbf{r}_1, \textbf{r}_2)=\hat{n}(\textbf{r}_1) \hat{n}(\textbf{r}_2) - \delta(\textbf{r}_1-\textbf{r}_2) \hat{n}(\textbf{r}_1)$ where $\hat{n}(\textbf{r})$ is the density operator. The range-separation parameter $\mu$ corresponds to an inverse distance controlling the range of the separation and the strength of the interaction at $r_{12} = 0$. For a given density, we will denote by $\Psi^\mu[n]$ the minimizing multideterminant wave function in Eq.~\eqref{lr_univ_fonc}. Inserting the decomposition of Eq.~\eqref{Fdecomp} into Eq.~\eqref{E0}, and recomposing the two-step minimization into a single one, leads to the following expression for the exact ground-state electronic energy \begin{equation} \label{min_rsdft} E_0= \min_{\Psi} \left\{ \left \langle\Psi|\hat{T}+\hat{W}_\mathrm{{ee}}^{\mathrm{lr},\mu}+\hat{V}_{\mathrm{ne}}|\Psi\right \rangle + \bar{E}^{\mathrm{sr},\mu}_{\mathrm{Hxc}}[n_\Psi]\right\}, \end{equation} where the minimization is done over normalized $N$-electron multideterminant wave functions, $\hat{V}_{\mathrm{ne}} = \int v_{\mathrm{ne}} (\textbf{r}) \hat{n}(\textbf{r}) \mathrm{d} \textbf{r}$, and $n_\Psi$ refers to the density of $\Psi$, i.e. $n_\Psi(\textbf{r})=\left \langle\Psi|\hat{n}(\textbf{r})|\Psi\right \rangle$. The minimizing multideterminant wave function $\Psi^\mu$ in Eq.~\eqref{min_rsdft} 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}_{\mathrm{ee}}^{\mathrm{lr},\mu}+\hat{V}_{\mathrm{ne}}+ \hat{\bar{V}}_{\mathrm{Hxc}}^{\mathrm{sr},\mu}[n_{\Psi^{\mu}}], \end{equation} where $\hat{\bar{V}}_{\mathrm{Hxc}}^{\mathrm{sr},\mu}[n]=\int \delta \bar{E}^{\mathrm{sr},\mu}_{\mathrm{Hxc}}[n]/\delta n(\textbf{r}) \, \hat{n}(\textbf{r}) \mathrm{d} \textbf{r}$ is the complementary short-range Hartree-exchange-correlation potential operator. Note that $\Psi^{\mu}$ is not the exact physical ground-state wave function but only an effective wave function. However, its density $n_{\Psi^{\mu}}$ is the exact physical ground-state density. Once $\Psi^{\mu}$ has been calculated, the exact electronic ground-state energy is obtained by \begin{equation} \label{E-rsdft} E_0= \braket{\Psi^{\mu}|\hat{T}+\hat{W}_\mathrm{{ee}}^{\mathrm{lr},\mu}+\hat{V}_{\mathrm{ne}}|\Psi^{\mu}}+\bar{E}^{\mathrm{sr},\mu}_{\mathrm{Hxc}}[n_{\Psi^\mu}]. \end{equation} Note that, for $\mu=0$, the long-range interaction vanishes, $w_{\mathrm{ee}}^{\mathrm{lr},\mu=0}(r_{12}) = 0$, and thus RS-DFT reduces to standard KS-DFT. For $\mu\to\infty$, the long-range interaction becomes the standard Coulomb interaction, $w_{\mathrm{ee}}^{\mathrm{lr},\mu\to\infty}(r_{12}) = 1/r_{12}$, and thus RS-DFT reduces to standard wave-function theory (WFT). \subsection{Approximations} \label{sec:approx} Provided that the exact short-range complementary functional is known and that the wave function $\Psi^{\mu}$ is developed in a complete basis set, the theory is exact and therefore provides the exact density and ground state energy whatever value of $\mu$ is chosen to split the universal Levy-Lieb functional (see Eqs. \eqref{Fdecomp} and \eqref{lr_univ_fonc}). Of course in practice, two kind of approximations must be made: i) the approximation of the wave function by using a finite basis set, ii) the approximation on the exact complementary functionals. As long as approximations are made, the theory is not exact anymore and might depend on $\mu$. Despite its $\mu$ dependence, approximated RS-DFT schemes provide potentially appealing features as: i) the approximated wave function $\Psi^{\mu}$ is an eigenvector of a non-divergent operator $\hat{H}^\mu[n_{\Psi^{\mu}}]$ (see Eqs. \eqref{H_mu} and \eqref{rs-dft-eigen-equation}) and therefore converges more rapidly with respect to the basis set~\cite{FraMusLupTou-JCP-15}, and also produces more compact wave function~\cite{GinPraFerAssSavTou-JCP-18} as it will be illustrated here, ii) the semi-local approximations for RS-DFT complementary functionals are usually better suited to describe the remaining short-range part of the electron-electron correlation. In this work, we use the short-range version of the Perdew-Burke-Ernzerhof (PBE)~\cite{PerBurErn-PRL-96} exchange and correlation functionals of Ref.~\onlinecite{GolWerStoLeiGorSav-CP-06} (see also Refs.~\onlinecite{TouColSav-JCP-05,GolWerSto-PCCP-05}) which takes the form \begin{eqnarray} \bar{E}^{\mathrm{sr},\mu,\textsc{pbe}}_{\mathrm{x/c}}[n] = \int \bar{e}_\mathrm{{x/c}}^\mathrm{sr,\mu,\textsc{pbe}}(n(\textbf{r}),\nabla n(\textbf{r})) \, \mathrm{d}\textbf{r}. \end{eqnarray} \subsection{RSDFT-CIPSI} It is possible to use DFT for short-range interactions and CIPSI for the long-range. This scheme has been recently implemented.\cite{GinPraFerAssSavTou-JCP-18} \begin{figure}[h] \centering \includegraphics[width=\columnwidth]{algorithm.pdf} \caption{Algorithm showing the generation of the RSDFT-CIPSI wave function} \end{figure} \section{Influence of the range-separation parameter on the fixed-node error} \label{sec:mu-dmc} \begin{table}[h] \caption{Fixed-node energies of the water molecule.} \label{tab:h2o-dmc} \centering \begin{tabular}{crlrl} \hline & \multicolumn{2}{c}{BFD-VDZ} & \multicolumn{2}{c}{BFD-VTZ}\\ $\mu$ & $N_{\text{det}}$ & E(DMC) & $N_{\text{det}}$ & E(DMC)\\ \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)$\\ \hline \end{tabular} \end{table} \begin{figure}[h] \centering \includegraphics[width=\columnwidth]{h2o-dmc.pdf} \caption{Fixed-node energies of the water molecule with variable values of $\mu$.} \label{fig:h2o-dmc} \end{figure} The water molecule was taken at the equilibrium geometry,\cite{Caffarel_2016} and RSDFT-CIPSI wave functions were generated with BFD pseudopotentials and the corresponding double-zeta basis set using multiple values of the range-separation parameter $\mu$. The convergence criterion for stopping the CIPSI calculation was set to 1~m$E_h$ on the PT2 correction. Then, these wave functions were used as trial wave functions for FN-DMC calculations, and the corresponding energies are shown in table~\ref{tab:h2o-dmc} and figure~\ref{fig:h2o-dmc}. Using FCI trial wave functions gives FN-DMC energies which are lower than the energies obtained with a single Kohn-Sham determinant: 3~m$E_h$ at the double-zeta level and 7~m$E_h$ at the triple-zeta level. Interestingly, with the double-zeta basis one can obtain a FN-DMC energy 2.5~m$E_h$ lower than the energy obtained with the FCI trial wave function, using the RSDFT-CIPSI with a range-separation parameter $\mu=1.75$. This can be explained by the inability of the basis set to properly describe short-range correlation, shifting the nodes from their optimal position. Using DFT to take account of short-range correlation frees the determinant expansion from describing short-range effects, and enables a better placement of the nodes. At the triple-zeta level, the short-range correlations can be better described, and the improvement due to DFT is insignificant. However, it is important to note that the same FN-DMC energy can be obtained with a CI expansion which is eight times smaller when sr-DFT is introduced. One can also remark that the minimum has been shifted towards the FCI, which is consistent with the fact that in the CBS limit we expect the minimum of the FN-DMC energy to be obtained for the FCI wave function, at $\mu=\infty$. \section{Computational details} \label{sec:comp-details} All the calculations were made using BFD pseudopotentials\cite{Burkatzki_2008} with the associated double, triple and quadruple zeta basis sets. CCSD(T) and DFT calculations were made with \emph{Gaussian09},\cite{g16} using an unrestricted Hartree-Fock determinant as a reference for open-shell systems. All the CIPSI calculations and range-separated CIPSI calculations were made with \emph{Quantum Package}.\cite{Garniron_2019,qp2_2020} Quantum Monte Carlo calculations were made with QMC=Chem,\cite{scemama_2013} in the determinant localization approximation.\cite{Zen_2019} In the determinant localization approximation, 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, the pseudopotential operator does not depend on the Jastrow factor, as it is the case in all-electron calculations. This improves the reproducibility of the results, as they depend only on parameters optimized in a deterministic framework. \section{Atomization energy benchmarks} \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. 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 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. Another important feature required to get accurate atomization energies is size-extensivity, since the numbers of correlated electrons in the isolated atoms are different from the number of correlated electrons in the molecule. In the context of selected CI calculations, when the variational energy is extrapolated to the FCI energy\cite{Holmes_2017} there is obviously no size-consistence error. But when the selected wave function is used for as a reference for post-Hartree-Fock methods or QMC calculations, there is a residual size-consistence error originating from the truncation of the determinant space. QMC calculations 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-consistence error can be reduced by choosing the number of selected determinants such that the sum of the PT2 corrections on the atoms is equal to the PT2 correction of the molecule, enforcing that the variational dissociation potential energy surface (PES) is parallel to the perturbatively corrected PES, which is an accurate estimate of the FCI PES.\cite{Giner_2015} Another source of size-consistence error in QMC calculation may originate 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