Two single-reference approaches to singlet biradicaloid problems: Complex, restricted orbitals and approximate spin-projection combined with regularized orbital-optimized Møller-Plesset perturbation theory

Strong correlation is usually associated with multiple open-shell electrons that are antiferromagnetically coupled into a low-spin state.^1–3 For instance, molecular magnets with multiple metal centers,^4–8 noninnocent ligands,^9,10 metalloenzymes,^11–13 and oligoacenes^14–17 exhibit strong correlation. Of such cases, singlet biradicaloids exhibit the simplest form of strong correlation.^18–24 As this is usually outside the scope of single-reference electronic structure methods, it is common to employ multiconfigurational methods.^25–27 A brute-force approach to treat this strong correlation is a complete active space self-consistent field (CASSCF) with an active space of two electrons in two orbitals (2e, 2o). However, CASSCF does not incorporate electron correlation outside the active space so subsequent dynamic correlation treatments^28–30 are necessary for quantitatively correct answers. A related single reference approach is to start from the triplet single determinant (i.e., M_S = 1) and flip a spin to access the M_S = 0 manifold, either using configuration interaction (CI)^5,31–34 or coupled-cluster (CC) via the equation of motion approach.^19,35,36

Alternatively, one could try to treat such systems using single-reference methods with the help of essential symmetry breaking. It is essential in the sense that the qualitative character of a single-determinant wavefunction is fundamentally wrong without essential breaking. A majority of essential symmetry breaking is spin-restricted (R) to spin-unrestricted (U) symmetry breaking, namely, spin-polarization. In the case of singlet biradicaloids, such essential symmetry breaking can be combined with Yamaguchi’s approximate spin-projection (AP) to produce spin-pure energies.^37–44 The applicability of AP is dependent on whether the underlying wavefunction contains only one contaminant. It is an exact projection only if there is one single contaminant. This sets a limit to ⟨$Ŝ$²⟩ of broken-symmetry M_S = 0 solutions to be effective for AP: 0.0 ≤ ⟨$Ŝ$²⟩ ≤ 2.0.

UHF is heavily spin-contaminated for most biradicaloids. For instance, this was observed by us in the heptazethrene dimer (HZD) where broken-symmetry UHF yields ⟨$Ŝ$²⟩ = 6.3 in the cc-pVDZ basis set.⁴⁵ The subsequent correlation treatment based on these UHF solutions via second-order Møller-Plesset (MP2) perturbation theory is not effective in removing such heavy spin contamination. It is possible to employ orbital optimized MP2 (OOMP2) as an attempt to produce a reference determinant with only essential symmetry breaking (i.e., ⟨$Ŝ$²⟩ ≈ 1.0). However, it is likely that OOMP2 produces a divergent solution or a restricted solution that is unphysically low in energy if not divergent.^46,47 As a solution to this problem, we employed regularized OOMP2 (κ-OOMP2) to treat HZD.⁴⁵ In contrast to our previous δ-OOMP2 (regularized with a constant level-shift), κ-OOMP2 determines the strength of regularization of individual correlation energy contributions depending on the orbital energy gap associated with them. κ-OOMP2, in turn, achieves both the recovery of Coulson-Fischer points⁴⁸ and favorable thermochemistry performance, which was found to be challenging for δ-OOMP2 to achieve.⁴⁷ Returning to the HZD example, κ-OOMP2 with unrestricted orbitals (κ-UOOMP2) produces ⟨$Ŝ$²⟩ = 1.2 which is well-suited for subsequent AP treatment. Generally speaking, κ-UOOMP2 with AP (AP+κ-UOOMP2) is a simple and robust way to treat biradicaloids, which captures both the static and dynamic correlation. We will further highlight this particular combination of AP and κ-UOOMP2 later in this work.

A rather rarer class of essential symmetry breaking, which is another focus of this work, is real, R to complex, R (cR) symmetry breaking. This is referred to as “complex-polarization” in this work. Complex polarization was known for many years in the context of some strongly correlated molecules such as O₂ (¹Δ_g).^49–57 Our group established its connection to generalized valence bond perfect pairing (GVB-PP)⁵⁸ using the complex pairing theorem.⁵⁹ When such solutions exist, complex restricted Hartree-Fock (cRHF) can indeed capture some aspects of GVB-PP and behaves qualitatively better than RHF. It was shown that the subsequent correlation treatment, cRMP2, yields quantitatively more accurate results than RMP2 for the systems examined in Ref. 59. Moreover, cRMP2 outperformed UMP2 especially when there is a strong mixing between singlet and triplet states.

