RSDFT-CIPSI-QMC/Manuscript/rsdft-cipsi-qmc.tex

\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}}
\newcommand{\eg}[1]{\textcolor{blue}{#1}}
\definecolor{darkgreen}{HTML}{009900}
\usepackage[normalem]{ulem}

\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}


\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}
\end{abstract}

\maketitle


\section{Introduction}
\label{sec:intro}

Within a finite one-electron basis, full configuration interaction (FCI) 
delivers only an approximate solution of the Schr\"odinger equation.
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{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 CIPSI 3-CLASS AND INITIATOR APPROXIMATION.

Diffusion Monte Carlo (DMC) is a 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.

\alert{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 there is no variational principle for the nodes.}

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.
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 on 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}
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}
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 can't 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 range-separated 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 $|\EPT| <
1$~mE$_h$.



\subsection{Range-separated DFT}
\label{sec:rsdft}

Following the seminal work of Savin,\cite{Savin_1996,Toulouse_2004}
the Coulomb electron-electron interaction is split into a short-range
(sr) and a long range (lr) interaction as
\begin{equation}
  \frac{1}{r_{ij}} = w_{\text{ee}}^{\text{lr}, \mu}(r_{ij}) + \qty(
  \frac{1}{r_{ij}} - w_{\text{ee}}^{\text{lr}, \mu}(r_{ij}) )
\end{equation}
where
\begin{equation}
  w_{\text{ee}}^{\text{lr},\mu}(r_{ij}) = \frac{\erf \qty( \mu\, r_{ij})}{r_{ij}}
\end{equation}
The main idea is to treat the short-range electron-electron
interaction with DFT, and the long range with wave function theory.
The parameter $\mu$ controls the range of the separation, and allows
to go continuously from the Kohn-Sham Hamiltonian ($\mu=0$) to
the FCI Hamiltinoan ($\mu = \infty$).

To rigorously connect wave function theory 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}^{\mathrm{lr},\mu}[n] + \bar{E}_{\mathrm{Hxc}}^{\mathrm{sr,}\mu}[n],
  \label{Fdecomp}
\end{equation}
where $n$ is a one-particle density,
$\mathcal{F}^{\mathrm{lr},\mu}$ is a long-range universal density
functional and $\bar{E}_{\mathrm{Hxc}}^{\mathrm{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} \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}
with $\hat{T}$ the kinetic energy operator, 
$\hat{W}_\mathrm{ee}^{\mathrm{lr}}$ the long-range
electron-electron interaction,
$n_\Psi$ the one-particle density associated with $\Psi$,
and $\hat{V}_{\mathrm{ne}}$ the electron-nucleus potential.
The minimizing multi-determinant 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}_{\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}$
is the complementary short-range Hartree-exchange-correlation
potential operator.
Once $\Psi^{\mu}$ has been calculated, the electronic ground-state
energy is obtained by
\begin{equation}
  \label{E-rsdft}
  E_0=  \mel{\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
range-separated DFT (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,
$w_{\mathrm{ee}}^{\mathrm{lr},\mu\to\infty}(r_{12}) = 1/r_{12}$, and
thus RS-DFT reduces to standard wave-function theory and $\Psi^\mu$ is
the FCI wave function.

\begin{figure*}
  \centering
  \includegraphics[width=0.7\linewidth]{algorithm.pdf}
  \caption{Algorithm showing the generation of the RS-DFT wave
    function.}
  \label{fig:algo}
\end{figure*}

Hence we have a continuous path connecting the KS determinant to the
FCI wave function, and as 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 using an increasing number of selected
determinants.

We can follow this path by performing FCI calculations using the
RS-DFT Hamiltonian with different values of $\mu$.  In this work, we
have used the CIPSI algorithm to peform approximate FCI calculations
with the RS-DFT Hamiltonians,\cite{GinPraFerAssSavTou-JCP-18}
$\hat{H}^\mu$ as shown in figure~\ref{fig:algo}.  In the outer loop
(red), a CIPSI selection is performed with a RS-Hamiltonian
parameterized using the current density.
An inner loop (blue) is introduced to accelerate the
convergence of the self-consistent calculation, in which the set of
determinants is kept fixed, and only the diagonalization of the
RS-Hamiltonian is performed iteratively with the updated density.
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.
Note that any range-separated post-Hartree-Fock method can be
implemented using this scheme by just replacing the CIPSI step by the
post-HF method of interest.

\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 (BFD-V$n$Z).
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 were made with \emph{Quantum
Package}.\cite{Garniron_2019,qp2_2020} We used 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}).
The convergence criterion for stopping the CIPSI calculations
was $\EPT < 1$~m\hartree{} $\vee \Ndet > 10^7$.
All the wave functions are eigenfunctions of the $S^2$ 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 Jatrow factor, as in all-electron
calculations. Simple Jastrow factors were used to reduce the
fluctuations of the local energy.




