Second-order perturbative correlation energy functional in the ensemble density-functional theory

We use atomic units [$e=\hbar=m_{e}=1/(4\pi\epsilon_{0})=1$] unless otherwise specified. The interacting system is described by the Schödinger equation:

$$\hat{H}\left|\Psi_{ik}\right>=\mathcal{E}_{i}\left|\Psi_{ik}\right>,$$

(1)

where the Hamiltonian $\hat{H}$ is $\hat{H}=\hat{T}+\hat{V}_{\rm ext}+\hat{V}_{\rm ee}$ with $\hat{T}$, $\hat{V}_{\rm ext}$ and $\hat{V}_{\rm ee}$ being the kinetic, external potential and electron-electron interaction potential operators, and $i$ denotes a set of degenerate states (‘multiplet’) and $k$ an individual state in the multiplet. The multiplets are ordered by energy with $i=0$ being the ground state. The ordering of wavefunctions within a multiplet is arbitrary but fixed. An ensemble in EDFT contains consecutive multiplets from the ground state to a highest-energy multiplet $I$. Each state in multiplet $i$ is assigned a weight $w_{i}$ satisfying $w_{i}\geq w_{j}$ for $i<j$.GOKb88 We require $\sum_{i=0}^{I}g_{i}w_{i}=1$ for simplicity, where $g_{i}$ is the degeneracy of multiplet $i$. We denote all weights as ${\{w\}}$ in the following.

The ensemble density and energy are defined as

$$n_{\{w\}}({\bf r})=\operatorname{tr}\{\hat{D}_{\{w\}}\hat{n}({\bf r})\}=\sum_{i=0}^{I}w_{i}\sum_{k=1}^{g_{i}}n_{ik}({\bf r}),$$

(2)

and

$$E_{\{w\}}=\operatorname{tr}\{\hat{D}_{\{w\}}\hat{H}\}=\sum_{i=0}^{I}w_{i}\sum_{k=1}^{g_{i}}\mathcal{E}_{i},$$

(3)

where the ensemble density matrix $\hat{D}_{\{w\}}$ is

$$\hat{D}_{\{w\}}=\sum_{i=0}^{I}w_{i}\sum_{k=1}^{g_{i}}\left|\Psi_{ik}\right>\left<\Psi_{ik}\right|,$$

(4)

and $n_{ik}({\bf r})=\left\langle\Psi_{ik}\left|\hat{n}({\bf r})\right|\Psi_{ik}\right\rangle$.

EDFT proves the one-to-one correspondence betweeen $n_{\{w\}}$ and the external potential given $\hat{V}_{\rm ee}$GOKb88 , so $E_{\{w\}}$ is a functional of $n_{\{w\}}$. One can define a non-interacting EKS system with the same $n_{\{w\}}$ as the interacting system, with the EKS orbitals satisfying

$$\left\{-\frac{1}{2}\nabla^{2}_{\bf r}+v_{\text{s},{\{w\}}}[n_{\{w\}}]({\bf r})\right\}\phi_{\mu,{\{w\}}}({\bf r})=\epsilon_{\mu,{\{w\}}}\phi_{\mu,{\{w\}}}({\bf r}).$$

(5)

The EKS density matrix $\hat{D}_{\text{s},{\{w\}}}$ is

$$\hat{D}_{\text{s},{\{w\}}}=\sum_{i=0}^{I}w_{i}\sum_{k=1}^{g_{i}}\left|\Phi_{ik,{\{w\}}}\right>\left<\Phi_{ik,{\{w\}}}\right|,$$

(6)

where $\left|\Phi_{ik,{\{w\}}}\right>$ is an EKS wavefunction corresponding to the $\left|\Psi_{ik}\right>$ state of the interacting system. The existence of an adiabatic connectionPY89 ; FNM03 between the interacting and EKS systems is assumed. $n^{\text{KS}}_{\{w\}}({\bf r})=\operatorname{tr}\{\hat{D}_{\text{s},{\{w\}}}\hat{n}({\bf r})\}=n_{\{w\}}({\bf r})$ if $v_{\text{s},{\{w\}}}$ is exact.