Our recent work illustrated a way to obtain such essential symmetry breaking with κ-OOMP2.⁶⁰ Therein we discussed how to remove artificial spin-polarization using κ-OOMP2 with complex, generalized (cG) orbitals. It is artificial because orbital optimization in the presence of dynamic correlation such as MP2 (or other approaches for approximate Brückner orbitals) may remove such symmetry breaking. Artificial symmetry breaking occurs at the HF level not due to the lack of ability to describe strong correlation but because of the lack of dynamic correlation treatment. In Ref. 60, we show that it is possible to distinguish artificial and essential symmetry breaking based on κ-OOMP2. Interested readers are referred to Ref. 60, and we will further review some aspects of this relevant to this work in Sec. II B. In addition to essential spin-symmetry breaking, it is also possible to explore essential complex-polarization within the κ-OOMP2 method, which will combine the strengths of cRMP2 and κ-OOMP2. Namely, κ-cROOMP2 is able to describe multireference systems whenever complex-polarization is relevant.

For general biradicaloid systems, it is natural to consider AP and cR methods as simple single-reference alternatives to multireference and spin-flip methods. In particular, these are far simpler to implement than typical multireference second-order perturbation theory.^28,30 Compared to AP, cR methods offer more straightforward formalisms for response theory. For example, cRMP2 has the identical response theory formalism to that of usual MP2 and there is no need to derive additional terms. The analytic nuclear derivatives of AP methods have been derived and implemented at the mean-field level,^38,39,42,43 but there has been no study on response theory of correlated wavefunction methods with AP. While the formal and practical simplicity of cR methods is very desirable, its limited applicability due to the rareness of cR solutions makes it less appealing.

In this work, we will explore several biradicaloid systems that exhibit cRHF solutions and discuss the applicability of κ-cROOMP2 and AP+κ-UOOMP2. In particular, we will compare κ-cROOMP2 and AP+κ-UOOMP2 in these systems and discuss the similarities and differences between them. For simplicity, we will limit our discussion to HF, MP2, and κ-OOMP2 although other variants of MP2 and OOMP2, such as spin-component scaled methods,^61,62 can also be combined with cR orbitals or AP.

We will use i, j, k, l, … to index occupied orbitals, a, b, c, d, … to index virtual orbitals, and p, q, r, s, … to index either of those two.

The complex restricted Hartree-Fock (cRHF) energy is given by

where P is the one-particle reduced density matrix (1PDM), H₀ is the one-electron Hamiltonian, J and K are the Coulomb and exchange matrices, respectively, and E_nuc is the nuclear repulsion energy. In cRHF, we allow the molecular orbital (MO) coefficient matrix C to be complex and as a result P may become complex. As mentioned in Ref. 59, a cRHF solution is “fundamentally complex” if and only if the norm of the imaginary part of P is nonzero.

The use of complex restricted (cR) orbitals for multireference problems has been known for many years in electronic structure theory,^49–57 but they have been rarely employed in practice. The major reason for this underappreciation is due to the rareness of genuine cR solutions. Small et al. established the connection between cRHF and GVB-PP, and as a result, we have a better understanding of why cR solutions are rare and when to expect them.⁵⁹ 

Within a single pair of electrons, the R to cR instability is driven by the energy lowering due to a PP-like [or CAS(2,2)-like] configuration. However, a cRHF wavefunction also necessarily contains an open-shell singlet (OSS)-like configuration which is usually energetically high. The competition between the PP-like contribution (energy-lowering vs R) and the OSS-like contribution (energy-raising vs R) determines the R to cR instability. When the PP stabilization is greater than the OSS energetic cost, we observe the R to cR instability. This is, however, not very common to observe, and this explains the rareness of cRHF solutions. As we will see, some singlet biradicaloids exhibit complex-polarization and therefore cRHF can serve as a faithful starting point for subsequent correlation treatments. The relative energetics between the PP-like terms and OSS-like terms change in the presence of correlation treatment. Therefore, it is reasonable to expect that some cRHF solutions are artificial and they would lead to cRMP2 energies that are much higher than RMP2. We will encounter an example that demonstrates this later in this paper.

It is useful to run internal stability analysis to ensure the local stability of cRHF solutions. We provide the electronic Hessian of the energy expression in Eq. (1) in the supplementary material.

The MP2 energy expression with cRHF orbitals reads

where i and j are the occupied spatial orbitals, a and b are the unoccupied spatial orbitals, $ia|jb$ represents the two-electron four-center integrals, and the spin-adapted amplitudes τ are

$Δijab$ is a positive energy denominator defined as

where ϵ_p denotes canonical orbital energies. Orbital optimization of Eq. (2) yields orbital-optimized MP2 (OOMP2). As mentioned in Sec. I, OOMP2 has two major issues that limit its applicability. First, as we optimize orbitals in the presence of correlation energy, $Δijab$ can become very small and the resulting energy can become nonvariational and even approach divergence.⁴⁶ Second, as a result, OOMP2 may unphysically prefer restricted solutions and remove the Coulson-Fischer point.⁴⁷ Our group has developed a regularization scheme which fixes these two major issues in OOMP2.⁴⁵ 

