224 lines
16 KiB
TeX
224 lines
16 KiB
TeX
%****************************************************************
|
|
\section{Self-consistent electron-nucleus cusp correction for molecular orbitals}
|
|
%****************************************************************
|
|
At short interparticle distances, the Coulombic singularity dominates all other terms and, near the two-particle coalescence point, the behaviour of the exact wave function $\Psi$ becomes independent of other details of the system \cite{eee15}.
|
|
In particular, early work by Kato \cite{Kato51, Kato57}, and elaborations by Pack and Byers Brown \cite{Pack66}, showed that, as one electron at $\bri$ approaches a nucleus of charge $Z_A$ at $\brA$, we have
|
|
\begin{equation}
|
|
\label{eq:eNcusp}
|
|
\eval{\pdv{\expval{\Psi(\bR)}}{\ri}}_{\ri = \rA} = -Z_A \left.\expval{\Psi(\bR)}\right|_{\ri = \rA},
|
|
\end{equation}
|
|
where $\left.\expval{\Psi(\bR)}\right|_{\ri = \rA}$ is the spherical average of $\Psi(\bR)$ about $\bri = \brA$.
|
|
|
|
To remove divergences in the local energy at the electron-nucleus (e-n) coalescence points, cusp conditions such as \eqref{eq:eNcusp} must be satisfied.\footnote{Note that the use of pseudopotentials would also remove the divergence at the e-n coalescence points as routinely done in QMC calculations.}
|
|
These divergences are especially harmful in DMC calculations, where they can lead to a large increase of the statistical variance, population-control problems and significant biases \cite{Drummond04}.
|
|
|
|
There are two possible ways to enforce the correct e-n cusp.
|
|
One way to do it is to enforce the e-n cusp within the Jastrow factor in Eq.~\eqref{eq:Psi-trial}.
|
|
This has the disadvantage of increasing the number of parameters in $\Js(\bR)$, and their interdependence can be tricky as one must optimise the large number of linear and non-linear parameters contained in $\Js(\bR)$ via a stochastic (noisy) optimization of the energy and/or its variance.
|
|
However, it is frequently done in the literature thanks to some recent progress \cite{Toulouse07, Umrigar07, Toulouse08}.
|
|
|
|
Here, we will follow an alternative path which consists in imposing the e-n cusp within the multideterminant expansion of Eq.~\eqref{eq:Psi-trial}.
|
|
However, because one usually employs Gaussian basis functions \cite{Boys69} (as in standard quantum chemistry packages), the MOs $\MO{i}(\br)$ are cuspless, i.e.~
|
|
\begin{equation}
|
|
\eval{\pdv{\expval{\MO{i}(\br)}}{r}}_{r = \rA} = 0.
|
|
\end{equation}
|
|
One solution would be to use a different set of basis functions \cite{McKemmish14} as, for instance, Slater basis functions \cite{Slater30, Bouferguene96}.
|
|
However, they are known to be troublesome, mainly due to the difficulty of calculating multicentric two-electron integrals which require expensive numerical expansions or quadratures.
|
|
Nevertheless, some authors \cite{Nemec10} have explored using wave functions built with Slater basis functions \cite{Lenthe03} while imposing the right e-n cusp afterwards.\footnote{It is also possible to enforce the correct e-n cusp during the SCF process although it is rarely done \cite{Galek05}.}
|
|
These types of calculations can be performed with an electronic structure package such as ADF \cite{ADF}.
|
|
However, as far as we know, it is hard to perform large-scale calculations with Slater basis functions and the virtual space is usually poorly described.
|
|
Moreover, Gaussian basis are usually of better quality than Slater-based ones due to the extensive knowledge and experience gathered by quantum chemists over the last fifty years while building robust, compact and effective Gaussian basis sets \cite{Huzinaga65, Hehre69, Hehre72, Schmidt79, Feller79, Almlof87, Almlof90, Bauschlicher93, Dunning89, Woon93, Peterson93, Jensen01, Jensen07, Jensen12}.
|
|
|
|
Conventional cusp correction methods usually replace the part of $\GTO{\mu}(\br)$ or $\MO{i}(\br)$ close to the nuclei within a cusp-correction radius by a polynomial or a spline function which fulfils Kato cusp conditions and ensures a well-behaved local energy \cite{Manten01, Ma05, Kussmann07, Per08}.
|
|
For atoms, one can also substitute Gaussian core orbitals by tabulated Slater-based ones \cite{Bunge93, Scemama14}.
|
|
In the same vein, Toulouse and Umrigar have fitted Slater basis functions with a large number of Gaussian functions and replaced them within the QMC calculation \cite{Toulouse08}.
|
|
However, it is hardly scalable for large systems due to its lack of transferability and the ever-growing number of primitive two-electron integrals to compute.
|
|
|
|
Here, we propose to follow a different, alternative route by augmenting conventional Gaussian basis sets with cusp-correcting Slater basis functions.
|
|
Mixed Gaussian-Slater basis sets have been already considered in the past with limited success due to the difficultly of computing efficiently mixed electron repulsion integrals \cite{Allen59, Silver70, Bacskay72a, Bacskay72b, Bacskay72c, Bacskay72d, Bugaets76}.
|
|
However, we will show that, because of the way we introduce the cusp correction, the integrals required here are not that scary.
|
|
For the sake of simplicity, we will focus on the HF formalism in the present study, although our scheme can also be applied to KS calculations.
|
|
|
|
%----------------------------------------------------------------
|
|
\subsection{Cusp-corrected orbitals}
|
|
%----------------------------------------------------------------
|
|
Assuming that the Jastrow factor does not contribute, a sufficient condition to ensure that $\PsiT$ (cf Eq.~\eqref{eq:Psi-trial}) fulfils the e-n cusp \eqref{eq:eNcusp} is that each (occupied and virtual) MO $\ccedMO{i}(\br)$ satisfies the e-n cusp at each nuclear position $\brA$:
|
|
\begin{equation}
|
|
\label{eq:eNcusp-orb}
|
|
\eval{\pdv{\expval*{\ccedMO{i}(\br)}}{r}}_{r = \rA} = -Z_A \expval*{\ccedMO{i}(\br)}\big|_{r = \rA}.
|
|
\end{equation}
|
|
We also assume that the basis functions have been already orthogonalised via the standard procedure \cite{Szabo}, i.e.~$\braket{\GTO{\mu}}{\GTO{\nu}} = \delta_{\mu \nu}$, where $\delta_{\mu \nu}$ is the Kronecker delta \cite{NISTbook}.
|
|
|
|
Here, we enforce the correct value of the e-n cusp by adding a cusp-correcting orbital to each MO:
|
|
\begin{equation}
|
|
\label{eq:tphi}
|
|
\ccedMO{i}(\br) = \phi_i(\br) + \POp \ccingMO{i}(\br),
|
|
\end{equation}
|
|
with
|
|
\begin{equation}
|
|
\ccingMO{i}(\br) = \sum_A^M \cSTO{A}{i} \, \STO{A}{i}(\br),
|
|
\end{equation}
|
|
where $M$ is the number of nuclear centres and
|
|
\begin{equation}
|
|
\label{eq:Slater}
|
|
\STO{A}{i}(\br) = \sqrt{\frac{\expSTOa_{i}^3}{\pi}} \exp[-\expSTOa_{i} \abs{\br - \brA}]
|
|
\end{equation}
|
|
is a $s$-type Slater function centred on nucleus $A$ with an orbital-dependent exponent $\expSTOa_{i}$.
|
|
In Eq.~\eqref{eq:tphi}, the projector
|
|
\begin{equation}
|
|
\POp = \IdOp - \sum_\mu \dyad{\GTO{\mu}}
|
|
\end{equation}
|
|
(where $\IdOp$ is the identity operator) ensures orthogonality between $\MO{i}(\br)$ and the cusp-correcting orbital $\ccingMO{i}(\br)$.
|
|
|
|
It is easy to show that ensuring the right e-n cusp yields the following linear system of equations for the coefficients $c_{Ai}$:
|
|
\begin{equation}
|
|
\label{eq:cAi-lineq}
|
|
\sum_B \Bigg[ - \frac{\delta_{AB}}{Z_A} \partial_r \STO{A}{i}(\brA) - \STO{B}{i}(\brA)
|
|
\\
|
|
+ \sum_\mu \DrOvEl{B}{\mu}{i} \GTO{\mu} (\brA) \Bigg] \cSTO{B}{i} = \MO{i}(\brA),
|
|
\end{equation}
|
|
where the explicit expression of the matrix elements $\DrOvEl{\mu}{A}{i} = \braket{\GTO{\mu}}{\STO{A}{i}}$ can be easily found using standard procedures and
|
|
\begin{equation}
|
|
\partial_r \STO{A}{i}(\brA) \equiv \eval{\pdv{\STO{A}{i}(\br)}{r}}_{r = \rA}.
|
|
\end{equation}
|
|
Equation \eqref{eq:cAi-lineq} can be easily solved using standard linear algebra packages, and provides a way to obtain a cusp-corrected orbital $\ccedMO{i}(\br)$ from a given MO $\MO{i}(\br)$.
|
|
For reasons that will later become apparent, we will refer to this procedure as a one-shot (OS) calculation.
|
|
Next, we are going to explain how one can optimise self-consistently the coefficients $\cSTO{A}{i}$.
|
|
|
|
%----------------------------------------------------------------
|
|
\subsection{Self-consistent dressing of the Fock matrix}
|
|
%----------------------------------------------------------------
|
|
So far, the coefficient $\cSTO{A}{i}$ have been set via Eq.~\eqref{eq:cAi-lineq}.
|
|
Therefore, they have not been obtained via a variational procedure as their only purpose is to enforce the e-n cusp.
|
|
However, they do depend on $\phi_i(\brA)$, hence on the MO coefficients $\cGTO{\mu}{i}$.
|
|
We will show below that one can optimise simultaneously the coefficients $\cSTO{A}{i}$ and $\cGTO{\mu}{i}$ by constructing an orbital-dependent effective Fock matrix.
|
|
|
|
As it is ultimately what we wish for, the key point is to assume that $\ccedMO{i}(\br)$ is an eigenfunction of the Fock operator $\FkOp$, i.e.
|
|
\begin{equation}
|
|
\label{eq:Jia}
|
|
\FkOp \ket*{\ccedMO{i}} = \DrMOeigval{i} \ket*{\ccedMO{i}}.
|
|
\end{equation}
|
|
Note that, even at convergence of a conventional HF or KS calculation, the equality \eqref{eq:Jia} is never fulfilled (unless the basis happens to span the exact orbital).
|
|
This under-appreciated fact has been used by Deng et al.~to design a measure of the quality of a MO \cite{Deng10}.
|
|
|
|
Next, we project out Eq.~\eqref{eq:Jia} over $\bra*{\chi_\mu}$ yielding
|
|
\begin{equation}
|
|
\label{eq:HFeq}
|
|
\sum_\nu F_{\mu \nu} \cGTO{\nu}{i} + \sum_A \cSTO{A}{i} \qty( \DrFkEl{\mu}{A}{i} - \sum_\lambda \FkEl{\mu}{\lambda} \DrOvEl{\lambda}{A}{i} )
|
|
= \DrMOeigval{i} \cGTO{\mu}{i},
|
|
\end{equation}
|
|
where
|
|
\begin{align}
|
|
\FkEl{\mu}{\nu} & = \mel{\GTO{\mu}}{\FkOp}{\GTO{\nu}},
|
|
&
|
|
\DrFkEl{\mu}{A}{i} & = \mel{\GTO{\mu}}{\FkOp}{\STO{A}{i}}.
|
|
\end{align}
|
|
In the general case, because we must use basis functions with non-zero derivatives at the nucleus, finding the matrix elements $\DrFkEl{\mu}{A}{i}$ is challenging and costly.
|
|
However, because we are interested in the e-n cusp, we have found that a satisfactory approximation is
|
|
\begin{equation}
|
|
\label{eq:approx}
|
|
\DrFkEl{\mu}{A}{i} - \sum_\lambda \FkEl{\mu}{\lambda} \DrOvEl{\lambda}{A}{i} \approx \DrHcEl{\mu}{A}{i} - \sum_\lambda \HcEl{\mu}{\lambda} \DrOvEl{\lambda}{A}{i}
|
|
\end{equation}
|
|
where
|
|
\begin{align}
|
|
\HcEl{\mu}{\nu} & = \mel{\GTO{\mu}}{\HcOp}{\GTO{\nu}},
|
|
&
|
|
\DrHcEl{\mu}{A}{i} & = \mel{\GTO{\mu}}{\HcOp}{\STO{A}{i}},
|
|
\end{align}
|
|
and $\HcOp$ is the core Hamiltonian.
|
|
Note that, in Eq.~\eqref{eq:approx}, it is important to use the same approximation for both terms ($\DrFkEl{\mu}{A}{i} \approx \DrHcEl{\mu}{A}{i}$ and $\FkEl{\mu}{\nu} \approx \HcEl{\mu}{\nu}$) in order to preserve the subtle balance between the two terms.
|
|
For the sake of clarity, the expression of the matrix elements $\DrHcEl{\mu}{A}{i}$ will not be reported here, but they can be obtained easily realising that these integrals only involve $s$-type Slater functions.
|
|
|
|
The eigenvalue problem given by Eq.~\eqref{eq:HFeq} can be recast as
|
|
\begin{equation}
|
|
\sum_\nu \DrFkOpEl{\mu}{\nu}{i} \cGTO{\nu}{i} = \DrMOeigval{i} \cGTO{\mu}{i},
|
|
\end{equation}
|
|
where we have ``dressed'' the diagonal of the Fock matrix
|
|
\begin{equation}
|
|
\label{eq:DressedFock}
|
|
\DrFkOpEl{\mu}{\nu}{i} =
|
|
\begin{cases}
|
|
\FkEl{\mu}{\mu} + \DrEl{\mu}{i},
|
|
&
|
|
\text{if $\mu = \nu$},
|
|
\\
|
|
\FkEl{\mu}{\nu},
|
|
&
|
|
\text{otherwise},
|
|
\end{cases}
|
|
\end{equation}
|
|
with
|
|
\begin{equation}
|
|
\label{eq:dressing}
|
|
\DrEl{\mu}{i} = \cGTO{\mu}{i}^{-1} \sum_A \cSTO{A}{i} \qty( \DrHcEl{\mu}{A}{i} - \sum_\lambda \HcEl{\mu}{\lambda} \DrOvEl{\lambda}{A}{i} ).
|
|
\end{equation}
|
|
The process is repeated until our convergence criteria is met, i.e.~the largest absolute value of the elements of the commutator $\bDrFkOp{i} \bP - \bP \bDrFkOp{i}$ is lower than a given threshold, where $\bDrFkOp{i}$ is the dressed Fock matrix (Eq.~\eqref{eq:DressedFock}) and $\bP$ is the density matrix with
|
|
\begin{equation}
|
|
P_{\mu \nu} = \sum_i^\occ \cGTO{\mu}{i} \cGTO{\nu}{i}.
|
|
\end{equation}
|
|
We will refer to this procedure as self-consistent dressing (SCD).
|
|
|
|
Similarly to the Perdew-Zunger self-interaction correction \cite{PZ81}, the orbitals $\ccedMO{i}(\br)$ are eigenfunctions of different Fock operators and therefore no longer necessarily orthogonal.
|
|
Practically, we have found that the e-n cusp correction makes the cusp-corrected MOs $\ccedMO{i}$ slightly non-orthogonal.
|
|
However, this is not an issue as, within QMC, one evaluates the energy via MC sampling which only requires the evaluation of the MOs and their first and second derivatives.
|
|
|
|
Obviously, as evidenced by Eq.~\eqref{eq:DressedFock}, when $\cGTO{\mu}{i}$ is small, the dressing of the Fock matrix is numerically unstable.
|
|
Therefore, we have chosen not to dress the Fock matrix if $\cGTO{\mu}{i}$ is smaller than a user-defined threshold $\tau$.
|
|
We have found that a value of $10^{-5}$ is suitable for our purposes, and we use the same value for the convergence threshold.
|
|
Moreover, we have found that setting \cite{Ma05}
|
|
\begin{equation}
|
|
\label{eq:Zeff}
|
|
\expSTOa_i = \frac{\MO{i}(\brA)}{\MOsA{i}(\brA)} Z_A
|
|
\end{equation}
|
|
(where $\MOsA{i}(\br)$ corresponds to the $s$-type components of $\MO{i}(\br)$ centred at $\brA$) yields excellent results for molecular systems.
|
|
In the case where $\MOsA{i}(\brA) = 0$, the MO is effectively zero at $\br = \brA$ and, therefore, does not need to be cusp corrected.
|
|
|
|
%----------------------------------------------------------------
|
|
\subsection{An illustrative example}
|
|
%----------------------------------------------------------------
|
|
%%% FIGURE 1 %%%
|
|
\begin{figure}
|
|
\centering
|
|
\includegraphics[width=0.6\linewidth]{../Chapter5/fig/fig1-en}
|
|
\caption{
|
|
\label{fig:H-EL}
|
|
Local energy $E_\text{L}(r)$ for various wave functions of the \ce{H} atom.
|
|
The cuspless wave function is obtained with the decontracted STO-3G Gaussian basis set (red curve), and the OS and SCD cusp-corrected wave functions (blue and orange curves respectively) are obtained using $\expSTOa_\text{H} = 1$.}
|
|
\end{figure}
|
|
%%% %%%
|
|
|
|
%%% TABLE 1 %%%
|
|
\begin{table}
|
|
\centering
|
|
\caption{
|
|
\label{tab:H}
|
|
Variational energy and its variance for various wave functions of the \ce{H} atom.
|
|
The energy is obtained with the decontraced STO-3G Gaussian basis set.
|
|
The OS and SCD cusp-corrected energies are obtained by adding a Slater basis function of unit exponent ($\expSTOa_\text{H} = 1$) to the Gaussian basis set.
|
|
The energy and variance at each iteration of the SCD process is also reported.}
|
|
\begin{tabular}{cccccc}
|
|
\hline
|
|
Basis & Cusp correction & Iteration & Energy & Variance \\
|
|
\hline
|
|
Gaussian & \cdash & & $-0.495741$ & $2.23 \times 10^{-1}$ \\
|
|
Mixed & OS & & $-0.499270$ & $4.49 \times 10^{-2}$ \\
|
|
Mixed & SCD & \#1 & $-0.499270$ & $4.49 \times 10^{-2}$ \\
|
|
& & \#2 & $-0.499970$ & $3.07 \times 10^{-6}$ \\
|
|
& & \#3 & $-0.500000$ & $4.88 \times 10^{-9}$ \\
|
|
\hline
|
|
\end{tabular}
|
|
\end{table}
|
|
%%% %%%
|
|
|
|
Let us illustrate the present method with a simple example!
|
|
For pedagogical purposes, we have computed the wave function of the hydrogen atom within a small Gaussian basis (decontracted STO-3G basis).
|
|
In Fig.~\ref{fig:H-EL}, we have plotted the local energy associated with this wave function as well as its OS and SCD cusp-corrected versions.
|
|
The numerical results are reported in Table \ref{tab:H}.
|
|
As expected, the ``cuspless'' local energy (red curve) diverges for small $r$ with a variational energy off by $4.3$ millihartree compared to the exact value of $-1/2$.
|
|
The OS cusp-correcting procedure which introduces a Slater basis function of unit exponent (but does not re-optimise any coefficients) cures the divergence at $r = 0$ and significantly improves (by roughly one order of magnitude) both the variational energy and the variance.
|
|
Moreover, we observe that the long-range part of the wave function is also improved compared to the Gaussian basis set due to the presence of the Slater basis function which has the correct asymptotic decay.
|
|
The SCD cusp-correcting procedure further improve upon the OS scheme, and we reach a variance lower than $10^{-8}$ after only 3 iterations.
|
|
|
|
Obviously, more tests are required to ensure the validity of the method but we believe that the present results are encouraging and we hope to pursuit this study in the near future.
|
|
|