Kohn-Sham Density Functional Theory with Complex, Spin-Restricted Orbitals: Accessing a New Class of Densities without the Symmetry Dilemma

TL;DR

This paper introduces complex, spin-restricted orbitals in Kohn-Sham DFT, enabling access to new densities and improving triplet-singlet gap predictions, while highlighting limitations of current functionals.

Contribution

It demonstrates that complex, spin-restricted orbitals can access new density classes and improve certain predictions, offering a novel approach to overcome the symmetry dilemma in DFT.

Findings

cRKS orbitals improve triplet-singlet gap predictions

RKS and UKS yield qualitatively wrong densities

Modern functionals show limited improvement with cRKS densities

Abstract

We show that using complex, spin-restricted orbitals (cR) in Kohn-Sham density functional theory (KS-DFT) allows one to access a new class of densities that is not accessible by either spin-restricted (RKS) or spin-unrestricted (UKS) orbitals. We further show that the real part of a cRKS density matrix can be non-idempotent when the imaginary part of the density matrix is not zero. Using cRKS orbitals shows significant improvements in the triplet-singlet gaps of the TS12 benchmark set for the SPW92, PBE, BLYP, and TPSS functionals. Moreover, it was shown that RKS and UKS yield qualitatively wrong charge densities and spin densities, respectively, leading to worse energetics. We demonstrate that the modern functionals SCAN, MN15-L, and B97M-V show surprisingly no improvement even with a qualitatively more accurate density from cRKS orbitals. To this end, our work not only paves the way to escape the symmetry dilemma whenever there exists a cRKS solution, but also suggests a new route to design better approximate XC functionals.

Figures

Tables

Table 1. Table 1: The deviation (kcal/mol) with respect to experimental values in Δ ​ E T-S ( = E S − E T ) annotated Δ subscript 𝐸 T-S absent subscript 𝐸 𝑆 subscript 𝐸 𝑇 \Delta E_{\text{T-S}}(=E_{S}-E_{T}) of the SPW92, PBE, BLYP, and TPSS functionals obtained with restricted (RKS), unrestricted (UKS), and complex, restricted (cRKS) orbitals. RMSD stands for root-mean-square-deviation and MSD stands for mean-signed-deviation.

	SPW92	PBE	BLYP	TPSS
	RKS
MSD	10.85	13.09	9.94	13.85
RMSD	11.46	13.53	10.34	14.42
	UKS
MSD	-13.70	-16.57	-17.61	-17.26
RMSD	14.46	17.61	18.64	18.19
	cRKS
MSD	-1.23	3.00	0.90	7.94
RMSD	2.19	3.41	1.91	8.63

Table 2. Table 2: The energy differences (kcal/mol) between cRKS and RKS solutions of \ce O2 (i.e., Δ ​ E = Δ ​ E cRKS − Δ ​ E RKS Δ 𝐸 Δ subscript 𝐸 cRKS Δ subscript 𝐸 RKS \Delta E=\Delta E_{\text{cRKS}}-\Delta E_{\text{RKS}} ) for each component in Eq. ( 13 ). Δ ​ E Coul = Δ ​ E J + Δ ​ E V Δ subscript 𝐸 Coul Δ subscript 𝐸 𝐽 Δ subscript 𝐸 𝑉 \Delta E_{\text{Coul}}=\Delta E_{J}+\Delta E_{V} , Δ ​ E X Δ subscript 𝐸 𝑋 \Delta E_{X} , and Δ ​ E C Δ subscript 𝐸 𝐶 \Delta E_{C} are the energy differences in classical Coulomb energy, exchange energy, and correlation energy, respectively.

	$\Delta E_T$	$\Delta E_V$	$\Delta E_J$	$\Delta E_{\text{Coul}}$	$\Delta E_X$	$\Delta E_C$	$\Delta E_{\text{XC}}$
SPW92	12.89	-45.70	11.65	-34.05	7.88	0.20	8.09
PBE	10.69	-39.57	7.06	-32.51	12.13	-1.38	10.76
BLYP	9.96	-36.43	4.09	-32.33	13.05	-0.38	12.67
TPSS	5.95	-26.95	-2.28	-29.22	18.45	-1.65	16.80

Table 3. Table 3: Same as Table 1 except we have the SCAN, MN15-L, and B97M-V functionals presented here. Note that MN15-L and B97M-V have no stable cRKS solutions for a half of the data points (i.e., \ce C, \ce NF, \ce NH, \ce NO-, \ce O2, and \ce O).

	SCAN	MN15-L	B97M-V
	RKS