\section{Influence of the range-separation parameter on the fixed-node
  error}
\label{sec:mu-dmc}
\begin{table}
  \caption{Fixed-node energies and number of determinants in the water
    molecule and the fluorine dimer with different trial wave functions.}
  \label{tab:h2o-dmc}
  \centering
  \begin{ruledtabular}
  \begin{tabular}{crlrl}
&         \multicolumn{2}{c}{BFD-VDZ}  &           \multicolumn{2}{c}{BFD-VTZ}  \\
$\mu$     &  $\Ndet$     & $\EDMC$           &  $\Ndet$        &  $\EDMC$           \\
\hline
&         \multicolumn{4}{c}{H$_2$O}   \\
$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)$   \\
&         \multicolumn{3}{c}{F$_2$}    \\
$0.00$    &  $23$        & $-48.419\,5(4)$ \\
$0.25$    &  $8$         & $-48.421\,9(4)$ \\
$0.50$    &  $1743$      & $-48.424\,8(8)$ \\ 
$1.00$    &  $11952$     & $-48.432\,4(3)$ \\
$2.00$    &  $829438$    & $-48.441\,0(7)$ \\
$5.00$    &  $5326459$   & $-48.445(2)$ \\
$\infty$  &  $8302442$   & $-48.437(3)$ \\
  \end{tabular}
  \end{ruledtabular}
\end{table}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{h2o-dmc.pdf}
  \caption{Fixed-node energies of the water molecule for different
    values of $\mu$.}
  \label{fig:h2o-dmc}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{f2-dmc.pdf}
  \caption{Fixed-node energies of difluorine for different
    values of $\mu$.}
  \label{fig:f2-dmc}
\end{figure}
The first question we would like to address is the quality of the
nodes of the wave functions $\Psi^{\mu}$ obtained with an intermediate
range separation parameter (\textit{i.e.} $0 < \mu  < +\infty$).
We generated trial wave functions $\Psi^\mu$ with multiple values of
$\mu$, and computed the associated fixed node energy keeping all the
parameters having an impact on the nodal surface fixed.
We considered two weakly correlated molecular systems: the water
molecule and the fluorine dimer, near their equilibrium
geometry\cite{Caffarel_2016}.