The orbital-energy-dependent regularization introduced in Ref. 45 modifies the two-electron integrals in the correlation energy contribution in Eq. (2). The resulting κ-cRMP2 energy expression reads

where the regularized amplitudes are

Orbital optimizing Eq. (5) defines the κ-OOMP2 method (in this case, κ-cROOMP2). It is immediately obvious that the correlation energy can no longer diverge even when $Δijab=0$⁠. Based on carbon-carbon single, double, and triple bond breaking, we showed that the Coulson-Fischer point is recovered. Combining recovery of Coulson-Fischer points with reasonable performance for a thermochemistry benchmark, κ = 1.45 was recommended for chemical applications.⁴⁵ We also showed that κ ∈ [1.0, 2.0] (which comfortably includes κ = 1.45 in the middle) yields only essential symmetry breaking and can remove artificial HF symmetry breaking in fullerenes.⁶⁰ This is because κ-OOMP2 describes dynamic correlation, but regularization has removed the inaccurate description of static correlation present in conventional MP2.

Distinguishing artificial and essential symmetry breaking based on κ-OOMP2 may seem arbitrary. However, in Ref. 60, we compared this diagnosis of strong correlation with other approaches such as natural orbital occupation numbers and more sophisticated coupled-cluster methods. All these three independent probes suggested that C₆₀ is not strongly correlated and C₃₆ is strongly correlated. As such, κ-OOMP2 can reliably probe the underlying symmetry breaking and answer whether it is artificial (not strongly correlated) or essential (strongly correlated). κ-OOMP2 will be used to probe essential symmetry breaking and strong correlation in another fullerene C₃₀ below.

The implementation of κ-cROOMP2 was accomplished closely following the spin-orbital implementation described in Ref. 45. We apply the resolution-of-the-identity approximation to $ia|jb$⁠,

where P and Q are auxiliary basis indices and we define the expansion coefficients of an occupied-virtual product $|jb$ as

The spin-adapted two-particle density matrix (2PDM) consists of two parts: one is the usual MP2 2PDM contribution,

and another is the modification due to the regularizer,

Similarly, the spin-adapted 1PDM also consists of two parts: the first is the usual MP2 1PDM contributions,

and the second is the modification from the regularizer,

where the definition of ω_ij and ω_ab follows

and

These spin-adapted quantities can be used to produce appropriate orbital gradients for orbital optimization. Interested readers are referred to Ref. 45 for more technical details. In passing, we mention that Eqs. (13) and (14) were computed via a one-dimensional Legendre quadrature previously,⁴⁵ but in the pseudocanonical basis, it can be done analytically as shown above.

We apply the frozen-core approximation to the systems considered in this paper. This adds orbital rotation parameters between the frozen core and occupied orbitals to the orbital optimization problem. We present the pertinent orbital gradient equations and explain some numerical difficulties we encountered with this optimization problem in the supplementary material.

The approximate spin-projection method proposed by Yamaguchi³⁷ has been widely used in a wide variety of strong correlation problems.^37–44 Its working equation is very simple and it usually takes at most two separate single point calculations for two different M_S values to perform the projection. When projecting a triplet state out of an M_S = 0 broken symmetry solution, one can use the following equation which is derived using ⟨S²⟩:

where the spin-coupling coefficient α is

There are multiple ways to obtain E_S=1 and $⟨S2⟩S=1$⁠. The simplest way is to use a high spin M_S = 1 calculation to obtain E_S=1 at the same level of theory as E_BS. Therefore, we need a total of two unrestricted calculations, M_S = 0 and M_S = 1. Evidently, if the singlet is heavily spin-contaminated, the above spin-coupling equation is no longer valid. Furthermore, we need a nearly spin-pure value of ⟨S²⟩ for the M_S = 1 state. As we shall see later, κ-UOOMP2 can accomplish these objectives.

We will study multiple biradicaloid systems that have one pair of electrons that exhibit essential complex-polarization or spin-polarization. In other words, the singlet ground state of these systems involve a pair of open-shell electrons. Throughout the examples presented below, we will see how κ-cROOMP2 and/or AP+κ-UOOMP2 can be used for these singlet biradicaloids and also compare their strengths and weaknesses.

All calculations were performed with a development version of Q-Chem.⁶³ For κ-OOMP2 methods, we took a stable HF solution as an initial set of orbitals unless mentioned otherwise. All plots were generated with Matplotlib,⁶⁴ and all molecular figures were generated with Chemcraft.⁶⁵ 

We will consider triplet-singlet gaps (ΔE_T-S = E_S − E_T) of atoms and diatomics whose ground state is triplet. Systems with a triplet ground state are likely to have a (near) degeneracy between the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) so these are also likely to have near-degenerate OSS-like and PP-like configurations. Therefore, for these molecules, there is a good chance for essential complex polarization to occur.