MSD	19.50	10.64	11.82
RMSD	19.82	11.16	12.22
	UKS
MSD	-16.24	-10.79	-13.69
RMSD	17.74	11.97	14.71
	cRKS
MSD	15.55	10.50	10.34
RMSD	16.39	11.00	11.15

Equations

$$\Pi (\mathbf{r})={\rho }_{\alpha }(\mathbf{r}){\rho }_{\beta }(\mathbf{r}).$$

$$m(\mathbf{r})=\rho (\mathbf{r}){\left(1-\frac{4\Pi (\mathbf{r})}{\rho (\mathbf{r}{)}^2}\right)}^{1/2}$$

$$\mathcal{L}[\mathbf{P}]=$$

$$+\mathrm{tr}(\mathbf{B}({\mathbf{P}}^{†}-\mathbf{P}))$$

$$\mathbf{P}={\mathbf{C}}_{\text{occ}}{\mathbf{C}}_{\text{occ}}^{†}$$

$$\mathbf{P}=\mathbf{X}+i\mathbf{Y}$$

$$\mathrm{tr}(\mathbf{X})$$

$$\mathrm{tr}(\mathbf{Y})$$

$$\mathbf{X}$$

$$\mathbf{Y}$$

$$\mathbf{X}$$

$$\mathbf{Y}$$

$$\rho (\mathbf{r})=2\sum _{\mu \nu }{\eta }_{\mu }(\mathbf{r}){\eta }_{\nu }(\mathbf{r})P_{\mu \nu }$$

$$E_{\text{KS-DFT}}[\mathbf{X}]=$$

$$+E_{XC}\left[\rho (\mathbf{r}),\nabla \rho (\mathbf{r}),\cdot \cdot \cdot \right]+E_{nn}$$

$$E_T\left[\mathbf{X}\right]=2\mathrm{tr}\left(\mathbf{XT}\right)$$

$$E_V\left[\rho (\mathbf{r})\right]=-\sum _I{Z}_I{\int }_{\mathbf{r}}\frac{\rho (\mathbf{r})}{||\mathbf{r}-{\mathbf{R}}_I|{|}_2}$$

$$E_J\left[\rho (\mathbf{r})\right]=\frac{1}{2}{\int }_{{\mathbf{r}}_1}{\int }_{{\mathbf{r}}_2}\frac{\rho ({\mathbf{r}}_1)\rho ({\mathbf{r}}_2)}{||{\mathbf{r}}_1-{\mathbf{r}}_2|{|}_2},$$

$$|\Psi ⟩=\mathcal{A}[\xi \xi \left(\alpha \beta \right)]=\mathcal{A}[\left(\Pi +\Omega \right)\left(\alpha \beta \right)]$$

$$\mathbf{X}=\left[\begin{array}{cc}{\cos }^2(\theta )& 0\\ 0& {\sin }^2(\theta )\end{array}\right]$$

Full text

Introduction In Kohn-Sham density functional theory (KS-DFT), it is clear that for systems with antiferomagnetically coupled open-shell electrons, spin-unrestricted KS-DFT (UKS) behaves better than spin-restricted KS-DFT (RKS) at the expense of spin-polarization.Gunnarsson and Lundqvist (1976); Seminario (1994); Davidson (1998) For finite-sized systems, this spin-polarization is a manifestation of the lack of strong (or static) correlation in a single-determinant wavefunction.Lykos and Pratt (1963); Perdew et al. (1995) This is known as a symmetry dilemma in KS-DFT, so called because exact solutions for finite systems should not break spin-symmetry yet the unphysical, spin-symmetry-broken approximate solutions often provide lower energies than symmetric approximate solutions.

As an attempt to understand this symmetry dilemma, Perdew and co-workers proposed the use of “on-top” pair density,Perdew et al. (1995, 1997) $\Pi(\mathbf{r})$, which is defined as

[TABLE]

Instead of working with spin density variables, $\rho_{\alpha}(\mathbf{r})$ and $\rho_{\beta}(\mathbf{r})$, the fundamental variables are now the charge density, $\rho(\mathbf{r})=\rho_{\alpha}(\mathbf{r})+\rho_{\beta}(\mathbf{r})$, and the magnetization density, $m(\mathbf{r})$, which is a function of $\rho(\mathbf{r})$ and $\Pi(\mathbf{r})$. One can accomplish the one-to-one mapping between $m(\mathbf{r})$ and $\Pi(\mathbf{r})$ following