From table~\ref{tab:h2o-dmc} and figures~\ref{fig:h2o-dmc}
and~\ref{fig:f2-dmc}, one can clearly observe that using FCI trial
wave functions ($\mu = \infty$) gives FN-DMC energies which are lower
than the energies obtained with a single Kohn-Sham 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 for water, and
a gain of $18 \pm 3$~m\hartree{} for F$_2$. Interestingly, using the
RS-DFT-CIPSI trial wave function with a range-separation parameter
$\mu=1.75$~bohr$^{-1}$ with the double-zeta basis one can obtain for
water a FN-DMC energy $2.6 \pm 0.7$~m\hartree{} lower than the energy
obtained with the FCI trial wave function.  This can be explained by
the inability of the basis set to properly describe the short-range
correlation effects, 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 placement of the nodes closer to the optimum.  In the case
of F$_2$, a similar behavior with a gain of $8 \pm 4$ m\hartree{} is
observed for $\mu\sim 5$~bohr$^{-1}$.
The optimal value of $\mu$ is larger than in the case of water, and this
is probably the signature of the fact that the average
electron-electron distance in the valence is smaller in F$_2$ than in
H$_2$O due to the larger nuclear charge shrinking the electron
density. At the triple-zeta level, the short-range correlations can be
better described by the determinant expansion, and the effect of 
sr-DFT on the trial wave function is insignificant on the
fixed-node energy.  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 slightly 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, i.e. at
$\mu=\infty$.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{overlap.pdf}
  \caption{Overlap of the RS-DFT CI expansion with the
    CI expansion optimized in the presence of a Jastrow factor.}
  \label{fig:overlap}
\end{figure}

This data confirms that RS-DFT-CIPSI can give improved CI coefficients
with small basis sets, similarly to the common practice of
re-optimizing the trial wave function in the presence of the Jastrow
factor. To confirm that the introduction of sr-DFT has an impact on
the CI coefficients similar to the Jastrow factor, 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 diagonalize
self-consistently the RS-DFT Hamiltonian 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$.  In figure~\ref{fig:overlap}, we plot the overlaps
$\braket{\Psi^J}{\Psi^\mu}$ obtained for water and the fluorine dimer.

In the case of H$_2$O, there is a clear maximum of overlap at
$\mu=1$~bohr$^{-1}$. This confirms that introducing short-range
correlation with DFT has the an impact on the CI coefficients similar to
the Jastrow factor. In the case of F$_2$, the Jastrow factor has
very little effect on the CI coefficients, as the overlap
$\braket{\Psi^J}{\Psi^{\mu=\infty}}$ is very close to
$1$. Nevertheless, a slight maximum is obtained for
$\mu=5$~bohr$^{-1}$.


\section{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<i} \frac{a\, r_{ij}}{1 + b\, r_{ij}}.
\end{equation}
The parameter
$a$ is determined by cusp conditions, and $b$ is obtained by energy
or variance minimization.\cite{Coldwell_1977,Umrigar_2005}
One can easily see that this parameterization of the two-body
interation is not size-consistent: the dissociation of a
diatomic molecule $AB$ with a parameter $b_{AB}$
will lead to two different two-body Jastrow factors, each
with its own optimal value $b_A$ and $b_B$. To remove the
size-consistency error on a PES using this ansätz for $J_\text{ee}$,
one needs to impose that the parameters of $J_\text{ee}$ are fixed:
$b_A = b_B = b_{AB}$.

When pseudopotentials are used in a QMC calculation, it is common
practice to localize the non-local part of the pseudopotential on the
complete wave function (determinantal component and Jastrow).
If the wave function is not size-consistent,
so will be the locality approximation. Within, the determinant
localization approximation,\cite{Zen_2019} the Jastrow factor is
removed from the wave function on which the pseudopotential is localized.
The great advantage of this approximation is that the FN-DMC energy
only depends on the parameters of the determinantal component. Using a
non-size-consistent Jastrow factor, or a non-optimal Jastrow factor will
not introduce an additional error in FN-DMC calculations, although it
will reduce the statistical errors by reducing the variance of the
local energy. Moreover, the integrals involved in the pseudo-potential
are computed analytically and the computational cost of the
pseudo-potential is dramatically reduced (for more detail, see
Ref.~\onlinecite{Scemama_2015}).

