\documentclass[10pt]{article}
\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}
\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
linkcolor = blue,
urlcolor = blue,
citecolor = blue,
anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
\usepackage{natbib}
\renewenvironment{abstract}
{{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
{\bigskip}
\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}
\usepackage{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\usepackage{float}
\iflatexml
\documentclass[%
aps,
jcp,%
amsmath,amssymb,
preprint,%
%longbibliography,
%author-year,%
%author-numerical,%
]{revtex4-1}
\fi
\usepackage{graphicx}
% Include figure files
\usepackage{dcolumn}
% Align table columns on decimal point
\usepackage{bm}
% bold math
\usepackage{braket}
% \usepackage{epstopdf} % Not (yet) supported by LaTeXML.
\newcommand{\mb}[1][]{\mathbf{#1}}
\newcommand{\mr}[1][]{\mathrm{#1}}
\begin{document}
\title{Embedded cluster density approximation for systems with nonlocal
electron correlation}
\author[1]{Chen Huang}%
\affil[1]{Florida State University}%
\vspace{-1em}
\date{\today}
\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup
\selectlanguage{english}
\begin{abstract}
Local correlation methods rely on the assumption that electronic correlation is nearsighted. In this work, we develop a method to alleviate this assumption. The first step is to approximately decompose the electron correlation to the nearsighted and farsighted components based on the wavelength decomposition of electron correlation by Langreth and Perdew. The nearsighted component is then calculated using the recently developed embedded cluster density approximation (ECDA) which is a local correlation method formulated in the context of density functional theory. The farsighted component is calculated based on the system's Kohn-Sham orbitals. The accuracy of this new method depends on the quality of the decomposition.
We examined the method's accuracy by patching the random phase approximation (RPA) correlation energy in a H$_{24}$ chain in which the electron correlation is highly nonlocal. This new method predicts bond stretching energies, RPA correlation potential, and Kohn-Sham eigenvalues in good agreement with the benchmarks. Our results demonstrate the importance of including the farsighted part of electron correlation for studying systems having nonlocal correlations.%
\end{abstract}%
\sloppy
\email{chuang3@fsu.edu}
\section*{Introduction}
Local correlation methods are promising for performing highly accurate electronic structure simulations for large systems.
They are often based on the assumption that electron correlation is nearsighted\hyperref[csl:1]{(Kohn, 1964}; \hyperref[csl:2]{Heine, 1980}; \hyperref[csl:3]{Kohn, 1996}; \hyperref[csl:4]{Prodan \& Kohn, 2005}; \hyperref[csl:5]{Prodan, 2006)}. This assumption is the foundation for developing various embedding methods, such as density-based embedding\hyperref[csl:6]{(Senatore \& Subbaswamy, 1986}; \hyperref[csl:7]{Cortona, 1991}; \hyperref[csl:8]{Wesolowski \& Warshel, 1993}; \hyperref[csl:9]{Govind et al., 1998}; \hyperref[csl:10]{Kohn \& Mattsson, 1998}; \hyperref[csl:11]{Neugebauer et al., 2005}; \hyperref[csl:12]{Jacob \& Visscher, 2006}; \hyperref[csl:13]{Roncero et al., 2008}; \hyperref[csl:14]{Goodpaster et al., 2010}; \hyperref[csl:15]{Huang et al., 2011}; \hyperref[csl:16]{Laricchia et al., 2010}; \hyperref[csl:17]{Fux et al., 2010}; \hyperref[csl:18]{Pavanello \& Neugebauer, 2011}; \hyperref[csl:19]{Fromager, 2015}; \hyperref[csl:20]{Zhu et al., 2016)},
density-matrix-based embedding\hyperref[csl:21]{(Nov{\'a}k et al., 2006}; \hyperref[csl:22]{Knizia \& Chan, 2012}; \hyperref[csl:23]{Manby et al., 2012}; \hyperref[csl:24]{Fornace et al., 2015}; \hyperref[csl:25]{Pernal, 2016}; \hyperref[csl:26]{Welborn et al., 2016}; \hyperref[csl:27]{Yu \& Carter, 2017)}, and self-energy-based embedding\hyperref[csl:28]{(Kananenka et al., 2015}; \hyperref[csl:29]{Chibani et al., 2016)}.
This assumption was also exploited to develop linear-scaling Hartree-Fock method\hyperref[csl:30]{(Schwegler \& Challacombe, 1996}; \hyperref[csl:31]{Ochsenfeld et al., 1998}; \hyperref[csl:32]{Wu et al., 2009)} and low-scaling correlated wave-function methods\hyperref[csl:33]{(Pulay, 1983}; \hyperref[csl:34]{Saebo \& Pulay, 1985}; \hyperref[csl:35]{Kirtman \& Dykstra, 1986}; \hyperref[csl:36]{Stoll, 1992}; \hyperref[csl:37]{Sch{\"u}tz et al., 1999}; \hyperref[csl:38]{Sch{\"u}tz \& Werner, 2000}; \hyperref[csl:22]{Knizia \& Chan, 2012)}. In the context of Kohn-Sham density function theory (KS-DFT)\hyperref[csl:39]{(Hohenberg \& Kohn, 1964}; \hyperref[csl:40]{Kohn \& Sham, 1965)}, this assumption was used to develop low-scaling algorithms for large KS-DFT calculations\hyperref[csl:41]{(Yang, 1991}; \hyperref[csl:42]{Li et al., 1993}; \hyperref[csl:43]{Mauri et al., 1993}; \hyperref[csl:44]{Ordej{\'{o}}n et al., 1993}; \hyperref[csl:45]{Goedecker, 1999)} and random phase approximation (RPA) correlation calculations,\hyperref[csl:46]{(K{\'a}llay, 2015)} and is also the basis for the region-dependent exchange-correlation (XC) energy functional.\hyperref[csl:47]{(Armiento \& Mattsson, 2002)}
Nevertheless, the assumption that electron correlation is nearsighted can fail in practice.
In this work, we develop a method to alleviate this assumption.
The first step is to approximately decomposing electron correlation to the nearsighted and farsighted components.
The accuracy of this new method then depends on the quality of the decomposition. We employ the wavelength-decomposition scheme developed by Langreth and Perdew\hyperref[csl:48]{(Langreth \& Perdew, 1975}; \hyperref[csl:49]{Langreth \& Perdew, 1977)} to decompose electron correlation to the short-wavelength (SW) and long-wavelength (LW) components.
The SW component is considered nearsighted, and the LW component is considered farsighted.
After the decomposition, we treat the SW component using the recently developed embedded cluster density approximation (ECDA)\hyperref[csl:50]{(Huang, 2018}; \hyperref[csl:51]{Huang, 2019)} which is a local correlation method formulated in the framework of KS-DFT.
For the LW component, we calculate it based on the total system.
ECDA shares the similar idea with the local density approximation\hyperref[csl:40]{(Kohn \& Sham, 1965)} (LDA). The major difference between them is that, instead of referencing to a uniform electron gas, in ECDA the XC energy density at a spatial point is calculated based on a cluster's electron density that encloses that point.
This new wavelength-decomposition-based ECDA is termed $\lambda$-ECDA.
In the framework of adiabatic connection fluctuation-dissipation theorem (ACFDT)\hyperref[csl:52]{(Harris \& Jones, 1974}; \hyperref[csl:48]{Langreth \& Perdew, 1975}; \hyperref[csl:53]{Gunnarsson \& Lundqvist, 1976}; \hyperref[csl:49]{Langreth \& Perdew, 1977)}, XC energy is formulated based on the system's linear response functions at different coupling constants. Let's consider a cluster embedded in a system. Since the linear response for a closed system is not nearsighted in general,\hyperref[csl:54]{(Fias et al., 2017)} inside the cluster, the system and the cluster's linear response functions are not expected to be very similar. Therefore, inside the cluster, the cluster's XC energy density may not be a good approximation to the system's XC energy density. This is the reason for ECDA's poor performance in systems having highly nonlocal correlation, as observed in this work.
The assumption underlying $\lambda$-ECDA is that the \textit{short-wavelength} components of the linear response functions (at different coupling constants) of the system and the cluster are similar inside the cluster. This assumption is motivated by previous observations that LDA\hyperref[csl:40]{(Kohn \& Sham, 1965)} was excellent at short wavelength for certain cases.\hyperref[csl:49]{(Langreth \& Perdew, 1977}; \hyperref[csl:55]{Langreth \& Perdew, 1980}; \hyperref[csl:56]{Burke et al., 1994)} However, we note that this short-wavelength assumption is not true in general.\hyperref[csl:56]{(Burke et al., 1994)}.
The paper is organized as follows. First, we give a description of ECDA and discuss the wavelength decomposition of the RPA correlation energy\hyperref[csl:57]{(Bohm \& Pines, 1951}; \hyperref[csl:48]{Langreth \& Perdew, 1975}; \hyperref[csl:49]{Langreth \& Perdew, 1977)} which is the high-level correlation energy to be patched in this work. The performance of $\lambda$-ECDA is then investigated on a H$_{24}$ chain. We demonstrate the importance of the LW RPA correlation for treating this system in which the electron correlation is highly nonlocal. An important parameter in $\lambda$-ECDA is the cutoff wavelength for the wavelength decomposition, we develop a practical method to determine it for given cluster sizes.
\section*{Theoretical methods}
\subsection*{The embedded cluster density approximation}
The $\lambda$-ECDA method developed in this work is based on the ECDA method formulated in Ref.\hyperref[csl:51]{(Huang, 2019)}. The major difference between them is that ECDA directly patches the clusters' XC energies, while $\lambda$-ECDA patches the SW components of the clusters' XC energies, with the LW component of the XC energy calculated based on the system's KS orbitals.
Therefore, ECDA calculations depend on the cluster sizes, while $\lambda$-ECDA calculations depend on both the cluster sizes and the cutoff wavelength used in the wavelength decomposition.
In ECDA, a system's total energy is defined as
\begin{equation}
\label{eq:etot}
E_{tot} = T_s + J + E_{ext} + E_{xc} - TS + E_p
\end{equation}
in which $T_s= - \sum_j f_j \int d\mb r' \psi_j^*(\mb r')\frac{1}{2}\nabla^2\psi_j(\mb r')$ is the system's KS kinetic energy with \{$\psi_j$\} and \{$f_j$\} being the system's KS orbitals and occupation numbers. In this work, the occupation numbers are assigned based on the Fermi-Dirac statistics.
$J=\frac{1}{2}\iint d\mb r d\mb r'\rho(\mb r)\rho(\mb r')/|\mb r - \mb r'|$ is the Hartree energy with $\rho(\mb r)$ being the system's electron density.
$E_{ext} = \int d\mb r \rho(\mb r) v_{ext}(\mb r)$ is the external energy and $v_{ext}(\mb r)$ is the system's external potential.
$T$ is the smearing temperature, and $S$ is the electronic entropy. $E_p$ is the penalty energy for regularizing the system's XC potential obtained from solving the optimized effective potential (OEP) equation.\hyperref[csl:58]{(Sharp \& Horton, 1953}; \hyperref[csl:59]{Talman \& Shadwick, 1976}; \hyperref[csl:60]{Sahni et al., 1982}; \hyperref[csl:61]{K{\"u}mmel \& Kronik, 2008)} In this work, the penalty energy is set to zero, since the system's OEP equation is solved exactly without employing regularization.
In ECDA, the system's XC energy $E_{xc}$ is calculated through patching the clusters' XC energy densities
\begin{equation}
\label{eq:xcpatch}
E_{xc} = \sum_{i=1}^{N_{\mr{atom}}} \int d\mb r w_i(\mb r) \varepsilon_{xc}^{clu,i}(\mb r).
\end{equation}
where $w_i(\mb r)$ is Becke's weight\hyperref[csl:62]{(Becke, 1993)} function defined on atom $i$. $N_{\mr{atom}}$ is the number of atoms. The sum is over all the atoms, since each atom has a cluster. In Ref.\hyperref[csl:51]{(Huang, 2019)}, the cluster's XC energy density, $\varepsilon_{xc}^{clu,i}(\mb r)$, is defined as\hyperref[csl:50]{(Huang, 2018)}
\begin{equation}
\label{eq:xc_eq}
\varepsilon_{xc}^{clu,i}(\mb r) =
\varepsilon_{xc,high}^{clu,i}(\mb r)
- \varepsilon_{xc,low}^{clu,i}(\mb r)
+ \varepsilon_{xc,low}^{tot}(\mb r),
\end{equation}
where $\varepsilon_{xc,low}^{clu,i}(\mb r) $ and $\varepsilon_{xc,low}^{tot}(\mb r)$ are the cluster and the system's XC energy densities calculated using a low-level XC energy functional, such as LDA.
The last two terms in Eq.~\ref{eq:xc_eq} represent the correction from the environment.
Different from the previous work, in this work we only patch the correlation energy, and $\varepsilon_{c}^{clu,i}(\mb r)$ is just set to the cluster's RPA correlation energy density without considering the correction from the environment, that is,
\begin{equation}
\varepsilon_{c}^{clu,i}(\mb r)=\varepsilon_{c,\mr{RPA}}^{clu,i}(\mb r).
\end{equation}
The system's exact exchange (EXX) energy is calculated based on the system's KS orbitals.
A key step in ECDA is to define the clusters' electron densities, based on which their XC energy densities are calculated. A cluster's electron density is defined by partitioning the system's electron density among it and its environment, such that,
\begin{equation}
\label{eq:match}
\rho_{tot}(\mb r) = \rho_{clu}(\mb r) + \rho_{env}(\mb r).
\end{equation}
To make Eq.~\ref{eq:match} to hold, an embedding potential\hyperref[csl:63]{(Cohen \& Wasserman, 2006}; \hyperref[csl:15]{Huang et al., 2011)} is attached to both the cluster and its environment. The environment then influences the cluster through the embedding potential.
By adjusting the embedding potential, Eq.~\ref{eq:match} holds.
In this work, the partitioning is performed using a modified Wu-Yang method as described in Appendix A.
The system's XC potential $v_{xc}(\mb r)$ is needed to perform ECDA calculations self-consistently. $v_{xc}$ is calculated by solving the OEP equation
\begin{equation}
\frac{\delta E_{xc}}{\delta v_s(\mb r)} =
\int d\mb r' v_{xc}(\mb r')
\chi_0(\mb r, \mb r'),
\end{equation}
where $\chi_0$ is the system's KS linear response function and $v_s(\mb r)$ is the system's KS potential. We then need $\delta E_{xc}/\delta v_s(\mb r)$.
Since $E_{xc}$ is obtained through patching clusters' XC energy densities, it depends on $v_s(\mb r)$ in two ways. First, keeping the embedding potential fixed, a cluster's electron density changes as $v_s(\mb r)$ changes. Second, the embedding potential also changes as $v_s(\mb r)$ changes. Detailed procedure for calculating $\delta E_{xc}/\delta v_s(\mb r)$ was given in Ref.\hyperref[csl:51]{(Huang, 2019)}.
In this work, we focus on patching RPA correlation energy, and the EXX energy is calculated based on the system's KS orbitals. In the context of ACFDT, EXX energy is defined as
\begin{equation}
\label{eq:exx}
E_{x}^{\mathrm{EXX}} =
- \sum_l\sum_{k} f_k c_{kl} \iint d \mb r d\mb r' \frac{\psi_k(\mb r)\psi_l(\mb r)\psi_l(\mb r')\psi_k(\mb r')}{|\mb r - \mb r'|},
\end{equation}
where $k$ and $l$ loop over all the orbitals, $c_{kl}=1+\mathrm{sgn}(e_k -e_l)$, and $\{e_k\}$ are the KS eigenvalues.\hyperref[csl:64]{(Harl et al., 2010)}
The RPA correlation energy is
\begin{eqnarray}
\label{eq1}
E_{c}^{\mathrm{RPA}} = && -\frac{1}{2\pi}
\int_0^\infty d\omega
\int_0^1 d\eta \iint d\mb r d \mb r'
v(\mathbf r, \mathbf r') \nonumber \\
&& \times [\chi_{\mr{RPA}}^\eta(\mb r,\mb r';i\omega) - \chi_0(\mb r, \mb r'; i\omega)],
\end{eqnarray}
where $\eta$ is the switching parameter to connect the non-interacting and fully interacting systems. $\chi_{\mr{RPA}}^\eta(\mb r,\mb r';i\omega)$ is the RPA linear response function for $\eta v(\mb r, \mb r')$ at the imaginary frequency $i\omega$.
$\chi_0(\mb r,\mb r';i\omega)$ is the KS linear response at the imaginary frequency $i\omega$. $\chi_{\mr{RPA}}^\eta(\mb r,\mb r';i\omega)$ depends on $\chi_0(\mb r,\mb r';i\omega)$ as
\begin{eqnarray}
\label{eq2}
\chi_{\mathrm{RPA}}^\eta(\mb r,\mb r';i\omega) = &&
\chi_0(\mb r,\mb r';i\omega)
+ \iint d \mb r_1 d \mb r_2
\chi_0(\mb r, \mb r_1; i\omega)\nonumber \\
&& \times \eta v(\mb r_1,\mb r_2)\chi_{\mathrm{RPA}}^\eta(\mb r_2, \mb r'; i\omega).
\end{eqnarray}
For the patching, we define the RPA correlation energy density as
\begin{eqnarray}
\label{eq:rpa_density}
\varepsilon_c^{\mr{RPA}}(\mb r) = &&
-\frac{1}{2\pi}\int_0^\infty d\omega
\int_0^1 d\eta
\int d\mb r' v(\mb r, \mb r') \nonumber \\
&& \times
[\chi_{\mr{RPA}}^\eta(\mb r, \mb r'; i\omega)
- \chi_0(\mb r, \mb r'; i\omega)].
\end{eqnarray}
\subsection*{Wavelength decomposition of the RPA correlation energy}
In $\lambda$-ECDA, only the SW component of the RPA correlation energy is patched. To decompose RPA correlation energy, we write it in the Fourier space\hyperref[csl:49]{(Langreth \& Perdew, 1977)}
\begin{eqnarray}
E_c^{\mr{RPA}} = &&
\int \frac{d^3 k}{(2\pi)^3} \frac{4\pi}{k^2}
\int_0^\infty d \omega \int_0^1 d \eta \nonumber \\
&& \times
[\tilde \chi^\eta_{\mr{RPA}}(\mb k, -\mb k; i\omega) -
\tilde \chi_0(\mb k, -\mb k; i\omega)]
\end{eqnarray}
where $4\pi/k^2$ is the Fourier transformation of the Coulomb potential, and
\begin{eqnarray}
&& \tilde \chi^\eta_{\mr{RPA}}(\mb k, -\mb k'; i\omega) = \iint d\mb r d\mb r' e^{i\mb k \mb r} e^{-i \mb k' \mb r'}\chi^\eta_{\mr{RPA}}(\mb r, \mb r'; i\omega) \\
&& \tilde \chi_0(\mb k, -\mb k'; i\omega) = \iint d\mb r d\mb r' e^{i\mb k \mb r} e^{-i \mb k' \mb r'}\chi_0(\mb r, \mb r'; i\omega).
\end{eqnarray}
To decompose $E_c^{\mr{RPA}}$ to the SW and LW components, $\tilde \chi_{\mathrm{RPA}}^\eta$ and $\tilde \chi_0$ are decomposed to the SW and LW components. This is equivalent to decompose the Coulomb potential into the SW and LW components, which is the approach we take here
\begin{equation}
\tilde {v}(k) = \tilde {v}_{\mr{sw}}(k) + \tilde {v}_{\mr{lw}} (k),
\end{equation}
where $\tilde{v}$ is the Fourier transformation of $v$. $\tilde{v}_{\mr{sw}}$ and $\tilde{v}_{\mr{lw}}$ are the SW and LW components defined using a Gaussian function
\begin{eqnarray}
\label{eq:vcsw}
&& \tilde {v}_{\mr{sw}}(k) = (1-e^{-(k/k_c)^2})\tilde {v}(k) \\
\label{eq:vclw}
&& \tilde {v}_{\mr{lw}}(k) = e^{-(k/k_c)^2}\tilde {v}(k).
\end{eqnarray}
$k_c$ is the cutoff wave-vector and the cutoff wavelength is
\begin{equation}
\lambda_c = \frac{2\pi}{k_c}.
\end{equation}
The advantage of using a Gaussian function to partition $\tilde{v}(k)$ is that for three-dimensional systems this is equivalent to partitioning the Coulomb potential into short-range and long-range components in the real space using the error function. Then we achieve the wavelength decomposition and the range separation at the same time. But, this is not the case for one-dimensional systems studied in this work, as discussed later.
The SW and LW components of the RPA correlation energy are defined by replacing $v_c$ in Eq.~\ref{eq1} with $v_{\mr{sw}}$ and $v_{\mr{lw}}$, respectively,
\begin{eqnarray}
E_{\mr{c,sw}}^{\mathrm{RPA}} = && -\frac{1}{2\pi}
\int_0^\infty d\omega
\int_0^1 d\eta \iint d\mb r d \mb r'
v_{\mr{sw}}(\mathbf r, \mathbf r') \nonumber \\
&& \times [\chi_{\mr{RPA}}^\eta(\mb r,\mb r';i\omega) - \chi_0(\mb r, \mb r'; i\omega)],
\end{eqnarray}
and
\begin{eqnarray}
E_{\mr{c,lw}}^{\mathrm{RPA}} = && -\frac{1}{2\pi}
\int_0^\infty d\omega
\int_0^1 d\eta \iint d\mb r d \mb r'
v_{\mr{lw}}(\mathbf r, \mathbf r') \nonumber \\
&& \times [\chi_{\mr{RPA}}^\eta(\mb r,\mb r';i\omega) - \chi_0(\mb r, \mb r'; i\omega)].
\end{eqnarray}
The energy density of the SW RPA correlation is defined by replacing $v$ in Eq.~\ref{eq:rpa_density} with $v_\mr{sw}$ as
\begin{eqnarray}
\label{eq:rpa_density_sw}
\varepsilon_{\mr{c,sw}}^{\mr{RPA}}(\mb r) = &&
-\frac{1}{2\pi}\int_0^\infty d\omega
\int_0^1 d\eta
\int d\mb r' v_{\mr{sw}}(\mb r, \mb r') \nonumber \\
&& \times
[\chi_{\mr{RPA}}^\eta(\mb r, \mb r'; i\omega)
- \chi_0(\mb r, \mb r'; i\omega)].
\end{eqnarray}
\subsection*{Wavelength decomposition of the Coulomb potential for one-dimensional systems}
In this work, the performance of $\lambda$-ECDA is examined in one-dimensional systems. To avoid the non-integrable singularity of Coulomb potential at $r=0$, we employ a softened Coulomb potential to describe electron-electron interaction
\begin{equation}
v(r) = \frac{1}{\sqrt{1+r^2}}.
\end{equation}
We first transform it to the Fourier space as
\begin{equation}
\tilde v(k) = \int_{-\infty}^{\infty}
dr \frac{\cos(kr)}{\sqrt{1+r^2}} =
2 K_0(|k|)
\end{equation}
where $K_0(r)$ is the modified Bessel function of the second kind. Due to using a real-space KS-DFT program, we work with the SW and LW components of $v(r)$ in the real space. The SW and LW parts of the Coulomb potential in real space are obtained as
\begin{eqnarray}
\label{eq:sw1d}
&& v_{\mr{sw}}(r) =
\int_{-\infty}^{\infty} \frac{dk}{2\pi} \cos(-kr) \tilde v(k) [1-e^{-(k/k_c)^2}] \\
\label{eq:lw1d}
&& v_{\mr{lw}}(r) = \int_{-\infty}^{\infty} \frac{dk}{2\pi} \cos(-kr) \tilde v(k) e^{-(k/k_c)^2}.
\end{eqnarray}
$v_{\mr{sw}}(r)$ and $v_{\mr{lw}}(r)$ are then used in our $\lambda$-ECDA calculations.
The decomposition of the Coulomb potential for $\lambda_c=6$ Bohr is illustrated in Fig.~\ref{fig:v_srlr}. We see that for one-dimensional systems there is no clear range separation in the real space. $v_{\mr{lw}}(r)$ is large near $r=0$.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/v-SR-LR/v-SR-LR}
\caption{{\label{fig:v_srlr} The SW and LW components of the Coulomb potential in (a) the real space and (b) the Fourier space for $\lambda_c=6$ Bohr.%
}}
\end{center}
\end{figure}
\section {Numerical Details}
All KS-DFT and ECDA calculations are performed using a one-dimensional, real-space KS-DFT code written in the MATLAB program.\hyperref[csl:65]{(2019)} An evenly-spaced H$_{24}$ chain is used for the calculations. A softened Coulomb potential $1/\sqrt{1+r^2}$ is used to describe the electron-electron, ion-electron, and ion-ion interactions. The nuclear charge of the pseudo-hydrogen is 1.2$e$ in order to have a strong ion-electron interaction. Each pseudo-hydrogen contains one electron.
A uniform grid is used to expand KS orbitals.
A grid spacing of 0.33 Bohr is used for all calculations and converges KS-DFT-LDA energy to 0.5 mHa.
The KS equation is solved by first replacing the Laplacian operator with five-point stencil. A smearing temperature of 0.05 eV is used for all calculations. This converts the KS equation to an eigenvalue problem which is solved using MATLAB.
Local density approximations for exchange and correlation energies for one-dimensional systems are taken from Refs.\hyperref[csl:66]{(Wagner et al., 2012}; \hyperref[csl:67]{Helbig et al., 2011)}.
EXX and RPA potentials are calculated by solving the OEP equation. For RPA correlation energy, its derivative with respect to the system's KS potential is calculated using the centered finite difference method by adding/subtracting 0.01 a.u. to/from the KS potential at each grid point. For EXX energy, its derivative with respect to system's KS potential is calculated analytically.\hyperref[csl:50]{(Huang, 2018)}
For RPA correlation, the integral over the frequency is performed using the Gauss-Legendre method with 15 points on a re-scaled mesh.\hyperref[csl:68]{(Olsen \& Thygesen, 2013)} The maximum frequency for the mesh is 15 a.u. This mesh converges H$_{24}$'s RPA correlation energy to 0.1 mHa. The integration over the coupling parameter $\eta$ is calculated using the Gauss-Legendre method on a regular mesh with 8 points. This converges the RPA correlation energy to $1\times10^{-6}$ Ha.
Becke's weight function used in this work is softer than that in Ref.\hyperref[csl:62]{(Becke, 1993)} by setting $k$ to 1 in the Eq.19 in Ref.\hyperref[csl:62]{(Becke, 1993)}.
\section*{Results and Discussion}
\subsection*{The importance of long-wavelength component of RPA correlation}
We investigate the performance of $\lambda$-ECDA using a uniformly spaced H$_{24}$ chain. The H-H bond length is 2.0 Bohr. In this system, the electron correlation is highly nonlocal, providing a good testbed for $\lambda$-ECDA.
The SW and LW components of the RPA correlation energy for different $\lambda_c$ are plotted in Fig.~\ref{fig:sr_lr}. The sum of them is the system's RPA energy (-0.7842 Hartree) which is marked by the dashed line. RPA correlation energies are calculated based on the LDA orbitals.
Fig.~\ref{fig:sr_lr} shows that the LW RPA correlation energy decays slowly at large $\lambda_c$, indicating that the LW linear response still contributes to the RPA correlation energy at large $\lambda_c$.
The LW RPA correlation energy is also important for the system's RPA correlation potential. In Fig.~\ref{fig:vrpa}, we plot the RPA correlation potentials due to the SW and LW RPA correlation energies. The sum of them is equal to the system's RPA correlation potential marked by the solid black curves.
The contribution from the LW component is significant even at large $\lambda_c$.
At $\lambda_c=30$ Bohr, the LW RPA correlation energy is 2.7 mHa per atom, which is 11\% of the system's RPA correlation energy. Fig.~\ref{fig:vrpa}(a) shows that the LW RPA correlation potential dominates inside the H$_{24}$ chain and only becomes weak at the ends of the chain.
At $\lambda_c=60$ Bohr (Fig.~\ref{fig:vrpa}(b)), the LW RPA correlation energy is 0.8 mHa per atom which is 3.5\% of total RPA correlation energy, but the LW RPA correlation potential is still comparable with the SW RPA correlation potential inside the chain. These results indicate that the LW RPA correlation cannot be ignored for systems having highly nonlocal electron correlation, even if its contribution to the total correlation energy is not very significant.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/H24-sr-lr-wvec/H24-sr-lr-wvec}
\caption{{\label{fig:sr_lr} The SW and LW components of the RPA correlation energies per atom of H$_{24}$ calculated at different cutoff wavelengths. The total RPA correlation energy per atom is denoted by the dashed line.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/H24-sr-lr-vrpa-lambda-30-60/H24-sr-lr-vrpa-lambda-30-60}
\caption{{\label{fig:vrpa} The SW and LW components of the RPA correlation potentials calculated at $\lambda=30$ Bohr and $\lambda=60$ Bohr.
The system's RPA correlation potential (denoted by ``total'') is shown for comparison. The potentials are calculated based on the LDA orbitals. Hydrogen atoms are marked by the black dots.%
}}
\end{center}
\end{figure}
\subsection*{Determination of the maximum cutoff wavelength}
For a given cluster sizes, we need $\lambda_c$ to be as large as possible, in order to have more correlation energy from the cluster calculations. On the other hand, we need to keep $\lambda_c$ small enough such that the SW correlation energy can be captured by the clusters.
In what follows, we develop a method to determine the maximum cutoff wavelength, denoted as $\lambda_c^{\mr{max}}$, for given cluster sizes. Here, the cluster sizes are denoted by $N_b$, which means that the atoms up to the $N_b^{th}$ nearest neighbors are used to define each cluster.
To find $\lambda_c^{\mr{max}}$, we define the ratio
\begin{equation}
\label{eq:R1}
R(\lambda_c; N_b)=\frac{E_{\mr{c,sw}}^{\mr{RPA,patched}}(\lambda_c;N_b)}
{E_{\mr{c,sw}}^{\mr{RPA}}(\lambda_c)},
\end{equation}
which measures the difference between the patched SW RPA correlation energy ($E_{\mr{c,sw}}^{\mr{RPA,patched}}$) and the system's SW RPA correlation energy ($E_{\mr{c,sw}}^{\mr{RPA}}$). For given $N_b$, $R$ is a function of $\lambda_c$.
Fig.~\ref{fig:ratio} shows $R$ for several $N_b$ values. One common feature of these curves is that, at small $\lambda_c$, they are relatively flat and are not equal to one, which means that the relative difference between $E_{\mr{c,sw}}^{\mr{RPA,patched}}$ and $E_{\mr{c,sw}}^{\mr{RPA}}$ is persistent even at small $\lambda_c$. This indicates that the SW component of the RPA correlation is not fully nearsighted at short wavelength and contains a small farsighted component that the clusters cannot fully capture.
Fortunately, this farsighted component only contributes to a small fraction of the system's RPA energy. In the following sections, we find that the failure to fully capture this farsighted component does not cause significant errors.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/H24-lambda/H24-lambda}
\caption{{\label{fig:ratio} The log-log plot of $R(\lambda_c; N_b)$ calculated for different $\lambda_c$ and $N_b$. The limiting $R$ values at $\lambda_c=\infty$ are marked by the dashed lines.%
}}
\end{center}
\end{figure}
In Fig.~\ref{fig:ratio}, another common feature of the curves is that they quickly drop to a plateau beyond certain $\lambda_c$. The reason for the drop is that a cluster cannot capture the SW RPA correlation energy with a wavelength much larger than its size.
Therefore, $\lambda_c^{\mr{max}}$ is determined as the point where $R$ starts to drop.
To estimate that point, we calculate $R$'s curvature for a given $N_b$
\begin{equation}
\kappa = \frac{R''}{(1+ R'^2)^{3/2}}
\end{equation}
where $R' = \partial R / \partial \lambda_c$ and $R''=\partial^2 R / \partial \lambda_c^2$.
The turning point is the one that has the smallest curvature. Taking $N_b=2$ for example, Fig.~\ref{fig:curv}(a) plots $R$ versus $\lambda_c$, whose curvature is shown in Fig.~\ref{fig:curv}(b).
$\lambda_c^{\mr{max}}$ gives the smallest curvature and is marked by the red crosses in both subplots.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/H24-curvature/H24-curvature}
\caption{{\label{fig:curv} (a) The log-log plot of $R$ versus $\lambda_c$ for $N_b=2$. (b) The curvature of $R(\lambda_c;N_b=2)$. $\lambda_c^{\mr{max}}$ are marked by the red crosses in both subplots.%
}}
\end{center}
\end{figure}
Using this method, we determine $\lambda_c^{\mr{max}}$ for different cluster sizes. Fig.~\ref{fig:max_lambda} shows that $\lambda_c^{\mr{max}}$ is generally proportional to clusters' sizes.
In Table~\ref{tb:maxlambda}, we compare the system's and the patched SW RPA correlation energies for different $N_b$, calculated using $\lambda_c^{\mr{max}}$. They are calculated based on LDA orbitals. $E_{\mr{c,sw}}^{\mr{RPA,patched}}$ agrees with $E_{\mr{c,sw}}^{\mr{RPA}}$ within 1 mHa per atom for all $N_b$ values, which is expected since $\lambda_c^{\mr{max}}$ is the point where $R$ starts to drop.
Table~\ref{tb:maxlambda} also shows that, as $\lambda_c^{\mr{max}}$ increases, the LW component of RPA correlation energy ($E_{\mr{c,lw}}^{\mr{RPA}}$) decreases.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/H24-opt-lambda/H24-opt-lambda}
\caption{{\label{fig:max_lambda} $\lambda_c^{\mr{max}}$ determined for different $N_b$ values.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{table}[]
\centering
\caption{{$\lambda_c^{\mr{max}}$ (in Bohr) determined for different $N_b$ and the corresponding SW and LW RPA correlation energies.
The patched SW RPA correlation energies are listed to assess the quality of patching. All energies are in Hartree and are calculated based on the LDA orbitals.}}
\label{tb:maxlambda}
\begin{tabular}{ccccc}
\hline\hline
$N_b$ & $\lambda_c^{\mr{max}}$
& $E_{\mr{c,sw}}^{\mr{RPA}}$ &
$E_{\mr{c,sw}}^{\mr{RPA,patched}}$ &
$E_{\mr{c,lw}}^{\mr{RPA}}$ \\ \hline
1 & 8.615 & -0.0126 & -0.0124 & -0.0118 \\
2 & 11.869 & -0.0157 & -0.0152 & -0.0088 \\
3 & 18.931 & -0.0193 & -0.0195 & -0.0051 \\
4 & 22.629 & -0.0204 & -0.0199 & -0.0041 \\
5 & 31.103 & -0.0219 & -0.0220 & -0.0026 \\
6 & 33.802 & -0.0222 & -0.0217 & -0.0023 \\
7 & 37.669 & -0.0225 & -0.0227 & -0.0019 \\ \toprule
\end{tabular}
\end{table}
Despite of its simplicity, the above method for determining $\lambda_c^{\mr{max}}$ is computationally infeasible in practice, since the cost for computing $E_{\mr{c,sw}}^{\mr{RPA}}$ will be high. To fix this issue, we define a new ratio
\begin{equation}
\label{eq:Rbar}
\bar R(\lambda_c;N_b,N_b')
=\frac{E_{\mr{c,sw}}^{\mr{RPA,patched}}(\lambda_c;N_b)}
{E_{\mr{c,sw}}^{\mr{RPA,patched}}(\lambda_c;N_b')},
\end{equation}
which relies on $E_{\mr{c,sw}}^{\mr{RPA,patched}}$
evaluated for two different cluster sizes $N_b$ and $N_b'$ with $N_b