[TABLE]

Evidently, inserting Eq. (1) into Eq. (2) yields a more familiar definition of the magnetization density, $m(\mathbf{r})=\rho_{\alpha}(\mathbf{r})-\rho_{\beta}(\mathbf{r})$. With $\rho(\mathbf{r})$ and $m(\mathbf{r})$, one can back out $\rho_{\alpha}(\mathbf{r})$ and $\rho_{\beta}(\mathbf{r})$ and these can be used to compute the exchange-correlation (XC) potential, $v_{\text{xc}}(\mathbf{r})$. This framework allows one to treat $m(\mathbf{r})$ as an auxiliary variable that is closely related to the on-top pair density. Therefore, it helps to explain spin-polarization in KS-DFT in terms of the on-top pair density of a wavefuntion. While this is certainly useful in understanding this symmetry dilemma, it is somewhat unsatisfying that one has to interpret the magnetization density merely as an auxiliary variable. Moreover, it does not remedy problems arising from evaluating $v_{\text{xc}}(\mathbf{r})$ with spin densities.

In this Letter, we describe a way to access a new class of charge and spin densities that is not accessible by spin-polarization. This is achieved by breaking time-reversal symmetry and complex symmetry of the KS-DFT determinant.Fukutome (1981); Stuber and Paldus (2003); Small et al. (2015); Lee and Head-Gordon (2019) We will refer to this symmetry breaking as “complex-polarization”. For semi-local functionals, we show that this complex-polarization does not pose any dilemma in obtaining charge and spin densities. We illustrate how this is achieved and demonstrate the numerical performance of popular semi-local functionals on chemical systems where complex-polarization is essential to obtain correct charge and spin densities.

Theory We start from the KS-DFT Lagrangian,Taube (2010)

[TABLE]

where $N$ is the number of electrons, $\mathbf{P}$ is a one-particle density matrix (1PDM), $\mu$, $\mathbf{A}$, and $\mathbf{B}$ are Lagrange multipliers for constraining the trace, idempotency, and hermiticity of $\mathbf{P}$, respectively. The idempotency of $\mathbf{P}$ guarantees the non-interacting nature of the KS-DFT problem. For simplicity, we also assume that the computational basis is orthogonal.

These three constraints are naturally imposed by the definition of $\mathbf{P}$,

[TABLE]

where $\mathbf{C}_{\text{occ}}$ is the occupied molecular orbital (MO) coefficient matrix. However, when minimizing Eq. (3) with respect to $\mathbf{P}$ directly as is done in linear-scaling DFT approaches, it is necessary to consider these constraints explicitly or impose them through density matrix purification techniques.Kussmann et al. (2013) We write $\mathbf{P}$ in terms of its real and imaginary components,

[TABLE]

We can then rewrite these constraints as follows: the trace condition leads to

[TABLE]

the idempotency becomes

[TABLE]

and the hermiticity yields

[TABLE]

Based on Eq. (11), we conclude that $\mathbf{Y}$ is antisymmetric and therefore it automatically satisfies Eq. (7). This is an important observation since $\rho(\mathbf{r})$ is computed from $\mathbf{P}$ following

[TABLE]

where $\eta_{\mu}(\mathbf{r})$ is a real-valued computational basis function. This real-valuedness is sufficient for studying bound states of molecules and indeed most basis sets for finite-sized systems are real-valued. Since $\eta_{\mu}(\mathbf{r})\eta_{\nu}(\mathbf{r})$ is a symmetric tensor, an antisymmetric tensor, $\mathbf{Y}$, does not contribute to $\rho(\mathbf{r})$. In other words, it is equivalent to write $\rho(\mathbf{r})=2\sum_{\mu\nu}\eta_{\mu}(\mathbf{r})\eta_{\nu}(\mathbf{r})X_{\mu\nu}$.