%\begin{squeezetable}
\begin{table}
  \caption{FN-DMC Energies of the fluorine atom and the dissociated fluorine
    dimer, and size-consistency error.}
  \label{tab:size-cons}
  \begin{ruledtabular}
    \begin{tabular}{cccc}
      $\mu$ &  F                &  Dissociated F$_2$  & Size-consistency error \\
      \hline
      0.00  &  $-24.188\,7(3)$  & $-48.377\,7(3)$     & $-0.000\,3(4)$ \\
      0.25  &  $-24.188\,7(3)$  & $-48.377\,2(4)$     & $+0.000\,2(5)$ \\
      0.50  &  $-24.188\,8(1)$  & $-48.376\,9(4)$     & $+0.000\,7(4)$ \\
      1.00  &  $-24.189\,7(1)$  & $-48.380\,2(4)$     & $-0.000\,8(4)$ \\
      2.00  &  $-24.194\,1(3)$  & $-48.388\,4(4)$     & $-0.000\,2(5)$ \\
      5.00  &  $-24.194\,7(4)$  & $-48.388\,5(7)$     & $+0.000\,9(8)$ \\
   $\infty$ &  $-24.193\,5(2)$  & $-48.386\,9(4)$     & $+0.000\,1(5)$ \\
    \end{tabular}
  \end{ruledtabular}
\end{table}

In this section, we make a numerical verification that the produced
wave functions are size-consistent for a given range-separation
parameter.
We have computed the energy of the dissociated fluorine dimer, where
the two atoms are at a distance of 50~\AA.  We expect that the energy
of this system is equal to twice the energy of the fluorine atom.
The data in table~\ref{tab:size-cons} shows that it is indeed the
case, so we can conclude that the proposed scheme provides
size-consistent FN-DMC energies for all values of $\mu$ (within
$2\times$ statistical error bars).


\subsection{Spin-invariance}

Closed-shell molecules often dissociate into open-shell
fragments. To get reliable atomization energies, it is important to
have a theory which is of comparable quality for open-shell and
closed-shell systems. A good test is to check that all the components
of a spin multiplet are degenerate.
FCI wave functions have this property and give degenrate energies with
respect to the spin quantum number $m_s$, but the multiplication by a
Jastrow factor introduces spin contamination if the parameters 
for the same-spin electron pairs are different from those
for the opposite-spin pairs.\cite{Tenno_2004}
Again, when pseudo-potentials are used this tiny error is transferred
in the FN-DMC energy unless the determinant localization approximation
is used.

Within DFT, the common density functionals make a difference for
same-spin and opposite-spin interactions. As DFT is a
single-determinant theory, the density functionals are designed to be
used with the highest value of $m_s$, and therefore different values
of $m_s$ lead to different energies.
So in the context of RS-DFT, the determinantal expansions will be
impacted by this spurious effect, as opposed to FCI.

\begin{table}
  \caption{FN-DMC Energies of the triplet carbon atom (BFD-VDZ) with
    different values of $m_s$.}
  \label{tab:spin}
  \begin{ruledtabular}
    \begin{tabular}{cccc}
      $\mu$ &  $m_s=1$          & $m_s=0$             & Spin-invariance error \\
      \hline
      0.00  &  $-5.416\,8(1)$  & $-5.414\,9(1)$     & $+0.001\,9(2)$ \\
      0.25  &  $-5.417\,2(1)$  & $-5.416\,5(1)$     & $+0.000\,7(1)$ \\
      0.50  &  $-5.422\,3(1)$  & $-5.421\,4(1)$     & $+0.000\,9(2)$ \\
      1.00  &  $-5.429\,7(1)$  & $-5.429\,2(1)$     & $+0.000\,5(2)$ \\
      2.00  &  $-5.432\,1(1)$  & $-5.431\,4(1)$     & $+0.000\,7(2)$ \\
      5.00  &  $-5.431\,7(1)$  & $-5.431\,4(1)$     & $+0.000\,3(2)$ \\
   $\infty$ &  $-5.431\,6(1)$  & $-5.431\,3(1)$     & $+0.000\,3(2)$ \\
    \end{tabular}
  \end{ruledtabular}
\end{table}

