diff --git a/Manuscript/srDFT_SC.tex b/Manuscript/srDFT_SC.tex index 8cd4ea4..3758a50 100644 --- a/Manuscript/srDFT_SC.tex +++ b/Manuscript/srDFT_SC.tex @@ -463,7 +463,7 @@ is the valence-only effective interaction and One would note the restrictions of the sums to the set of active orbitals in Eqs.~\eqref{eq:fbasis_val} and \eqref{eq:twordm_val}. It is also noteworthy that, with the present definition, $\wbasisval$ still tends to the usual Coulomb interaction as $\Bas \to \CBS$. -\subsection{General form of the complementary density functional} +\subsection{General form of the complementary functional} \label{sec:functional} \subsubsection{Generic approximate form} @@ -531,12 +531,12 @@ Second, $\efuncdenpbe{\argebasis}$ correctly vanishes for systems with uniformil \end{equation} This property is doubly guaranteed by i) the choice of setting the effective interaction $\wbasis$ at $\infty$ for a vanishing pair density [see Eq.~\eqref{eq:wbasis}] leading to $\mu(\br{}) \to \infty$ and thus a vanishing $\ecmd(\argecmd)$ according to Eq.~\eqref{eq:lim_muinf}, and ii) the fact that $\ecmd(\argecmd)$ vanishes anyway when the on-top pair density vanishes [see Eq.~\eqref{eq:lim_n2}]. -\subsection{Requirements on the complementary density functional for strong correlation} +\subsection{Requirements on the complementary functional for strong correlation} \label{sec:requirements} -An important requirement for any electronic-structure method is size consistency, \ie$\,$ the additivity of the energies of non-interacting fragments, which is mandatory to avoid any ambiguity in computing interaction energies. When two subsystems \ce{A} and \ce{B} dissociate in closed-shell systems, as in the case of weak intermolecular interactions for instance, spin-restricted Hartree-Fock (RHF) is size-consistent. When the two subsystems dissociate in open-shell systems, such as in covalent bond breaking, it is well known that the RHF approach fails and an alternative is to use a complete-active-space self-consistent-field (CASSCF) wave function which, provided that the active space has been properly chosen, leads to additive energies. +An important requirement for any electronic-structure method is size consistency, \ie, the additivity of the energies of non-interacting fragments, which is mandatory to avoid any ambiguity in computing interaction energies. When two subsystems \ce{A} and \ce{B} dissociate in closed-shell systems, as in the case of weak intermolecular interactions for instance, spin-restricted Hartree-Fock (RHF) is size-consistent. When the two subsystems dissociate in open-shell systems, such as in covalent bond breaking, it is well known that the RHF approach fails and an alternative is to use a complete-active-space self-consistent-field (CASSCF) wave function which, provided that the active space has been properly chosen, leads to additive energies. -Another important requirement is spin-multiplet degeneracy, \ie$\,$ the independence of the energy with respect to the $S_z$ component of a given spin state, which is also a property of any exact wave function. Such a property is also important in the context of covalent bond breaking where the ground state of the supersystem $\ce{A + B}$ is generally of lower spin than the corresponding ground states of the fragments (\ce{A} and \ce{B}) which can have multiple $S_z$ components. +Another important requirement is spin-multiplet degeneracy, \ie, the independence of the energy with respect to the $S_z$ component of a given spin state, which is also a property of any exact wave function. Such a property is also important in the context of covalent bond breaking where the ground state of the supersystem $\ce{A + B}$ is generally of lower spin than the corresponding ground states of the fragments (\ce{A} and \ce{B}) which can have multiple $S_z$ components. \subsubsection{Spin-multiplet degeneracy} @@ -572,13 +572,13 @@ We refer the interested reader to the {\SI} for a detailed proof of the latter s In the case where the two subsystems \ce{A} and \ce{B} dissociate in closed-shell systems, a simple RHF wave function ensures this property, but when one or several covalent bonds are broken, a properly chosen CASSCF wave function is sufficient to recover this property. The underlying active space must however be chosen in such a way that it leads to size-consistent energies in the limit of dissociated fragments. -\subsection{Actual approximations used for the complementary density functional} +\subsection{Actual approximations used for the complementary functional} \label{sec:def_func} %\subsubsection{Definition of the protocol to design functionals} As the present work focuses on the strong-correlation regime, we propose here to investigate only approximate functionals which are $S_z$ independent and size-consistent in the case of covalent bond breaking. Therefore, the wave functions $\psibasis$ used throughout this paper are CASSCF wave functions in order to ensure size consistency of all local quantities. The difference between two flavors of functionals are only due to the type of i) spin polarization, and ii) on-top pair density. -Regarding the spin polarization that enters into $\varepsilon_{\text{c}}^{\text{PBE}}(\argepbe)$, two different types of $S_z$-independent formulations are considered: i) the \textit{effective} spin polarization $\tilde{\zeta}$ defined in Eq.~\eqref{eq:def_effspin-0} and calculated from the CASSCF wave function, and ii) a \textit{zero} spin polarization. In the latter case, the functional is referred as to ``SU'' which stands for ``spin unpolarized''. +Regarding the spin polarization that enters into the function $\varepsilon_{\text{c}}^{\text{PBE}}(\argepbe)$, two different types of $S_z$-independent formulations are considered: i) the \textit{effective} spin polarization $\tilde{\zeta}$ defined in Eq.~\eqref{eq:def_effspin-0} and calculated from the CASSCF wave function, and ii) a \textit{zero} spin polarization. In the latter case, the functional is referred as to ``SU'' which stands for ``spin unpolarized''. Regarding the on-top pair density entering in Eq.~\eqref{eq:def_beta}, we use two different approximations. The first one is based on the uniform electron gas (UEG) and reads \begin{equation} @@ -587,14 +587,14 @@ Regarding the on-top pair density entering in Eq.~\eqref{eq:def_beta}, we use tw \end{equation} where the pair-distribution function $g_0(n)$ is taken from Eq.~(46) of Ref.~\onlinecite{GorSav-PRA-06}. As the spin polarization appears in Eq.~\eqref{eq:def_n2ueg}, we use the effective spin polarization $\tilde{\zeta}$ of Eq.~\eqref{eq:def_effspin-0} in order to ensure $S_z$ independence. Thus, $\ntwo^{\text{UEG}}$ will depend indirectly on the on-top pair density of the CASSCF wave function through $\tilde{\zeta}$. When using $\ntwo^{\text{UEG}}(n,\tilde{\zeta})$ in a functional, we will refer to it as ``UEG''. -Another approach to approximate the exact on-top pair density consists in using directly the on-top pair density of the CASSCF wave function. Following the work of some of the present authors, \cite{FerGinTou-JCP-18,GinSceTouLoo-JCP-19} we introduce the extrapolated on-top pair density +The second approach to approximate the exact on-top pair density consists in using directly the on-top pair density of the CASSCF wave function. Following the work of some of the present authors, \cite{FerGinTou-JCP-18,GinSceTouLoo-JCP-19} we introduce the extrapolated on-top pair density \begin{equation} \label{eq:def_n2extrap} \ntwoextrap(\ntwo,\mu) = \qty( 1 + \frac{2}{\sqrt{\pi}\mu} )^{-1} \; \ntwo, \end{equation} which directly follows from the large-$\mu$ extrapolation of the exact on-top pair density derived by Gori-Giorgi and Savin\cite{GorSav-PRA-06} in the context of RSDFT. When using $\ntwoextrap(\ntwo,\mu)$ in a functional, we will simply refer it as ``OT'', which stands for "on-top". -We then define four functionals: +We then define four complementary functionals: \begin{itemize} @@ -621,7 +621,7 @@ We then define four functionals: \bar{E}^\Bas_{\pbeontns} = \int \d\br{} \,\denr \ecmd(\argrpbeontns). \end{equation} \end{itemize} -The performance of each of these four functionals is tested below. +The performance of each of these four functionals is tested in the following. %DFT: BLACK BOX and not CASSCF @@ -642,18 +642,7 @@ The performance of each of these four functionals is tested below. \end{figure*} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{Computational details} -The purpose of the present paper being to investigate the performance of the density-based basis-set correction in regimes of both weak and strong correlation, we study the potential energy curves up to the dissociation limit of a \ce{H10} chain with equally-spaced atoms and the \ce{N2}, \ce{O2}, and \ce{F2} diatomics. For a given basis set, in order to compute the approximate exact ground-state energy in Eq.~\eqref{eq:e0approx}, one needs an approximation to both the FCI energy $\efci$ and the basis-set correction $\efuncbasisFCI$. - -In the case of the \ce{N2}, \ce{O2}, and \ce{F2} molecules for the aug-cc-pVDZ and aug-cc-pVTZ basis sets, approximations to the FCI energies are obtained using frozen-core selected CI calculations followed by the extrapolation scheme proposed by Umrigar and coworkers (see Refs.~\onlinecite{HolUmrSha-JCP-17, SceGarCafLoo-JCTC-18, LooSceBloGarCafJac-JCTC-18, SceBenJacCafLoo-JCP-18, LooBogSceCafJac-JCTC-19, QP2} for more details). All these calculations are performed with the latest version of \textsc{QUANTUM PACKAGE}, \cite{QP2} and will be labeled as exFCI in the following. In the case of \ce{F2}, the correlation energy extrapolation by intrinsic scaling \cite{BytNagGorRue-JCP-07} (CEEIS) are used as reference for the cc-pVXZ (X $=$ D, T, and Q) basis set. The estimated exact potential energy curves are obtained from experimental data \cite{LieCle-JCP-74a} for the \ce{N2} and \ce{O2} molecules, and from extrapolated CEEIS calculations in the case of \ce{F2}. For all geometries and basis sets, the error with respect to the exact FCI energies are estimated to be of the order of $0.5$ mHa. For the \ce{N2}, \ce{O2}, and \ce{F2} molecules, we also performed single-point exFCI calculations in the aug-cc-pVQZ basis set at the equilibrium geometry to obtain reliable estimates of the FCI/CBS dissociation energy. -In the case of the \ce{H10} chain, the approximation to the FCI energies together with the estimated exact potential energy curves are obtained from the data of Ref.~\onlinecite{h10_prx} where the authors performed MRCI+Q calculations with a minimal valence active space as reference (see below for the description of the active space). - -Regarding the complementary density functional, we first perform full-valence CASSCF calculations with the GAMESS-US software\cite{gamess} to obtain the wave function $\psibasis$. Then, all density-related quantities involved in the functional [density $n(\br{})$, spin polarization $\zeta(\br{})$, reduced density gradient $s(\br{})$, and on-top pair density $n_2(\br{})$] together with the local range-separation function $\mu(\br{})$ of Eq.~\eqref{eq:def_mur} are calculated with this full-valence CASSCF wave function. The CASSCF calculations have been performed with the following active spaces: (10e,10o) for \ce{H10}, (10e,8o) for \ce{N2}, (12e,8o) for \ce{O2}, and (14e,8o) for \ce{F2}. - -Also, as the frozen-core approximation is used in all our selected CI calculations, we use the corresponding valence-only complementary density functionals (see Subsec.~\ref{sec:FC}). Therefore, all density-related quantities exclude any contribution from the $1s$ core orbitals, and the range-separation function follows the definition given in Eq.~\eqref{eq:def_mur_val}. - -Regarding the computational cost of the present approach, it should be stressed (see {\SI} for additional details) that the basis set correction represents, for all systems and basis sets studied here, a much smaller computational cost than any of the selected CI calculations. We thus believe that this approach is a significant step towards the routine calculation of near-CBS energetic quantities in strongly correlated systems. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{table*} @@ -707,10 +696,23 @@ Regarding the computational cost of the present approach, it should be stressed \end{ruledtabular} \fnt[1]{From Ref.~\onlinecite{h10_prx}.} \fnt[2]{From the extrapolated valence-only non-relativistic calculations of Ref.~\onlinecite{BytLaiRuedenJCP05}.} -\fnt[3]{CEEIS calculations obtained from non-relativistic calculations of Ref.~\onlinecite{BytNagGorRue-JCP-07}.} +\fnt[3]{From the CEEIS results obtained from non-relativistic calculations of Ref.~\onlinecite{BytNagGorRue-JCP-07}.} \label{tab:extensiv_closed} \end{table*} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsection{Computational details} + +The purpose of the present paper being to investigate the performance of the density-based basis-set correction in regimes of both weak and strong correlation, we study the potential energy curves up to the dissociation limit of a \ce{H10} chain with equally-spaced atoms and the \ce{N2}, \ce{O2}, and \ce{F2} diatomics. For a given basis set, in order to compute the ground-state energy in Eq.~\eqref{eq:e0approx}, one needs an approximation to both the FCI energy $\efci$ and the basis-set correction $\efuncbasisFCI$. + +In the case of the \ce{N2}, \ce{O2}, and \ce{F2} molecules for the aug-cc-pVDZ and aug-cc-pVTZ basis sets, approximations to the FCI energies are obtained using frozen-core selected CI calculations followed by the extrapolation scheme proposed by Holmes \textit{et al.} (see Refs.~\onlinecite{HolUmrSha-JCP-17, SceGarCafLoo-JCTC-18, LooSceBloGarCafJac-JCTC-18, SceBenJacCafLoo-JCP-18, LooBogSceCafJac-JCTC-19, QP2} for more details). All these calculations are performed with the latest version of \textsc{QUANTUM PACKAGE}, \cite{QP2} and will be labeled as exFCI in the following. In the case of \ce{F2}, the correlation energy extrapolated by intrinsic scaling (CEEIS) \cite{BytNagGorRue-JCP-07} are used as reference for the cc-pVXZ (X $=$ D, T, and Q) basis set. The estimated exact potential energy curves are obtained from experimental data \cite{LieCle-JCP-74a} for the \ce{N2} and \ce{O2} molecules, and from extrapolated CEEIS calculations in the case of \ce{F2}. For all geometries and basis sets, the error with respect to the exact FCI energies are estimated to be of the order of $0.5$ mHa. For the \ce{N2}, \ce{O2}, and \ce{F2} molecules, we also performed single-point exFCI calculations in the aug-cc-pVQZ basis set at the equilibrium geometry to obtain reliable estimates of the FCI/CBS dissociation energy. +In the case of the \ce{H10} chain, the approximation to the FCI energies together with the estimated exact potential energy curves are obtained from the data of Ref.~\onlinecite{h10_prx} where the authors performed MRCI+Q calculations with a minimal valence active space as reference (see below for the description of the active space). + +Regarding the complementary functional, we first perform full-valence CASSCF calculations with the GAMESS-US software\cite{gamess} to obtain the wave function $\psibasis$. Then, all density-related quantities involved in the functional [density $n(\br{})$, spin polarization $\zeta(\br{})$, reduced density gradient $s(\br{})$, and on-top pair density $n_2(\br{})$] together with the local range-separation parameter $\mu(\br{})$ of Eq.~\eqref{eq:def_mur} are calculated with this full-valence CASSCF wave function. The CASSCF calculations have been performed with the following active spaces: (10e,10o) for \ce{H10}, (10e,8o) for \ce{N2}, (12e,8o) for \ce{O2}, and (14e,8o) for \ce{F2}. + +Also, as the frozen-core approximation is used in all our selected CI calculations, we use the corresponding valence-only complementary density functionals (see Subsec.~\ref{sec:FC}). Therefore, all density-related quantities exclude any contribution from the $1s$ core orbitals, and the range-separation parameter follows the definition given in Eq.~\eqref{eq:def_mur_val}. + +Regarding the computational cost of the present approach, it should be stressed (see {\SI} for additional details) that the basis-set correction represents, for all systems and basis sets studied here, a much smaller computational cost than any of the selected CI calculations. +%We thus believe that this approach is a significant step towards the routine calculation of near-CBS energetic quantities in strongly correlated systems. \subsection{H$_{10}$ chain} @@ -718,7 +720,7 @@ Regarding the computational cost of the present approach, it should be stressed The study of the \ce{H10} chain with equally distant atoms is a good prototype of strongly-correlated systems as it consists in the simultaneous breaking of 10 covalent $\sigma$ bonds which all interact with each other. Also, being a relatively small system, benchmark calculations at near-CBS values can be obtained (see Ref.~\onlinecite{h10_prx} for a detailed study of this system). We report in Fig.~\ref{fig:H10} the potential energy curves computed using the cc-pVXZ (X $=$ D, T, and Q) basis sets for different levels of approximation. The computation of the atomization energies $D_0$ for each level of theory is reported in Table \ref{tab:d0}. A general trend that can be observed from these data is that, in a given basis set, the quality of the potential energy curves are globally improved by adding the basis-set correction, independently of the approximation level of $\efuncbasis$. Also, no erratic behavior is found when stretching the bonds, which shows that the present procedure (\textit{i.e.} the determination of the range-separation parameter and the definition of the functionals) is robust when reaching the strong-correlation regime. -In other words, smooth potential energy surfaces are obtained with the present basis-set correction. +In other words, smooth potential energy curves are obtained with the present basis-set correction. More quantitatively, the values of $D_0$ are within chemical accuracy (\ie, an error below $1.4$ mHa) from the cc-pVTZ basis set when using the $\pbeontXi$ and $\pbeontns$ functionals, whereas such an accuracy is not even reached at the standard MRCI+Q/cc-pVQZ level of theory. Analyzing more carefully the performance of the different types of approximate density functionals, the results show that $\pbeontXi$ and $\pbeontns$ are very similar (the maximal difference on $D_0$ being 0.3 mHa), and that they give slightly more accurate results than $\pbeuegXi$. These findings provide two important clues on the role of the different physical ingredients used in these functionals: i) the explicit use of the on-top pair density coming from the CASSCF wave function [see Eq.~\eqref{eq:def_n2extrap}] is preferable over the use of the UEG on-top pair density [see Eq.~\eqref{eq:def_n2ueg}] which is somehow understandable, and ii) removing the dependency on any kind of spin polarization does not lead to significant loss of accuracy providing that one employs a qualitatively correct on-top pair density. The latter point is crucial as it shows that the spin polarization in density-functional approximations essentially plays the same role as the on-top pair density. @@ -726,7 +728,6 @@ This could have significant implications for the construction of more robust fam Finally, the reader would have noticed that we did not report the performance of the $\pbeuegns$ functional as its performance are significantly inferior than the three other functionals. The main reason behind this comes from the fact that $\pbeuegns$ has no direct or indirect knowledge of the on-top pair density of the system. Therefore, it yields a non-zero correlation energy for the totally dissociated \ce{H10} chain even if the on-top pair density is vanishingly small. This necessary lowers the value of $D_0$. Therefore, from hereon, we simply discard the $\pbeuegns$ functional. -\subsection{Dissociation of diatomics} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %\begin{figure*} @@ -779,9 +780,11 @@ Finally, the reader would have noticed that we did not report the performance of \end{figure*} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsection{Dissociation of diatomics} + The \ce{N2}, \ce{O2} and \ce{F2} molecules are complementary to the \ce{H10} system for the present study as the level of strong correlation in these diatomics also increases while stretching the bond similarly to the case of \ce{H10}. In addition, these molecules exhibit more important and versatile types of weak correlations due to the larger number of electrons. Indeed, the short-range correlation effects are known to play a strong differential effect on the computation of $D_0$, while the shape of the curve far from the equilibrium geometry is governed by dispersion interactions which are medium to long-range weak-correlation effects. The dispersion forces in \ce{H10} play a minor role in the PES due to the much smaller number of near-neighboring electron pairs compared to \ce{N2}, \ce{O2} or \ce{F2}. Also, \ce{O2} has a triplet ground state and is therefore a good candidate for checking the spin-polarization dependence of the various functionals proposed here. -We report in Figs.~\ref{fig:N2} and \ref{fig:O2} the potential energy curves of \ce{N2} and \ce{O2} computed at various approximation levels using the aug-cc-pVDZ and aug-cc-pVTZ basis sets. Figure \ref{fig:F2} reports the potential energy surface of \ce{F2} using the cc-pVXZ (X $=$ D, T, and Q) basis set. The value of $D_0$ for each level of theory is reported in Table \ref{tab:d0}. +We report in Figs.~\ref{fig:N2} and \ref{fig:O2} the potential energy curves of \ce{N2} and \ce{O2} computed at various approximation levels using the aug-cc-pVDZ and aug-cc-pVTZ basis sets. Figure \ref{fig:F2} reports the potential energy curve of \ce{F2} using the cc-pVXZ (X $=$ D, T, and Q) basis set. The value of $D_0$ for each level of theory is reported in Table \ref{tab:d0}. Just as in \ce{H10}, the quality of $D_0$ is globally improved by adding the basis-set correction and it is remarkable that $\pbeontXi$ and $\pbeontns$ provide again very similar results. The latter observation confirms that the dependency on the on-top pair density allows one to remove the dependency of any kind of spin polarization for a quite wide range of electron density and also for open-shell systems like \ce{O2}. More quantitatively, an error below 1.0 mHa on the estimated exact valence-only $D_0$ is found for \ce{N2}, \ce{O2}, and \ce{F2} with the aug-cc-pVTZ basis set using the $\pbeontns$ functional, whereas such a feat is far from being reached within the same basis set at the near-FCI level. In the case of \ce{F2} it is clear that the addition of diffuse functions in the double- and triple-$\zeta$ basis sets strongly improves the results, a result that can be anticipated due to the strong breathing-orbital effect induced by the ionic valence bond forms in this molecule. \cite{HibHumByrLen-JCP-94} It should be also noticed that when reaching the aug-cc-pVQZ basis set for \ce{N2}, the quality of $D_0$ slightly deteriorates for the $\pbeontXi$ and $\pbeontns$ functionals, but it remains nevertheless more accurate than the estimated FCI $D_0$ and very close to chemical accuracy.