377 lines
17 KiB
TeX
377 lines
17 KiB
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}}
|
|
\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}
|
|
|
|
|
|
\section{Range-separated DFT in a nutshell}
|
|
\label{sec:rsdft}
|
|
|
|
\subsection{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}
|
|
|
|
|
|
|
|
|
|
|
|
\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{g09} 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{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}.
|
|
|
|
|
|
\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<i} \frac{a r_{ij}}{1 + b r_{ij}}.
|
|
\end{equation}
|
|
$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
|
|
heteroatomic diatomic molecule $AB$ with a parameter $b_{AB}$
|
|
will lead to two different two-body Jastrow factors on each atom, each
|
|
with its own optimal value $b_A$ and $b_B$. To remove the
|
|
size-consistence 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 pseudopotential on the complete wave
|
|
function. If the wave function is not size-consistent, so will be the
|
|
locality approximation. Recently, the determinant localization
|
|
approximation was introduced.\cite{Zen_2019} This approximation
|
|
consists in removing the Jastrow factor from the wave function on
|
|
which the pseudopotential is localized.
|
|
The great advantage of this approximation is that the FN-DMC energy
|
|
within this approximation only depends on the parameters of the
|
|
determinantal component. Using a size-inconsistent Jastrow factor, or
|
|
a non-optimal Jastrow factor will not introduce an additional
|
|
size-consistence error in FN-DMC calculations, although it will
|
|
reduce the statistical errors by reducing the variance of the local energy.
|
|
|
|
|
|
The energy computed within density functional theory is extensive, 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, the convergence of the FCI energies to the CBS
|
|
limit will be slower because of the description of short-range electron
|
|
correlation with atom-centered functions, but ultimately the exact
|
|
energy will be reached.
|
|
|
|
|
|
The 55 molecules of the benchmark for the Gaussian-1
|
|
theory\cite{Pople_1989,Curtiss_1990} were chosen to test the quality
|
|
of the RSDFT-CIPSI trial wave functions for energy differences.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
%%---------------------------------------
|
|
\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}
|
|
|
|
\end{document}
|