In this section, we investigate the impact of the spin contamination
due to the short-range density functional on the FN-DMC energy. We have
computed the energies of the carbon atom in its triplet state
with BFD pseudo-potentials and the corresponding double-zeta basis
set. The calculation was done with $m_s=1$ (3 $\uparrow$ electrons
and 1 $\downarrow$ electrons) and with $m_s=0$ (2 $\uparrow$ and 2
$\downarrow$ electrons).

The results are presented in table~\ref{tab:spin}.
Although using $m_s=0$ the energy is higher than with $m_s=1$, the
bias is relatively small, more than one order of magnitude smaller
than the energy gained by reducing the fixed-node error going from the single
determinant to the FCI trial wave function. The highest bias, close to
2~m\hartree, is obtained for $\mu=0$, but the bias decreases quickly
below 1~m\hartree{} when $\mu$ increases. As expected, with $\mu=\infty$
there is no bias (within the error bars), and the bias is not
noticeable with $\mu=5$~bohr$^{-1}$.




\subsection{Benchmark}

The 55 molecules of the benchmark for the Gaussian-1
theory\cite{Pople_1989,Curtiss_1990} were chosen to test the quality
of the RS-DFT-CIPSI trial wave functions for energy differences.



%\begin{squeezetable}
\begin{table*}
  \caption{Mean absolute error (MAE), mean signed errors (MSE) and
    standard deviations (RMSD) obtained with the different methods and
    basis sets.}
  \label{tab:mad}
  \begin{ruledtabular}
    \begin{tabular}{ll rrr rrr rrr}
Method           &  \(\mu\)     &  \phantom{}  &  VDZ-BFD                &  \phantom{}  &  \phantom{}  &  VTZ-BFD    &  \phantom{}  &  \phantom{}  &  VQZ-BFD    &  \phantom{}  \\
\phantom{}       &  \phantom{}  &  MAE         &  MSE                    &  RMSD          &  MAE         &  MSE        &  RMSD          &  MAE         &  MSE        &  RMSD          \\
\hline
PBE              &  0           &  5.02        &  -3.70                  &  6.04        &  4.57        &  1.00       &  5.32        &  5.31        &  0.79       &  6.27        \\
BLYP             &  0           &  9.53        &  -9.21                  &  7.91        &  5.58        &  -4.44      &  5.80        &  5.86        &  -4.47      &  6.43        \\
PBE0             &  0           &  11.20       &  -10.98                 &  8.68        &  6.40        &  -5.78      &  5.49        &  6.28        &  -5.65      &  5.08        \\
B3LYP            &  0           &  11.27       &  -10.98                 &  9.59        &  7.27        &  -5.77      &  6.63        &  6.75        &  -5.53      &  6.09        \\
\hline
CCSD(T)          &  \(\infty\)  &  24.10       &  -23.96                 &  13.03       &  9.11        &  -9.10      &  5.55        &  4.52        &  -4.38      &  3.60        \\
\hline
RS-DFT-CIPSI      &  0           &  4.53        &  -1.66                  &  5.91        &  6.31        &  0.91       &  7.93        &  6.35        &  3.88       &  7.20        \\
\phantom{}       &  1/4         &  5.55        &  -4.66                  &  5.52        &  4.58        &  1.06       &  5.72        &  5.48        &  1.52       &  6.93        \\
\phantom{}       &  1/2         &  13.42       &  -13.27                 &  7.36        &  6.77        &  -6.71      &  4.56        &  6.35        &  -5.89      &  5.18        \\
\phantom{}       &  1           &  17.07       &  -16.92                 &  9.83        &  9.06        &  -9.06      &  5.88        &  ---         &  ---        &  ---         \\
\phantom{}       &  2           &  19.20       &  -19.05                 &  10.91       &  ---         &  ---        &  ---         &  ---         &  ---        &  ---         \\
\phantom{}       &  5           &  22.93       &  -22.79                 &  13.24       &  ---         &  ---        &  ---         &  ---         &  ---        &  ---         \\
\phantom{}       &  \(\infty\)  &  23.62       &  -23.48                 &  12.81       &  ---         &  ---        &  ---         &  ---         &  ---        &  ---         \\
\hline
DMC@RS-DFT-CIPSI  &  0           &  4.61(34)    &  -3.62(\phantom{0.}34)  &  5.30        &  3.52(19)    &  -1.03(19)  &  4.39        &  3.16(26)    &  -0.12(26)  &  4.12        \\
\phantom{}       &  1/4         &  4.04(37)    &  -3.13(\phantom{0.}37)  &  4.88        &  3.39(77)    &  -0.59(77)  &  4.44        &  2.90(25)    &  0.25(25)   &  3.74        \\
\phantom{}       &  1/2         &  3.74(35)    &  -3.53(\phantom{0.}35)  &  4.03        &  2.46(18)    &  -1.72(18)  &  3.02        &  2.06(35)    &  -0.44(35)  &  2.74        \\
\phantom{}       &  1           &  5.42(29)    &  -5.14(\phantom{0.}29)  &  4.55        &  4.38(94)    &  -4.24(94)  &  5.11        &   ---        &   ---       &   ---        \\
\phantom{}       &  2           &  5.98(83)    &  -5.91(\phantom{0.}83)  &  4.79        &   ---        &   ---       &   ---        &   ---        &   ---       &   ---        \\
\phantom{}       &  5           &  6.18(84)    &  -6.13(\phantom{0.}84)  &  4.87        &   ---        &   ---       &   ---        &   ---        &   ---       &   ---        \\
\phantom{}       &  \(\infty\)  &  7.38(1.08)  &  -7.38(1.08)            &  5.67        &   ---        &   ---       &   ---        &   ---        &   ---       &   ---        \\
\phantom{}       &  Opt.        &  5.85(1.75)  &  -5.63(1.75)            &  4.79        &   ---        &   ---       &   ---        &   ---        &   ---       &   ---        \\
    \end{tabular}
  \end{ruledtabular}