We will compare HF, MP2, and κ-OOMP2 methods with different types of orbitals for treating the singlet ground state of the following molecules: C, NF, NH, NO⁻, O₂, O, PF, PH, S₂, S, Si, and SO. The reference triplet-singlet gaps as well as the equilibrium bond length of diatomics for each electronic state were taken from the NIST Chemistry WebBook.⁶⁶ Individual references for these experimental values and geometries are given in the supplementary material. We validated the experimental gaps against near-exact full configuration interaction calculations using the heat-bath algorithm developed by Holmes and co-workers.⁶⁷ The theoretical estimation lies within 1.0 kcal/mol of the experimental values and we provide these data in the supplementary material. This dataset will be referred to as the “TS12” set for the rest of this manuscript.

In benchmarking HF, MP2, and κ-OOMP2 methods, we employed the aug-cc-pVQZ basis set^68,69 along with its auxiliary basis set.⁷⁰ The frozen core approximation was used for all correlated wavefunction calculations. Unrestricted orbitals are used for the triplet state (M_S = 1).

For the molecules in the TS12 set, using real, restricted orbitals for the singlet ground state is fundamentally incorrect as it cannot capture the biradicaloid character of the singlet ground state. UHF orbitals are heavily spin-contaminated as the singlet ground state is a strong biradicaloid. This is well illustrated in Table I. The M_S = 0 states exhibit ⟨S²⟩ = 1.0 which indicates nearly perfect singlet biradicals. The M_S = 1 states are more or less spin-pure which validates the use of UHF orbitals for M_S = 1 states. Therefore, UHF and UMP2 are expected to perform very poorly on this test set. However, all these ⟨S²⟩ values are very well-suited for the AP approach. Therefore, one may expect that AP+UMP2 and AP+κ-UOOMP2 perform similarly well. We will see whether these predictions are indeed true in the TS12 set.

First, we discuss HF and MP2 with real, restricted (R) and real, unrestricted orbitals (U). The results of these methods are presented in Table II. Based on the mean-signed-deviation (MSD) of each method, it is evident that restricted orbitals overestimate the gap, whereas unrestricted orbitals underestimate the gap. This suggests that the singlet ground state of these molecules is too high in energy when described by R orbitals and too low in energy when described by U orbitals. This is expected for RHF because the closed-shell electronic structure produced by R orbitals should be less stable than an open-shell one. It is also expected for UHF, as the triplet ground state is lower in energy than the singlet ground state, triplet-singlet spin contamination lowers the energy of M_S = 0 unrestricted state. With the MP2 level of correlation, these failures of R and U orbitals do not disappear. RMP2 has an RMSD of 11.60 kcal/mol, and UMP2 has an RMSD of 12.42 kcal/mol.

How does κ-OOMP2 change this conclusion? As long as R or U orbitals are employed, very similar behavior is observed. As it is typical for RMP2 to overestimate correlation energies for singlet biradicaloids, we expect κ-ROOMP2 to produce larger triplet-singlet gaps than those of RMP2. This is mainly due to the regularization which is more effective on the singlet states here. Since none of the systems exhibit artificial spin-symmetry breaking (as presented in Table I), it is expected that κ-UOOMP2 methods do not significantly change the energetics of these systems.

In Table III, we see that the κ-ROOMP2 gaps are all greater than the RMP2 gaps in Table II, which confirms our prediction. For κ-UOOMP2, the gaps are all within 2 kcal/mol from those of UMP2, except O₂. In O₂, the difference between these two methods is 12.90 kcal/mol. This is due to the underlying artificial reflection spatial symmetry breaking in addition to the essential spin symmetry breaking in the UHF M_S = 0 solution. The artificial symmetry breaking is removed with κ-UOOMP2, while the essential one still persists.

We discuss whether these unrestricted states serve as reasonable bases to apply AP as well as whether cR orbitals can improve these catastrophic failures of HF, MP2, and κ-OOMP2 with R and U orbitals. The results of cR and AP methods are presented in Table IV. Neither cRHF nor AP+UHF produces satisfying results due to the lack of dynamic correlation. Moreover, cRHF and AP+UHF show significant differences in all molecules (the smallest difference is 4.21 kcal/mol and the largest one is 15.75 kcal/mol!). With MP2, cRMP2 is quite satisfying in that it has an RMSD of 1.64 kcal/mol with an MSD of −0.21 kcal/mol. The TS12 set can indeed be described properly with cR orbitals. On the other hand, the performance of AP+UMP2 is somewhat disappointing as it is poorer than cRMP2. In particular, an error of 29.34 kcal/mol in the case of O₂ is a striking outlier. This is due to spatial symmetry breaking in UHF M_S = 0 which cannot be fixed by UMP2 but can be fixed by κ-UOOMP2. Other than O₂, we observe a non-negligible difference (5.92 kcal/mol) in S₂ which is also caused by spatial symmetry breaking in the UHF solution. All the other molecules exhibit 2–3 kcal/mol differences between these two methods.