The key insight here is that by having non-zero $\mathbf{Y}$, the idempotency constraint on $\mathbf{X}$ can be relaxed to Eq. (8) and Eq. (9). In other words, $\mathbf{X}$ does not have to satisfy $\mathbf{X}=\mathbf{X}^{2}$ as long as $\mathbf{Y}$ is non-negligible. Such a density matrix $\mathbf{P}$ with $\operatorname{tr}(\mathbf{Y}^{T}\mathbf{Y})\neq 0$ is referred to as a “fundamentally complex” density matrix. It is fundamentally complex in the sense that no unitary rotation of a real-valued MO coefficient matrix in the complex plane can represent a fundamentally complex density matrix.Small et al. (2015) Consequently, the use of fundamentally complex density matrices results in a broader class of densities which may come from a non-idempotent density matrix $\mathbf{X}$. It is not obvious whether lifting this idempotency constraint on $\mathbf{X}$ would always yield a lower energy solution. In fact, the energy lowering turns out to be quite rare to observe and systems with a high point group symmetry tend to exhibit this energy lowering.

For semi-local functionals, $E_{\text{KS-DFT}}$ depends only on $\mathbf{X}$ and reads

[TABLE]

where the kinetic energy is defined as

[TABLE]

with $\mathbf{T}$ is the kinetic energy matrix in the computational basis ($T_{\mu\nu}=\langle\mu|-\frac{1}{2}\nabla^{2}|\nu\rangle$), the nuclear-electron attraction energy is

[TABLE]

with $I$ denoting the nuclei (or ions), the electron-electron repulsion energy reads

[TABLE]

$E_{XC}$ denotes the XC energy, and $E_{nn}$ is the nuclear-nuclear repulsion energy. We emphasize that as far as Eq. (13) is concerned there is no symmetry breaking associated with $\mathbf{X}$ and $\rho$. Furthermore, there is no auxiliary variable introduced in the energy evaluation.

Beyond permitting access to densities not describable by RKS and even UKS, one other point should be mentioned. Just as $\mathbf{X}$ is non-idempotent, the cRKS determinant is inherently multiconfigurational (MC) as shown in ref. 9. A single complex orbital $\xi$ is parametrized by $\theta$ to have an arbitrary linear combination of real $\eta$ and imaginary $\bar{\eta}$ orbitals, $\xi=\cos(\theta)\eta-i\sin(\theta)\bar{\eta}.$ A two-electron closed-shell determinant made of a complex spatial orbital $\xi_{i}$ follows

[TABLE]

where $\mathcal{A}$ is the antisymmetrizer including a normalization factor, $\Pi=\cos^{2}(\theta)\eta\eta-\sin^{2}(\theta)\bar{\eta}\bar{\eta}$, and $\Omega=-i\sin(\theta)\cos(\theta)(\eta\bar{\eta}+\bar{\eta}\eta)$. These two spatial wavefunctions, $\Pi$ and $\Omega$, highlight the MR character of a cRKS single-determinant $|\Psi\rangle$. $\Pi$ is a two-configuration wavefunction (lower energy) and $\Omega$ is a open-shell wavefunction (higher energy). The competition between energy-lowering via $\Pi$ and an energy penalty through $\Omega$ determines whether an RKS electron pair complexifies to cRKS. In the $\eta$ basis, the corresponding 2x2 matrix, $\mathbf{X}$ is

[TABLE]

$\theta$ = 0 represents an RKS density matrix whereas $\theta=\pi/4$ yields two half-occupied orbitals. Depending on $\theta$, it is possible to obtain a non-idempotent $\mathbf{X}$. While the hidden MC form in the cRKS determinant is limited, its real 1PDM (for $E_{T}$) and its $\rho(\mathbf{r})$ (for $E_{J}$ and $E_{V}$) can all be consistently evaluated using existing XC functionals. Thus cRKS can be viewed as a limited case of MC-DFT without the formal challengesMiehlich et al. (1997); Gräfenstein et al. (1998); Kurzweil et al. (2009) and practical double-counting problemsMiehlich et al. (1997); Gusarov et al. (2004); Khait and Hoffmann (2004); Li Manni et al. (2014) normally associated with that field.

Results We evaluate the numerical performance of RKS, UKS, and cRKS on the recently developed TS12 benchmark set.Lee and Head-Gordon (2019) This data set contains the experimental singlet-triplet gaps of 12 atoms and molecules: C, NF, NH, NO-, O2, O, PF, PH, S2, S, Si, and SO. The ground states of these molecules are triplets. The lowest singlet states for each system are singlet biradicals and exhibit a spin-restricted Hartree-Fock (RHF) to complex RHF (cRHF) instability. This instability is driven by the underlying point group symmetry which gives rise to the degeneracy between the highest occupied MO (HOMO) and the lowest unoccupied MO (LUMO).