\end{table*}
%\end{squeezetable}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{g2-dmc.pdf}
  \caption{Errors in the DMC atomization energies with the different
    trial wave functions. Each dot corresponds to an atomization
    energy.
    The boxes contain the data between first and third quartiles, and
    the line in the box represents the median. The outliers are shown
    with a cross.}
  \label{fig:g2-dmc}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{g2-ndet.pdf}
  \caption{Number of determinants in the different trial wave
    functions. Each dot corresponds to an atomization energy.
    The boxes contain the data between first and third quartiles, and
    the line in the box represents the median. The outliers are shown
    with a cross.}
  \label{fig:g2-ndet}
\end{figure}

The number of determinants in the wave functions are shown in
figure~\ref{fig:n2-ndet}. For all the calculations, the stopping
criterion of the CIPSI algorithm was $\EPT < 1$~m\hartree{} or $\Ndet >
10^7$.
For FCI, we have given extrapolated values at $\EPT\rightarrow 0$.
At $\mu=0$ the number of determinants is not equal to one because
we have used the natural orbitals of a first CIPSI calculation.
So the Kohn-Sham determinant is expressed as a linear combination of
determinants built with natural orbitals. Note that it is possible to add
an extra step to the algorithm to compute the natural orbitals from the
RS-DFT/CIPSI wave function, and re-do the RS-DFT/CIPSI calculation with
these orbitals to get an even more compact expansion. In that case, we would
have obtained the Kohn-Sham orbitals with $\mu=0$, and the solution would have
been the KS single determinant.

We could have obtained
single-determinant wave functions by using the natural orbitals of a first


\section{Conclusion}