Orbital optimization in the presence of MP2 yields significantly better AP results, but κ-cROOMP2 produces slightly worse results than cRMP2. The slight degradation in the performance of cRMP2 in κ-cROOMP2 shows an interesting trend. All data points show larger triplet-singlet gaps with κ-cROOMP2 than with cRMP2. This indicates that there may be some overcorrelation problems with cRMP2 which is being regularized by κ-cROOMP2. Given the substantially better performance of κ-cROOMP2 compared to its R and U versions, this result is still very encouraging. Moreover, we emphasize that it is only κ-OOMP2 orbitals that yield quantitatively similar results between cR and AP approaches by harnessing only essential symmetry breaking.

For the rest of this work, we will further numerically show the quantitative similarity between AP+κ-UOOMP2 and κ-cROOMP2 beyond model systems.

There are not so many chemical systems for which cR methods can be a useful alternative to standard multireference methods. Any systems involving singlet oxygen [O₂(¹Δ_g)] are good candidates. In particular, singlet oxygen appears frequently in reactions in biological systems.⁷¹ An example that we will study here is the reaction between an amino acid, cysteine (Cys), and singlet oxygen. Cys is one of the five amino acids that are susceptible to singlet oxygen attack.⁷² Because of the multireference nature of singlet oxygen, studying reactivity of Cys is challenging for single-reference methods. As shown in Sec. III A, O₂(¹Δ_g) exhibits essential complex polarization. Therefore, this is an interesting example for comparing AP and cR approaches.

Lu et al. studied the reactivity of Cys ions with O₂(¹Δ_g) using Yamaguchi’s AP.⁷³ As mentioned earlier, in the case of singlet oxygen, the only spin contaminant is the triplet ground state. Therefore, AP is well-suited for this case. What Lu and co-workers found is that the reactivity of Cys ions with singlet oxygen is much smaller than that of neutral Cys. This was shown by a high activation barrier along a reactive pathway.

We will study a reaction between deprotonated Cys ([Cys-H]⁻) and singlet oxygen. Although there are multiple local minima geometries available, we investigated the lowest energy geometries from among those which Lu and co-workers reported. The molecular geometries of the precursor and transition state are shown in Fig. 1. Lu and co-workers optimized the geometries at the level of B3LYP with the 6-31+G(d) basis set with restricted orbitals.

The precursor in Fig. 1 has a substantial open-shell character due to the presence of singlet oxygen, but the transition state (TS) is a closed-shell molecule because of the formation of a persulfoxide. It is possible that the geometry optimization of the precursor [Fig. 1(a)] may produce a qualitatively wrong geometry when performed with restricted orbitals. We independently investigated this using unrestricted orbitals and could not find a local minimum similar to Fig. 1(a). A precise determination of the precursor geometry would be interesting to study in the future using cR orbitals or AP methods.

Nonetheless, for present purposes, we studied this system using the RB3LYP geometries from those of Lu and co-workers for single point cR and AP calculations. We employed the cc-pVTZ basis set⁶⁸ and the associated auxiliary basis set.⁷⁰ For the computational efficiency, the frozen core approximation was used for correlated wavefunction calculations. The goal of our study is to demonstrate the power of cR orbitals in comparison to AP methods (and conventional R and U orbitals) for the open-shell singlet precursor geometry.

We first discuss how different types of HF methods perform in predicting the reaction energy barrier [i.e., E(TS)—E(precursor)]. We compare the use of R, U, and cR orbitals for the precursor. The precursor RHF energy should be much higher than cRHF, whereas the UHF energy should be too low since the triplet contaminant is much more stable. Therefore, a back-of-the-envelope estimation for the relative energy barrier ordering is RHF < cRHF < UHF. This is indeed supported by numerical results presented in Table V.

The relative activation energy ordering will change based on the subsequent correlation treatment. For instance, UHF orbitals are heavily spin-contaminated so the subsequent UMP2 correlation energy will be underestimated, which then leads to a substantially smaller energy barrier than for UHF. Similarly, RHF orbitals should also lead to somewhat high energy when combined with MP2, which then yields a smaller energy barrier than that of cRMP2. Therefore, it is expected that the relative energy barrier ordering of MP2 methods is cR > R > U.

In Table VI, we observe the following trend instead: cR ≈ R > U. It is perhaps surprising that cR and R produce more or less the same energy barriers. There is quite strong complex polarization within a pair of electrons which led to a substantial energy lowering at the HF level. Evidently, despite the poor RHF reference, RMP2 recovers more correlation energy than cRMP2, perhaps because of overcorrelating small gap contributions.

Finally, we note that there is a significant energy difference between AP+UMP2 and cRMP2 similar to the O₂ triplet-singlet gap result observed in Sec. III A. However, this is mainly due to the qualitative difference between cR and U solutions. Spin-contamination often drives artifacts in the spin-density distribution which cannot be easily fixed by a posteriori spin projection methods. However, this can potentially be fixed by orbital optimizing in the presence of correlation as we shall see.