Unlike in ground-state DFT, $\left|\Phi_{ik,{\{w\}}}\right>$ is not necessarily a single Slater determinant of EKS orbitals since the EKS degenercies are in general greater than or equal to the corresponding interacting ones.PYTB14 We require $\Phi_{ik,{\{w\}}}$ to have the same spatial and spin symmetries as $\Psi_{ik}$PYTB14 to distinguish the otherwise degenerate EKS states by linearly combine the EKS Slater determinants:

$$\left|\Phi_{ik,{\{w\}}}\right>=\sum_{p=1}^{\tilde{g}_{\tilde{i}}}C_{ikp}\left|\tilde{\Phi}_{\tilde{i}p,{\{w\}}}\right>,$$

(7)

where $\tilde{i}$ and $\tilde{g}_{\tilde{i}}$ denote the EKS multiplet and its degeneracy corresponding to interacting multiplet $i$, $\tilde{\Phi}_{\tilde{i}p,w}$ denotes the EKS Slater determinant $p$ in EKS multiplet $\tilde{i}$, and $C_{ikp}$ is the linear combination coefficient.

$E_{\{w\}}$ can be decomposed as

$$\begin{split}E_{\{w\}}[n]=&T_{\text{s},{\{w\}}}[n]+V_{\rm ext}[n]+E_{\text{Hx},{\{w\}}}[n]+E_{\text{c},{\{w\}}}[n]\\ =&\operatorname{tr}\{\hat{D}_{\text{s},{\{w\}}}\hat{T}\}+\int\mathrm{d}^{3}r\,n({\bf r})v_{\rm ext}({\bf r})\\ &+\operatorname{tr}\{\hat{D}_{\text{s},{\{w\}}}\hat{V}_{\rm ee}\}+E_{\text{c},{\{w\}}}[n].\end{split}$$

(8)

The ensemble Hartree-exchange energy is

$$E_{\text{Hx},{\{w\}}}=\operatorname{tr}\{\hat{D}_{\text{s},{\{w\}}}{\hat{V}_{\rm ee}}\}=\sum_{i=0}^{I}w_{i}\sum_{k=1}^{g_{i}}\left\langle\Phi_{ik,{\{w\}}}\left|\hat{V}_{\rm ee}\right|\Phi_{ik,{\{w\}}}\right\rangle,$$

(9)

which can be evaluated using the Slater-Condon rulesL14 . The ensemble correlation energy $E_{\text{c},{\{w\}}}$ is defined by Eq. (8) and must be approximated in practice. The EKS potential $v_{\text{s},{\{w\}}}({\bf r})$ is determined variationally by $\delta E_{\{w\}}[n]/\delta n({\bf r})=0$ as

$$\begin{split}v_{\text{s},{\{w\}}}[n]({\bf r})&=v_{\rm ext}({\bf r})+v_{\text{Hx},{\{w\}}}[n]({\bf r})+v_{\text{c},{\{w\}}}[n]({\bf r})\\ &=v_{\rm ext}({\bf r})+\frac{\delta E_{\text{Hx},{\{w\}}}[n]}{\delta n({\bf r})}+\frac{\delta E_{\text{c},{\{w\}}}[n]}{\delta n({\bf r})}.\end{split}$$

(10)

EDFT calculates excitation energies either by subtracting ensemble energies or by taking derivative of $E_{\{w\}}$ with respect to the weightsGOK88 ; DF19 , and both requires solving Eq. (5) self-consistently. The DEC method of our previous workYPBU17 can be used to calculate excitation energies without extra self-consistent calculations other than the ground-state one. The DEC method uses a special type of ensemble defined as