We have seen that introducing short-range correation via
a range-separated Hamiltonian in a full CI expansion yields improved
nodes, especially with small basis sets. The effect is similar to the
effect of re-optimizing the CI coefficients in the presence of a
Jastrow factor, but without the burden of performing a stochastic
optimization.
The proposed procedure provides a method to optimize the
FN-DMC energy via a single parameter, namely the range-separation
parameter $\mu$. The size-consistency error is controlled, as well as the
invariance with respect to the spin projection $m_s$.
Finding the optimal value of $\mu$ gives the lowest FN-DMC energies
within basis set. However, if one wants to compute an energy
difference, one should not minimize the
FN-DMC energies of the reactants independently. It is preferable to
choose a value of $\mu$ for which the fixed-node errors are well
balanced, leading to a good cncellation of errors. We found that a
value of $\mu=0.5$~bohr${^-1}$ is the value where the errors are the
smallest. Moreover, such a small value of $\mu$ gives extermely
compact wave functions, making this recipe a good candidate for
accurate calcultions of large systems with a multi-reference character.



%%---------------------------------------
\begin{acknowledgments}
An award of computer time was provided by the Innovative and Novel
Computational Impact on Theory and Experiment (INCITE) program. This
research has used resources of the Argonne Leadership Computing
Facility, which is a DOE Office of Science User Facility supported
under Contract DE-AC02-06CH11357.  AB, was supported by the
U.S. Department of Energy, Office of Science, Basic Energy Sciences,
Materials Sciences and Engineering Division, as part of the
Computational Materials Sciences Program and Center for Predictive
Simulation of Functional Materials.
\end{acknowledgments}


\bibliography{rsdft-cipsi-qmc}

\begin{enumerate}
 \item Total energies and nodal quality:
  \begin{itemize}
   \item Facts: KS occupied orbitals closer to NOs than HF
   \item Even if exact functional, complete basis set, still approximated nodes for KS
   \item KS -> exponentially fast convergence (as HF) with basis because of non divergence of effective KS potential (citer le papier de Gill)
   \item With correlation consistent basis set, FCI nodes (which include correlation) are better than KS
   \item With FCI, good limit at CBS ==> exact energy
   \item But slow convergence with basis set because of divergence of the e-e interaction not well represented in atom centered basis set
   \item Exponential increase of number of Slater determinants
   \item Cite papiers RS-DFT: there exists an hybrid scheme combining fast convergence wr to basis set (non divergent basis set) and short expansion in SCI (cite papier Ferté)
   \item Question: does such a scheme provide better nodal quality ?
   \item In RS-DFT we cannot optimize energy with $\mu$ , but in FNDMC
   \item Factual stuffs: with optimal $\mu$, lower FNDMC energy than HF/KS/FCI
    \begin{itemize}
     \item  less determinants $\Rightarrow$ large systems
     \item  only one parameter to optimize $\Rightarrow$ deterministic
     \item $\Rightarrow$ reproducible
    \end{itemize}
   \item with the optimal $\mu$:
    \begin{itemize}
     \item Direct optimization of FNDMC with one parameter
     \item Do we improve energy differences ?
     \item system dependent
     \item basis set dependent: $\mu \rightarrow \infty$ when $\mathcal{B}\rightarrow \text{CBS}$
     \item large wave functions
    \end{itemize}
    \item Invariance with $m_s$

  \end{itemize}
\end{enumerate}



\end{document}


% * Recouvrement avec Be
% * Tester sr-LDA avec H2O-DZ
% * Manu doit faire des programmes pour des plots de ensite a 1 et 2
%   corps le long des axes de liaison, et l'integrale de la densite a
%   2 corps a coalescence.
% 1 Manu calcule Be en cc-pvdz tous electrons: FCI -> NOs -> FCI ->
%   qp edit -n 200
% 2 Manu calcule qp_cipsi_rsh avec mu = [ 1.e-6 , 0.25, 0.5, 1.0, 2.0, 5.0, 1e6 ]
% 3 Manu fait tourner es petits programmes
% 4 Manu envoie a toto un tar avec tous les ezfio
% 5 Toto optimise les coefs en presence e jastrow
% 6 Toto renvoie a manu psicoef
% 7 Manu fait tourner ses petits programmes avec psi_J