In Table VII, we present the activation barrier obtained using various types of κ-OOMP2 methods and two popular, combinatorially optimized density functional theory (DFT) methods (ωB97X-V⁷⁴ and ωB97M-V⁷⁵). First, we note that κ-cROOMP2 and AP+κ-UOOMP2 predict a barrier within 1 kcal/mol from each other. This is because the κ-UOOMP2 M_S = 0 state no longer has any artificial symmetry breaking and produces a solution with only essential spin-symmetry breaking. κ-ROOMP2 is similar to κ-OOMP2 with cR or AP despite the lack of open-shell character in the wavefunction. The R to cR instability at the κ-OOMP2 level causes an energy lowering of only about 2 kcal/mol. κ-UOOMP2 overestimates the gap by a factor of 2 compared to the corresponding AP results.

To see how well κ-cROOMP2 and AP+κ-UOOMP2 perform, we also compare this with two DFT methods. Without AP, both DFT methods with U orbitals predict the barrier too high. With AP, ωB97X-V predicts a barrier of 15.80 kcal/mol while ωB97M-V predicts a barrier of 11.03 kcal/mol. There is a quite significant functional dependence on the barrier height with the AP prescription (this may be related to the fact that ⟨$Ŝ$²⟩ cannot be rigorously evaluated: the expectation value of the KS determinant is used). A barrier height of about 10 kcal/mol was obtained with κ-cROOMP2, AP+κ-cROOMP2, and AP+ωB97M-V and an even higher height with AP+ωB97X-V. All of these suggest that the reactivity of [Cys-H]⁻ with O₂(¹Δ_g) is moderate at room temperature. In passing, we note that a higher level benchmark data would be desirable, using more sophisticated and computationally demanding methods such as the equation of motion spin-flip coupled-cluster with singles and doubles (EOM-SF-CCSD).³⁵ This will be interesting to study in the future.

Fullerenes are an interesting class of molecular clusters that are made solely of carbon atoms. They all form intriguing cage structures and often are stable enough to be experimentally synthesized. C₃₀ is one of the smaller fullerenes and it has been quite challenging to isolate C₃₀ experimentally due to its instability. It was pointed out in several experimental⁷⁶ and theoretical^77,78 studies that the highest symmetry structure, D_5h, is highly reactive. This particular molecular geometry is presented in Fig. 2.

The molecular geometries used in this work are optimized with unrestricted B97M-V⁷⁹ for each M_S state, with D_5h geometry within the cc-pVDZ basis set.⁶⁸ As this system is a biradicaloid, the geometry of the M_S = 0 state may require special care, but for simplicity, we employed unrestricted calculations. The ⟨S²⟩ value of each state with this particular functional is 1.020 and 2.013, respectively. We provide the geometries of C₃₀ used in this work in the supplementary material. Details about the geometries will not alter the qualitative conclusion we are drawing in this section as long as the underlying point group symmetry is D_5h.

Jiménez-Hoyos and co-workers reported the existence of complex generalized HF (cGHF) solutions for C₃₀(D_5h) and concluded that C₃₀ is a polyradicaloid based on the cGHF solution.⁸⁰ Due to its pronounced strong correlation, it is not surprising to observe symmetry breaking at the HF level. However, one may wonder if breaking every symmetry from RHF to cGHF is essential since UHF is sufficient for most singlet biradicaloid systems.

We have developed a computational strategy which can identify artificial symmetry breaking at the HF level using κ-OOMP2 with cG orbitals.⁶⁰ We scan over a range of κ values (i.e., the regularization strength) and compute the critical regularization strength, κ_c, to break/restore a given symmetry. Symmetry breaking with κ_c ∈ [0.0, 1.0] is categorized as artificial symmetry breaking, κ_c ∈ [1.0, 2.0] is essential symmetry breaking, and symmetry restoration for κ_c > 2.0 may be considered to be artificial restoration (i.e., too little symmetry breaking). The symmetry landscape of C₃₀ will help us to identify the character of essential symmetry breaking in this system.

We obtained the symmetry breaking landscape of C₃₀ within the 6-31G basis set⁸¹ along with the cc-pVDZ auxiliary basis set.⁷⁰ The frozen core approximation was used for computational efficiency. We focused on three symmetry breaking parameters: the spin expectation value ⟨$Ŝ$²⟩, the noncollinearity order parameter μ,⁸² and the fundamental complexification measure ξ.^59,60