The correct lowest singlet states of these systems should singly-occupy both HOMO and LUMO and thereby obtain charge and spin densities that obey the underlying point group symmetry. Even with the exact XC functional, RKS orbitals are qualitatively wrong in these cases as they doubly-occupy the HOMO, breaking the point group symmetry of the system and leading to a broken-symmetry charge density. Since the exact XC functional yields not only exact energy but also exact charge density, accessing a different class of density matrix other than those from RKS orbitals is necessary. A similar phenomenon was first pointed out by Pople, Gill, and Handy in the context of spin-polarization in open-shell systems.Pople et al. (1995) Indeed, UKS can achieve these single occupations by breaking spin-symmetry. However, we will show that spin-polarization is not the way to access the right charge and spin densities in the systems in TS12.

We investigate a local density approximation (LDA) functional, SPW92,Hohenberg and Kohn (1964); Kohn and Sham (1965); Slater (1963); Perdew and Wang (1992) two generalized gradient approximation (GGA) functionals, BLYPBecke (1988); Lee et al. (1988) and PBE,Perdew et al. (1996) and four meta-GGA (mGGA) functionals, TPSS,Tao et al. (2003) SCAN,Sun et al. (2015) MN15-L,Haoyu et al. (2016) and B97M-V.Mardirossian and Head-Gordon (2015) We used aug-cc-pVQZ basis set,Dunning (1989); Kendall et al. (1992) 99 points for XC radial quadrature and 590 points for XC angular quadrature. We only report KS-DFT solutions found to be locally stable based on stability analysis.Bauernschmitt and Ahlrichs (1996); Sharada et al. (2015) The triplet ground states are computed with UKS ($M_{S}=1$) and we focus on the symmetry dilemma of the singlet states. All calculations were performed with a development version of Q-Chem.Shao et al. (2015)

In Table 1, we present the root-mean-square-deviation (RMSD) and mean-signed deviation (MSD) of the relatively old XC functionals SPW92, BLYP, PBE, and TPSS. Using RKS orbitals clearly overestimates the triplet-singlet gap as restricted orbitals cannot describe the open-shell nature of the lowest singlet state of these systems. These RKS energies for singlet states are too high so that the gap is overestimated, as evidenced by the positive MSDs. Using UKS orbitals has the opposite problem. Namely, the gaps are all underestimated. This is due to the undesired mixing of the singlet and triplet states (i.e., spin-contamination). Since the triplet state is lower in energy, the resulting UKS $M_{S}=0$ state energy is too low. The negative MSDs are due to this spin-contamination. The use of cRKS orbitals improves both the MSD and RMSD by more than a factor of 5 in the case of SPW92, PBE and BLYP. The improvement in TPSS is less impressive although it is still a factor of 2 improvement.

This significant improvement is simply due to obtaining $\rho$ from a non-idempotent $\mathbf{X}$ without spin-polarization. Indeed, in all systems in TS12, the resulting $\mathbf{X}$ has exactly two eigenvalues of 0.5 when obtained from these four XC functionals. These two eigenvalues represent the single-occupancy of HOMO and LUMO which is qualitatively correct for these systems.

To understand what drives the energy lowering from RKS to cRKS, we decompose the energy difference between RKS and cRKS into individual energy contributions. We present this energy decomposition for \ceO2 in Table 2. The qualitative behavior discussed here holds for other systems as well. All four XC functionals show an increase in the kinetic energy as well as the XC energy. The increase in the XC energy is driven by the increase in the exchange energy. The correlation energy shows only a small change. The RKS to cRKS instability is driven by the classical Coulomb energy. In particular, a significant energy lowering in the electron-nuclear attraction energy is the driving force of this instability in all systems in the TS12 set.