$$w_{i}(\mathtt{w})=\left\{\begin{array}[]{cl}\frac{1-\mathtt{w}(M_{I}-g_{0})}{g_{0}}&i=0,\\ \mathtt{w}&i\neq 0,\end{array}\right.$$

(11)

where $M_{I}=\sum_{i}g_{i}$ is the total number of states in the ensemble, and $\mathtt{w}\in[0,1/M_{I}]$ is the weight parameter. The excitation energy $\omega_{I}=\mathcal{E}_{I}-\mathcal{E}_{0}$ can be written as a correction to the ground-state KS excitation energy:

$${\omega_{I}=\omega_{I}^{\text{KS}}+\frac{1}{g_{I}}\left.\frac{d}{dw}(E_{\text{Hxc},I,\mathtt{w}}-E_{\text{Hxc},I-1,\mathtt{w}})\right|_{\mathtt{w}=0}.}$$

(12)

The ensemble density reduces to the ground-state density at $\mathtt{w}=0$, so Eq. (12) only requires a self-consistent ground-state KS calculation.

GLPTGL93 ; GL94 is developed for ground-state DFT as a systematic approach for developing orbital-dependent xc functionals. The perturbation is along the adiabatic connection with the KS Hamiltonian as the zeroth order and density fixed. The perturbation is the difference between the interacting and KS Hamiltonians $\hat{V}_{\rm ee}-\hat{V}_{\scriptscriptstyle\rm HXC}$, so one can extract approximations to $E_{\scriptscriptstyle\rm XC}$ and $v_{\scriptscriptstyle\rm XC}$ from the perturbative corrections. GLPT differs from the Rayleigh-Schrödinger perturbation theory since the perturbation is dependent on the coupling strength, and perturbative corrections to the density vanish at all orders. GLPT yields the exact exchange (EXX) approximation which is widely used in both DFT and TDDFT.

We apply GLPT on EDFT in this paper. GLPT has to be applied along a different adiabatic connection with the ensemble density fixed, with the EKS Hamiltonian as the zeroth-order Hamiltonian. The Schrödinger equation with the electron-electron interaction scaled by the coupling constant $\lambda\in[0,1]$ is

$${(\hat{T}+\hat{V}_{\lambda,{\{w\}}}+\lambda\hat{V}_{\rm ee})}\left|\Phi_{ik,\lambda,{\{w\}}}\right>=\mathcal{E}_{i,\lambda,{\{w\}}}\left|\Phi_{ik,\lambda,{\{w\}}}\right>,$$

(13)

where $\hat{V}_{\lambda,{\{w\}}}=\sum_{i}v_{\lambda,{\{w\}}}({\bf r}_{i})$, and the existence of $v_{\lambda,{\{w\}}}$ is assumed. The scaled ensemble energy is

$$\begin{split}E_{\lambda,{\{w\}}}&=\operatorname{tr}\{\hat{D}_{\lambda,{\{w\}}}\hat{H}_{\lambda}\}=\sum_{i=0}^{I}w_{i}g_{i}\mathcal{E}_{i,\lambda,{\{w\}}}\\ &=T_{\text{s},\lambda,{\{w\}}}+\int\mathrm{d}^{3}r\,n_{\lambda,{\{w\}}}({\bf r})\left[v_{\rm ext}({\bf r})+v_{\text{Hxc},\lambda,{\{w\}}}({\bf r})\right],\end{split}$$

(14)

where $\hat{D}_{\lambda,{\{w\}}}=\sum_{i=0}^{I}w_{i}\sum_{k=1}^{g_{i}}\left|\Phi_{ik,\lambda,{\{w\}}}\right>\left<\Phi_{ik,\lambda,{\{w\}}}\right|$. Treating $\lambda\hat{V}_{\rm ee}+\hat{V}_{\lambda,{\{w\}}}-\hat{V}_{\text{s},{\{w\}}}$ as the perturbation and expanding $v_{\lambda,{\{w\}}}({\bf r})$ as $\sum_{p=0}^{\infty}\lambda^{p}v^{(p)}_{\{w\}}({\bf r})$, we obtain the formula for the ensemble PT2 correlation energy:

$$E^{\text{PT2}}_{\text{c},{\{w\}}}=\sum_{i=0}^{I}w_{i}\sum_{k=1}^{g_{i}}\left\langle\Phi_{ik,{\{w\}}}\left|\hat{V}_{\rm ee}+\hat{V}^{(1)}_{\{w\}}\right|\Phi^{(1)}_{ik,{\{w\}}}\right\rangle.$$

(15)

Details of the derivations is available in Appendix A.

Degenerate perturbation theoryD61 ; L14 requires that the zeroth-order wavefunctions of an EKS multiplet diagonalizes $\hat{H}^{\prime}$, and this is satisfied by the EKS wavefunctions in Eq. (7) since they are eigenfunctions of the spatial and spin symmetry operators that commute with $\hat{H}$. $\left|\Phi^{(1)}_{ik,{\{w\}}}\right>$ in Eq. (15) is then

$$\left|\Phi^{(1)}_{ik,{\{w\}}}\right>=\sum_{\tilde{j}\neq\tilde{i}}^{\infty}\sum_{q=1}^{\tilde{g}_{\tilde{j}}}\frac{\left\langle\tilde{\Phi}^{(0)}_{\tilde{j}q,{\{w\}}}\left|\hat{V}_{\rm ee}+\hat{V}^{(1)}_{\{w\}}\right|\Phi^{(0)}_{ik,{\{w\}}}\right\rangle}{\mathcal{E}^{\text{KS}}_{i,{\{w\}}}-\mathcal{E}^{\text{KS}}_{j,{\{w\}}}}\left|\tilde{\Phi}^{(0)}_{\tilde{j}q,{\{w\}}}\right>,$$

(16)

where $\mathcal{E}^{\text{KS}}_{i,{\{w\}}}$ is the EKS energy of EKS multiplet $\tilde{i}$.

The OEP equations for the perturbative corrections of the potential can be derived by requiring the perturbative corrections of the density to vanishGL93 . However, these equations involve weighted sums over states and are more difficult to solve than their counterparts in ground-state DFT. Approximations to $v^{(1)}_{\{w\}}$ can be derived and are provided in the supplemental materialsupplemental . While approximations to $v^{(2)}_{\{w\}}$ can be derived similarly in principle, it is more cumbersome as it would contain functional derivatives of $v^{(1)}_{\{w\}}$, so self-consistent EKS calculation with ensemble PT2 correlation can be unhandy. We also encounter divergences related to the ensemble derivative discontinuityL95 in the OEP calculations. We therefore choose to use the DEC method to assess the performance of the ensemble PT2 correlation in this paper to avoid solving the EKS equation self-consistently.

The orbital-dependent Eqs. (9) and (15) are implicit density functionals. In the DEC method, these equations represent $E_{\text{Hx},\mathtt{w}}[n^{\text{KS}}_{\mathtt{w}}]$ and $E^{\text{PT2}}_{\text{c},\mathtt{w}}[n^{\text{KS}}_{\mathtt{w}}]$ instead of $E_{\text{Hx},\mathtt{w}}[n]$ and $E^{\text{PT2}}_{\text{c},\mathtt{w}}[n]$ required by Eq. (12). Since the $\mathtt{w}$-dependences of the functional and of the EKS density are inseparable in these equations, we calculate the $\mathtt{w}$-derivatives in Eq. (12) byYPBU17

$${\left.\frac{dE_{\text{Hxc},I,\mathtt{w}}}{d\mathtt{w}}\right|_{\mathtt{w}=0}=\mathcal{D}\left(E_{\text{Hxc},I,\mathtt{w}}\right)-\int\mathrm{d}^{3}r\,v_{\scriptscriptstyle\rm HXC}({\bf r})\mathcal{D}\left(n^{\text{KS}}_{I,\mathtt{w}}\right),}$$

(17)

where $v_{\scriptscriptstyle\rm HXC}$ is the ground-state Hxc potential, and the shorthand $\mathcal{D}$ means

	$\displaystyle\mathcal{D}\left(E_{\text{Hxc},I,\mathtt{w}}\right)$	$\displaystyle=\left.\frac{\partial E_{\text{Hxc},I,\mathtt{w}}[n^{\text{KS}}_{I,\mathtt{w}}[\{\phi_{\mu}\}]]}{\partial\mathtt{w}}\right|_{\begin{subarray}{c}\mathtt{w}=0\\ \{\phi_{\mu}\}=\{\phi^{\text{KS}}_{\mu}\}\end{subarray}},$		(18)
	$\displaystyle\mathcal{D}\left(n^{\text{KS}}_{I,\mathtt{w}}\right)$	$\displaystyle=\left.\frac{\partial n^{\text{KS}}_{I,\mathtt{w}}[\{\phi_{\mu}\}]({\bf r})}{\partial\mathtt{w}}\right|_{\begin{subarray}{c}\mathtt{w}=0\\ \{\phi_{\mu}\}=\{\phi^{\text{KS}}_{\mu}\}\end{subarray}},$		(19)

with $\{\phi_{\mu}\}$ and $\{\phi^{\text{KS}}_{\mu}\}$ being a complete set of one-electron orbitals and the ground-state KS orbitals respectively, and $n^{\text{KS}}_{I,\mathtt{w}}[\{\phi_{\mu}\}]$ being the EKS density of the ensemble defined by Eq. (11) with $I$ being the highest energy multiplet constructed with orbitals $\{\phi_{\mu}\}$.

The orbitals are held fixed for the derivatives in Eq. (17) to avoid the otherwise required self-consistent EKS calculations, since the $\mathtt{w}$-dependence of $\left|\tilde{\Phi}_{\tilde{i}p,\mathtt{w}}\right>$ and $E^{\text{KS}}_{i,\mathtt{w}}$ originates from that of the EKS orbitals. $\mathcal{D}(E_{\text{Hxc},I,\mathtt{w}})$ is then the direct differentiation of Eqs. (9) and (15). Its PT2 part is

$$\mathcal{D}\left(E^{\text{PT2}}_{\text{c},I,\mathtt{w}}\right)=\sum_{i=0}^{I}\mathcal{D}(w_{i})\sum_{k=1}^{g_{i}}\sum_{\tilde{j}\neq\tilde{i}}^{\infty}\sum_{q=1}^{\tilde{g}_{\tilde{j}}}\frac{\left|\mathcal{M}_{\tilde{j}q}^{ik}\left(\hat{V}_{\rm ee}-\hat{V}_{\scriptscriptstyle\rm HX}\right)\right|^{2}}{\mathcal{E}^{\text{KS}}_{i,\mathtt{w}}-\mathcal{E}^{\text{KS}}_{j,\mathtt{w}}}\\ +\frac{1}{g_{0}}\sum_{k=1}^{g_{0}}\sum_{\tilde{j}\neq\tilde{0}}^{\infty}\sum_{q=1}^{\tilde{g}_{\tilde{j}}}\frac{\left\{\mathcal{M}_{\tilde{j}q}^{ik}\left(\hat{V}_{\rm ee}-\hat{V}_{\scriptscriptstyle\rm HX}\right)^{*}\mathcal{M}_{\tilde{j}q}^{0k}\left[-\mathcal{D}\left(\hat{V}_{\text{Hx},\mathtt{w}}\right)\right]+\text{c.c.}\right\}}{\mathcal{E}^{\text{KS}}_{0}-\mathcal{E}^{\text{KS}}_{j}}$$

(20)

where $\mathcal{D}(w_{i})=dw_{i}(\mathtt{w})/d\mathtt{w}|_{\mathtt{w}=0}$, $\mathcal{D}(\hat{V}_{\text{Hx},\mathtt{w}})=\partial\hat{V}_{\text{Hx},\mathtt{w}}[n_{\text{s},\mathtt{w}}[\{\phi_{\mu}\}]]/\partial\mathtt{w}|_{\{\phi_{\mu}\}=\{\phi^{\text{KS}}_{\mu}\},\mathtt{w}=0}$, and the notation $\mathcal{M}$ means

	$\displaystyle\mathcal{M}_{\tilde{j}q}^{ik}\left(\hat{V}_{\rm ee}-\hat{V}_{\scriptscriptstyle\rm HX}\right)$	$\displaystyle=\left\langle\tilde{\Phi}_{\tilde{j}q}\left|\hat{V}_{\rm ee}-\hat{V}_{\scriptscriptstyle\rm HX}\right|\Phi_{ik}\right\rangle,$		(21)
	$\displaystyle\mathcal{M}_{\tilde{j}q}^{0k}\left[-\mathcal{D}\left(\hat{V}_{\text{Hx},\mathtt{w}}\right)\right]$	$\displaystyle=\left\langle\tilde{\Phi}_{\tilde{j}q}\left|-\mathcal{D}\left(\hat{V}_{\text{Hx},\mathtt{w}}\right)\right|\Phi_{ik}\right\rangle.$		(22)

We ignore the second term of Eq. (20) in the calculations in this paper since it is usually small (Appendix B). This approximation allows us to greatly simplify the calculation by using a bi-ensemble comprised of the ground state and the $I$-th multiplet, which is equivalent to the full ensemble due to the functional form of the first term of Eq. (20).YPBU17 Eq. (20) apparently indicates that the number of orbitals required for convergence determines the scaling of the computational cost.

We solve the interacting Schrödinger equation for 1D model systems by direct diagonalization on a grid. All 1D model systems have two electrons so that the exact $v_{\scriptscriptstyle\rm X}(x)$ can be obtained as $-v_{\scriptscriptstyle\rm H}(x)/2$. The second term of Eq. (20) vanishes for two-electron systems, so the 1D PT2 calculations are not approximated. Numerically exact KS ground state is obtained by inverting the ground-state KS equationYTPB14 .

We use an 1D Hooke’s atom as the model system for double excitation MZCB04 ; YPBU17 . The system has $v_{\rm ext}(x)=x^{2}/2$ and contact interaction $v_{\rm ee}(x,x^{\prime})=0.2\delta(x-x^{\prime})$, where $\delta$ is the Dirac $\delta$ function. The doubly excited states in this system are close to singly excited states due to the weak electron-electron interaction. We use a large $20001\times 20001$ grid with grid-point spacing $0.001$ and $x\in[-10,10]$ to ensure that the effect of the grid boundary on all orbitals is negligible. Table 1 lists the errors in the excitation energies.

$\Delta\omega_{I}^{\text{EEXX}}$ shows the error in excitation energies due to the lack of correlation. With PT2 correlation, EEXX+PT2 in Table 1 improves the accuracy of excitation energies in most cases, including both doubly-excited states. Due to the perturbative nature of the ensemble PT2 correlation, however, improvement in accuracy is not guaranteed, as seen in the first excited state.

The correlation contribution of the two terms of Eq. (17) can be obtained from Table 1 by subtracting column 1(EEXX) from columns 2(EEXX+$E_{\scriptscriptstyle\rm C}^{\text{PT2}}$) and 3(EEXX+$v_{\scriptscriptstyle\rm C}$), which shows that the correlation effect in the 1D Hooke’s atom is dominated by the first term of Eq. (17). This dominance allows us to estimate the PT2 errors in excitation energies with this first term alone while ignoring the effect of using exact $v_{\scriptscriptstyle\rm C}$: comparing $\omega_{I}^{\text{exact}}-\omega_{I}^{\text{EEXX}+v_{\scriptscriptstyle\rm C}}$ (the effect of this term of the exact functional) and $\omega_{I}^{\text{EEXX}+E_{\scriptscriptstyle\rm C}^{\text{PT2}}}-\omega_{I}^{\text{EEXX}}$ (that of the PT2 approximation) yields $163\%$, $26\%$, $19\%$, $19\%$, $62\%$ for the 5 states in Table 1. These large errors for PT2 indicate that higher-order terms are probably needed to properly describe the ensemble correlation energy.

We use a system with two potential wells in a box as the model system for charge-transfer excitation. The external potential is

$$v_{\rm ext}(x)=\left\{\begin{array}[]{ll}\infty&x\in(-\infty,0]\cup[6.5,\infty),\\ 0&x\in(0,1)\cup(5,6.5),\\ 20&x\in[1,5].\end{array}\right.$$

(23)

and the electron-electron interaction is the soft-Coulomb potential

$$v_{\rm ee}(x,x^{\prime})=\frac{1}{\sqrt{(x-x^{\prime})^{2}+\alpha}}.$$

(24)

We set $\alpha=1$ in order to carry out adiabatic TDDFT calculations with the 1D local-density approximation (LDA)HFCV11 kernel for comparison. We orbtain the numerically exact KS system on a $1301\times 1301$ grid with grid-point spacing $0.005$ and $x\in[0,6.5]$. Excitations to the first (triplet) and second (singlet) excited states are charge-transfer excitations, both of which correspond to one electron transferring from the right potential well to the left one. The distance between the two wells lead to negligible overlap between the involved KS orbitals, so the TDDFT couplings between the initial and final states vanish for local or semi-local approximated xc kernels. The corresponding TDDFT excitation energies would be very close to the KS ones.

We list the error in the KS excitation energy in Table 2. Since the singlet-triplet splitting between the first two states is very small, Table 2 only lists the calculation result for the first excited state. The correlation effect dominates the correction to $\omega_{I}^{\text{KS}}$ for the charge-transfer box as well. Both DEC/EEXX and TDDFT/ALDA fail to correct the KS excitation energy due to non-overlapping orbitals leading to vanishing matrix elements. Unlike the 1D Hooke’s atom, the correlation contribution from the $\mathtt{w}$-dependence of the orbital-dependent $E^{\text{PT2}}_{\text{c},\mathtt{w}}$ [first term of Eq. (17)] only changes the result slightly due to either vanishing matrix elements or large energy differences in Eq. (20). The contribution from $v_{\scriptscriptstyle\rm C}$ dominates the correlation effect. This suggests that the two correlation contributions in Eq. (17) can be regarded as representing the local and the non-local correlation effects, respectively.

We also calculate the excitation energies of an 1D flat box, where both terms of Eq. (17) contribute significantly to the correlation effect (Appendix C). We find that the inclusion of ensemble PT2 correlation improves the exchange-only EEXX excitation energies for most of the states of the tested 1D systems. However, the PT2 correction to EEXX is not always in the correct direction. The magnitude of improvement also varies a lot for different states.

We carry out DEC calculations on He and Be atoms using numerically exact KS potentialsUG93 ; UG94 , and compare the excitation energies to the experimental valuesNIST_ASD . Since we solve the ground-state KS orbitals on a radial grid, the higher-energy orbitals are not well represented due to tails being truncated at the grid boundary. This leads to the ordering of higher-energy orbitals being different from that of the exact KS system (see He PT2 convergece curve in supplemental materialsupplemental ). Nevertheless, this should not affect the PT2 results since the orbitals form a complete basis set with respect to the grid, and convergence is achieved fairly quickly for PT2 (Sec. III.3).

Table 3 lists the errors in the excitation energies for He and Be. The performance of the ensemble PT2 correlation energy functional is disappointing in atoms. Unlike 1D systems, EEXX+PT2 excitation energies are in general worse than those of EEXX despite a few exceptions. This includes the He atom as well, which is a two-electron system as the 1D model systems. The ground-state PT2 correlation is also found to be problematic when used directlyMWY05 , and the failure is attributed to the algebraic structure of KS orbitals. However, this is unlikely the reason for the failure seen in Table 3 since the DEC method does not involve self-consistent ensemble PT2 calculations.

The reason for larger error of EEXX+PT2 might be related to the two terms of Eq. (17) being incompatible. For EEXX, these two terms are close in magnitude and have different signsYPBU17 . Due to their cancellation, the resulting corrections to KS excitation energies have much smaller magnitudes. If these two terms are incompatible, the error can be large since Eq. (17) subtracts two large numbers, which can be seen from the EEXX+$v_{\scriptscriptstyle\rm C}$ results. This might be the case for EEXX+PT2 shown in Table 3, where we use the exact $v_{\scriptscriptstyle\rm C}$ in Eq. (17) instead of the OEP $v_{\scriptscriptstyle\rm C}^{\text{PT2}}$. As an example to demonstrate the effect of incompatible $E_{\text{c},\mathtt{w}}$ and $v_{\scriptscriptstyle\rm C}$, using the LDA correlation potentialPW92 in an EEXX calculation of the first excitation energy of He yields an error of -55.98 mH, which is much greater than the one in Table 3. Since the correlation is dominated by one of the two terms for the 1D Hooke’s atom and charge-transfer box, the error due to this incompatibility is small (but the error due to PT2 itself can still be large), so the overall accuracy of EEXX+PT2 only depends on the dominant term. This problem is a disadvantage of orbital-dependent Hxc energy functionals and does not affect explicit density functionals. This failure also indicates that the OEP $v_{\scriptscriptstyle\rm C}^{\text{PT2}}$ is not close to the exact $v_{\scriptscriptstyle\rm C}$ for atoms, which agrees with the large relative errors due to PT2 in the 1D Hooke’s atom, showing that higher-order terms are probably needed for correlation.

EEXX+PT2(no single) results for the 1D Hooke’s atom and 1D charge-transfer box are close to the DEC/EEXX+PT2 results, but the error is larger for atoms. Despite this, the error introduced by this approximation is still small comparing with the excitation energies (see supplemental materialsupplemental and Ref. NIST_ASD ), suggesting that the ground-state approximation of neglecting singly-excited matrix elements is also applicable in EDFT.

Although the $\tilde{j}$- and $q$-sums in Eq. (20) sums over all KS states, these sums can be truncated by only including a finite number of KS orbitals, since the contributions from higher-energy orbitals vanish due to the KS energy differences in the denominators. We check the convergence of the EEXX+PT2 excitation energies with respect to the number of KS orbitals, which are plotted in Fig. 1 for the 1D Hooke’s atom.

Fig. 1 shows that convergence is achieved for all EEXX+PT2 excitation energies of the 1D Hooke’s atom with only a few KS orbitals. Higher-energy orbitals have a small impact on excitation energies as the KS energy differences become large. Due to the near-degeneracies between KS states $(2,2)$ and $(1,3)$ and between $(2,3)$ and $(1,4)$, the changes of the $I=2$ and $I=4$ excitation energies are big when orbitals 3 and 4 are included, respectively. The excitation energies of the 1D charge-transfer box and 1D flat box (see supplemental materialsupplemental ) also converge quickly with respect to the number of orbitals, since the orbital energies of these systems increase rapidly.

Fig. 2 plots the EEXX+PT2 excitation energies of the Be atom with different numbers of KS orbitals included. The He atom have similar trends (see supplemental materialsupplemental ). Comparing with 1D cases, higher-energy orbitals can affect the excitation energies significantly since the orbital energies form a Rydberg series. More orbitals are needed to achieve convergence in 3D atoms.

The excitation energies in Fig. 2 only change notably when orbitals with certain angular momentum quantum numbers are included, since otherwise many of the newly included matrix elements in Eq. (20) would vanish due to symmetry. Therefore one can achieve convergence with fewer orbitals for atoms by ignoring the orbitals whose angular momentum quantum numbers differ too much from those of the orbitals in the KS ensemble. As an example, the EEXX+PT2 excitation energy of the ¹S(1s²2s3s) state of Be calculated with all orbitals (upto 6g) only differ from that calculated with only s and p orbitals (upto 6p) by 0.007 mH.