In Fig. 3, we see that κ_c = 1.40 for μ, κ_c = 2.70 for ⟨S²⟩, and there is no obvious symmetry restoration for ξ. Compared to our previous work on characterizing other fullerenes such as C₆₀ and C₃₆, this landscape is more complex than the well-known biradicaloid C₃₆. Between κ = 1.40 and κ = 2.70, cU solutions are found. It is interesting that for κ > 2.70, cR solutions are most stable and there are no U or cU solutions. R solutions are commonly observed in the weak regularization regime, κ > 2.0, but cR solutions are quite unusual to observe. It turns out that this complex symmetry breaking in κ-cROOMP2 exists for all κ values as shown with the purple dashed line in Fig. 3. With κ-cGOOMP2 with κ < 2.70, κ-cROOMP2 solutions are higher in energy than other spin symmetry broken solutions. This is why these solutions are only observed with very weak regularization in the landscape. Based on these results, we conclude that the symmetry breaking of ⟨S²⟩, ξ, and μ is essential and this molecule is strongly correlated.

We also computed the triplet-singlet gap of C₃₀ using HF, MP2, and κ-OOMP2 methods with multiple types of orbitals. At κ = 1.45 with the cc-pVDZ basis set, we found a κ-cUOOMP2 solution (with ⟨S²⟩ = 2.0) when started from a cGHF solution. This triplet κ-cUOOMP2 solution was found to be almost exactly degenerate with a triplet κ-UOOMP2 solution. Therefore, for the remaining discussion, we employed κ-UOOMP2 for the triplet state.

The computed triplet-singlet gaps are presented in Table VIII, which are obtained with the cc-pVTZ basis set. Without correlation, RHF and cRHF predict very large gaps with a triplet ground state whereas UHF predicts a large gap with a singlet ground state. The MP2 correction on top of these reference states all prefers the singlet state with a significant spin gap. This is a qualitative failure of the MP2-level correlation treatment.

κ-OOMP2 provides a significant improvement over the MP2 results. κ-ROOMP2 predicts the sign of the gap correctly with a gap of 5.36 kcal/mol. κ-cROOMP2 yields a slightly smaller gap than κ-ROOMP2 and the energy lowering from complex polarization is only about 1.43 kcal/mol. κ-UOOMP2 yields almost a perfect open-shell solution (i.e., ⟨S²⟩ ≈ 1.0), so AP+κ-UOOMP2 is effective for this system. AP+κ-UOOMP2 predicts more or less the same gap as κ-cROOMP2 and the difference between two is only 0.18 kcal/mol. In terms of the triplet-singlet gap, all of the κ-OOMP2 approaches predict the biradicaloid character of C₃₀.

Although the triplet-singlet gap from the R methods is similar to the cR methods, the use of R orbitals breaks the spatial symmetry (D_5h) of C₃₀. This is evident when looking at the Mulliken population of individual carbon atoms. To illustrate this, we present the Mulliken population of the five carbon atoms in the top pentagon of C₃₀ in Fig. 2. Obviously, they are all equivalent due to the D_5h symmetry, but using restricted or unrestricted orbitals breaks this symmetry as shown in Table IX. Thus, geometry optimization with other methods than cR methods will likely break this spatial symmetry. This is not because of the Jahn-Teller distortion, but because of the artificial spatial symmetry breaking present at the electronic level. In passing, we mention that orbital-optimizing the spin-projected energy in Eq. (17) could potentially yield qualitatively better density than κ-UOOMP2.⁴¹ 

In summary, in this example, we showed that κ-cROOMP2 is better suited than AP+κ-UOOMP2 in the presence of high point group symmetry such as D_5h. Although they both yield similar energies, the underlying wavefunction breaks spatial symmetry if not treated with cR orbitals. This highlights the unique utility of electronic structure methods with cR orbitals whenever complex polarization is essential.

Although organic triplet biradicals are very rare to isolate due to their normally high reactivity, there have been some reports of synthesizing the stable ones.^24,84 Indeed, many stable singlet biradicaloids are stable because of some closed shell character.²² Since triplet biradicals lack in any closed shell character, it is difficult to observe stable ones. Gallagher and co-workers synthesized an organic biradical with a triplet ground state.⁸³ This biradical has quite robust stability compared to usual biradicals and survives at 140°C without significant decomposition. Experimentally, the triplet-singlet gap of this molecule was measured to be about 0.5 kcal/mol. However, such a small gap allows for a thermal mixture of singlet and triplet states as the temperature is raised to ambient conditions and above.

Gallagher and co-workers suggested a modification to this synthesized complex and hypothesized a triplet ground state, aiming for a larger triplet-singlet gap than 0.5 kcal/mol.⁸³ The structure of this proposed molecule is presented in Fig. 4. They supported their claim using UB3LYP/6-31G(d,p) calculations which yielded a gap of 3.5 kcal/mol for this newly suggested complex. Our goal is to confirm whether this hypothesis is correct using κ-cROOMP2 and/or AP+κ-UOOMP2. We studied this system within the cc-pVDZ basis set⁶⁸ and its auxiliary basis set⁷⁰ with the frozen-core approximation and the geometries were taken from Ref. 83 which were optimized at the UB3LYP/6-31G(d,p) level. This proposed system was recently synthesized and characterized with ΔE_TS ≥ 1.7 kcal/mol.⁸⁵ 