The results for the more modern mGGA functionals SCAN, MN15-L, and B97M-V, are presented in Table 3. The RKS and UKS results are qualitatively similar to what was observed in the above four functionals. Namely, the overestimation in gaps was observed with RKS while the underestimation in gaps was observed with UKS. SCAN with RKS orbitals shows a particularly poor behavior compared to all the other functionals examined in this work. SCAN becomes comparable to other functionals when used with UKS orbitals. The cRKS results with these three modern functionals are rather surprising due to a qualitatively different behavior from that of the four functionals discussed before. The use of cRKS orbitals for these functionals does not improve the quantitative energetics. Furthermore, MN15-L and B97M-V exhibit no stable cRKS solutions for \ceC, \ceNF, \ceNH, \ceNO-, \ceO2, and \ceO. Additionally, MN15-L shows occupation numbers of (0.61, 0.39) and (0.78, 0.22) for \cePF and \cePH, respectively. Non-idempotent $\mathbf{X}$ is essential to describe the open-shell nature of singlet ground states in these systems. Therefore, with this improved density matrix and $\rho$ therefrom, approximate XC functionals should perform better. The poor performance of these modern mGGA functionals with cRKS orbitals suggests that data such as the TS12 set may be useful for XC functional development.

Lastly, we present density plots to compare the qualitative differences between densities from RKS and cRKS orbitals. As mentioned throughout this Letter, RKS orbitals break the point group spatial symmetry of systems in the TS12 set whereas cRKS orbitals preserve the symmetry.foo The symmetry breaking at the orbital level results in the symmetry breaking in $\rho$. We shall illustrate this point by looking at $\rho$ represented on a real-space grid.

In Figure 1, we present the real-space density of \ceO2 computed with BLYP. The discussion also applies to other systems in the TS12 set and real-space densities are qualitatively similar across different XC functionals. \ceO2 is cylindrically symmetric around its interatomic axis. Therefore, qualitatively correct density should exhibit this symmetry. This is found to be true with the cRKS density in Figure 1 (a). However, in Figure 1 (b), the RKS density is elongated along y-axis and overall not cylindrically symmetric. This is due to the spatial symmetry breaking caused by doubly occupying HOMO which mixes two low-lying singlet states, ${}^{1}\Delta_{g}$ and ${}^{1}\Sigma_{g}^{+}$.Small et al. (2015) The UKS $\alpha$-spin density also breaks the cylindrical symmetry to a smaller extent than does RKS and is elongated along $x$-axis. The UKS charge density (i.e., the sum of $\alpha$ and $\beta$ density) is only $x$-$y$ symmetric and breaks the cylindrical symmetry. The $x$-$y$ symmetry is because the UKS $\beta$-spin density is rotated 90∘ from the $\alpha$ density (i.e., elongated along $y$ axis).

The qualitatively improved UKS charge density (compared to the RKS one) is attractive though it is still qualitatively wrong. The on-top pair density interpretation provides a good way to understand the spin-symmetry breaking in this case, but there is no quantitative benefit. The resulting energetics from the UKS spin-densities are far from chemical accuracy (1 kcal/mol) as it is clear in Tables 1 and 2. In particular, due to the spatial symmetry breaking in each spin density and charge density, the KS-DFT XC energy evaluated with these spin densities are inadequate to reliably estimate the triplet-singlet gaps. In contrast to this, we observe a qualitatively correct charge density (and also a spin density which is a half of the charge density) with the cRKS determinant. At the same time, the resulting cRKS triplet-singlet gaps are improved as it is shown in Table 1.

Conclusions In this Letter, we showed that it is possible to access a different class of (both charge and spin) densities that is not possible to obtain within either RKS or UKS. This is achieved by using cRKS to obtain the density. Although the cRKS determinant breaks time-reversal symmetry and complex symmetry, the resulting densities do not exhibit this symmetry breaking. Furthermore, based on the triplet-singlet test set (TS12), we showed that the cRKS charge densities follow the point group symmetry while RKS does not. This allows cRKS to improve the quantitative accuracy of some XC functionals by a factor of 5. We also showed that the UKS charge density is qualitatively incorrect and the resulting energies are far from chemical accuracy due to spin-contamination. Lastly, we note that modern mGGA functionals (SCAN, MN15-L, and B97M-V) do not show any significant improvements even when correct densities from cRKS are used. Even with the exact XC functional, one needs a cRKS determinant to obtain qualitatively correct charge and spin densities for the systems considered here. This Letter suggests that these modern functionals might lack some aspects of the exact XC functional. We hope that our study provides a new class of data that can be used to assess, and possibly inform the design of new approximate XC functionals. The key benefit is a route to escape the symmetry dilemma whenever complex-polarization is relevant.

Acknowledgement This research was supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We thank Fionn Malone, Yuezhi Mao, and Narbe Mardirossian for helpful comments. The raw data for TS12 (Tables 1, 2, 3) is available in the Supplemental Material.