Unlike other examples presented above, there are no obvious symmetry constraints that give rise to a R to cR instability in this system. This is why it is interesting that there is a R to cR instability at the HF level (see Table X). However, this complex polarization turns out to be artificial and κ-OOMP2 with κ = 1.45 yields only a restricted solution. Therefore, in this case, κ-cROOMP2 is not applicable, whereas AP+κ-UOOMP2 is well-suited.

In Table X, the triplet-singlet gap of this system is presented. At the HF level, none of the orbital types predict small enough gaps to be considered to be a biradical. RHF and cRHF states are nearly degenerate, and thus, the complex polarization is not as strong as other examples presented before. UHF exhibits striking spin-symmetry breaking and predicts a much smaller spin gap than RHF and cRHF.

The MP2 treatment on top of these reference HF determinants does not improve these poor energetics. There is about a 20 kcal/mol energy difference between RMP2 and cRMP2, and RMP2 is lower in energy than cRMP2. This may indicate artificial complex polarization which indeed turns out to be the case in this system (vide infra). UMP2 removes a large portion of the spin contamination present at the HF level, but it still is heavily spin-contaminated. As a result, it predicts the sign of the gap incorrectly.

κ-UOOMP2 predicts a reasonably small gap with satisfying spin contamination for the singlet state (⟨S²⟩ ≈ 1.0) and almost no spin contamination for the triplet state. With the AP scheme, the gap is predicted to be 2.97 kcal/mol. This supports the original hypothesis⁸³ made by experimentalists that this system has a gap larger than 0.5 kcal/mol, with a triplet ground state. This is also in agreement with the recent experiment which studied this system.⁸⁵Finally, we note that κ-cROOMP2 collapses to a real, restricted solution and yields a substantially larger gap (36.23 kcal/mol) because this method does not adequately describe the strongly correlated singlet.

In summary, in this example, AP+κ-UOOMP2 successfully describes the biradicaloid character of the singlet state in the molecule, whereas κ-cROOMP2 cannot describe such character because there is no cR solution at the κ-OOMP2 level.

In this work, we examined two single-reference approaches based on regularized orbital-optimized Møller-Plesset perturbation theory (κ-OOMP2) that exploit essential symmetry breaking to describe singlet biradicaloids. Combined with Yamaguchi’s approximate projection (AP), unrestricted κ-OOMP2 (κ-UOOMP2) offers a way to access almost spin-pure singlet energies. Alternatively, complex, restricted κ-OOMP2 (κ-cROOMP2) can describe biradicaloid character if there is complex polarization. We compared these two methods over a variety of systems: a total of 12 triplet-singlet gaps in the TS12 set, the barrier height of a reaction between a cysteine ion and a singlet oxygen molecule, the C₃₀(D_5h) fullerene, and finally an organic biradical with a triplet ground state. We summarize the major conclusions from these numerical experiments as follows:

1. Without orbital optimization at the MP2 level, Hartree-Fock (HF) orbitals tend to exhibit artificial symmetry breaking in singlet biradicaloids. In the case of cRHF, this is sometimes reflected in spurious charge distribution of molecules, whereas it often manifests as heavy spin contamination (and commonly also spurious charge distribution) in UHF. In such cases, we recommend κ-OOMP2 which is an electronic structure tool that removes most artificial symmetry breaking and yields orbitals with only essential symmetry breaking.

2. κ-cROOMP2 is recommended whenever there is essential complex polarization. This is due to the fact that κ-UOOMP2 manifests not only spin-symmetry breaking but also spatial symmetry breaking which cannot be purified with the AP scheme.

3. When there is no essential complex polarization but only essential spin polarization, AP+κ-UOOMP2 is recommended. cR solutions are rare in nature and it is difficult to observe them with systems without point group symmetry. Therefore, the applicability of AP+κ-UOOMP2 is broader than that of κ-cROOMP2.

Strong correlation is a difficult problem to solve and there is no universal approach to it other than brute-force approaches such as complete active space methods.²⁷ However, at least for two-electron strong correlation problems studied here, either κ-cROOMP2 or AP+κ-UOOMP2 can be a single-reference electronic structure method that correctly describes strong correlation character. It will be interesting to apply these tools to a broader range of chemical systems along with more developments on their response theory such as excited states and analytic nuclear gradients in the future. The presented approaches, which use cR orbitals or AP, can be extended to higher order single-reference correlation methods such as coupled-cluster with singles and doubles (CCSD) and third-order Møller-Plesset (MP3) perturbation theory.

See supplementary material for the discussion of the frozen core and frozen virtual approximation in κ-OOMP2, the theoretical reference data of TS12 set, and the Cartesian coordinate of C₃₀.

This work was supported by a subcontract from MURI Grant No. W911 NF-14-1-0359. J.L. is thankful to Soojin Lee for consistent encouragement.