Does computational chemistry get the molecular orbitals of dioxygen wrong?

The molecular orbitals of $\ce{O2}$ are typically shown as follows, with every orbital filled by a spin alpha (or up) and spin beta (or down) electron (except the homo levels which are singly occupied). 


[![O2 Molecular Orbitals](https://i.stack.imgur.com/ltbAl.jpg)](https://i.stack.imgur.com/ltbAl.jpg)


Calculations using for instance Hartree-Fock and density functional theory however produce alpha and beta orbitals with different energies - no degeneracy of alpha and beta orbitals. For instance, the B3LYP/6-311G\*\* HOMO/LUMO levels (taken from the [computational chemistry database](https://cccbdb.nist.gov/)) are as follows: alpha HOMO (-0.31186), alpha LUMO (0.19755), beta HOMO (-0.46536), beta LUMO (-0.11702). Clearly the HOMO levels are not degenerate, at odds with what we expect from the above molecular orbital diagram. Similar mismatch between the alpha and beta orbitals exists for orbitals below the HOMO, as well as using other levels of theory. 


Is computational chemistry getting the orbitals of $\ce{O2}$ wrong? Am I missing something? Conceptually I can't grasp that alpha and beta orbitals have different energies since at least for the case of filled orbitals when the electrons pair up I expect them to have the same energy. 




A triplet species, as your MO scheme shows, has two unpaired alpha spin electrons, which is why you don't have to look at the beta orbitals at all but you have to look, if the $\alpha$-HOMO and $\alpha$-HOMO-1 are equivalent or not.


After optimizing the geometry of triplet oxygen using the method you states, `B3LYP 6-311G**`, we end up with a bond length of $\pu{1.2058 Å}$ and are given the following energies:



```
 **********************************************************************

            Population analysis using the SCF density.

 **********************************************************************

 Orbital symmetries:
 Alpha Orbitals:
       Occupied  (SGG) (SGU) (SGG) (SGU) (PIU) (PIU) (SGG) (PIG)
                 (PIG)
       Virtual   (SGU) (PIU) (PIU) (SGU) (SGG) (SGG) (PIG) (PIG)
                 (SGU) (PIU) (PIU) (DLTG) (DLTG) (DLTU) (DLTU)
                 (SGG) (PIG) (PIG) (SGU) (PIU) (PIU) (SGG) (PIG)
                 (PIG) (SGU) (SGG) (SGU)
 Beta  Orbitals:
       Occupied  (SGG) (SGU) (SGG) (SGU) (SGG) (PIU) (PIU)
       Virtual   (PIG) (PIG) (SGU) (PIU) (PIU) (SGU) (SGG) (SGG)
                 (PIG) (PIG) (SGU) (PIU) (PIU) (DLTG) (DLTG) (SGG)
                 (DLTU) (DLTU) (PIG) (PIG) (SGU) (PIU) (PIU) (SGG)
                 (PIG) (PIG) (SGU) (SGG) (SGU)
 The electronic state is 3-SGG.
 Alpha  occ. eigenvalues --  -19.28009 -19.27992  -1.31462  -0.83929  -0.56862
 Alpha  occ. eigenvalues --   -0.56862  -0.55311  -0.31192  -0.31192   <------
 Alpha virt. eigenvalues --    0.19733   0.66242   0.66242   0.69937   0.71138
 Alpha virt. eigenvalues --    0.77370   0.78062   0.78062   1.23343   2.23120
 Alpha virt. eigenvalues --    2.23120   2.53472   2.53472   2.80188   2.80188
 Alpha virt. eigenvalues --    2.83167   3.40715   3.40715   3.75295   4.50288
 Alpha virt. eigenvalues --    4.50288   4.64528   4.80485   4.80485   5.97229
 Alpha virt. eigenvalues --   49.48209  49.67839
  Beta  occ. eigenvalues --  -19.24954 -19.24923  -1.25730  -0.75086  -0.51327
  Beta  occ. eigenvalues --   -0.46530  -0.46530
  Beta virt. eigenvalues --   -0.11708  -0.11708   0.24099   0.70896   0.70896
  Beta virt. eigenvalues --    0.71935   0.72389   0.79954   0.83108   0.83108
  Beta virt. eigenvalues --    1.27043   2.27711   2.27711   2.62147   2.62147
  Beta virt. eigenvalues --    2.87800   2.89802   2.89802   3.46348   3.46348
  Beta virt. eigenvalues --    3.78751   4.58054   4.58054   4.66543   4.88322
  Beta virt. eigenvalues --    4.88322   6.00448  49.51237  49.70844

```

You can clearly see, that the $\alpha$-HOMO and $\alpha$-HOMO-1 are degenerate with both an energy of $\pu{-0.31192 Eh}$. If you look at both MOs you can easily see why, because they are qualitatively the same combination of 2p AOs but with different orientation.


[![Both singly occupied alpha orbitals of triplet oxygen](https://i.stack.imgur.com/R19jJ.jpg)](https://i.stack.imgur.com/R19jJ.jpg)


As symmetry was turned on, I guess, this needs to be the case, but the result stays the same, if it is turned off.


So, does computational chemistry get the molecular orbitals of dioxygen wrong? No, the MOs are right.




---


**But actually you have a different question**: Why are the energies of alpha and beta MOs in unrestricted open shell calculations not identical?


There are two methods to treat open shell systems. One is restricted open (RO\*) and the other is unrestricted (U\*). To not expand it too much, the main difference is, that as restricted closed shell methods, the restricted open shell methods force the alpha and beta orbitals to share the same MOs. Then, of course, they also would have the same energies.


But the energies that you show are from an unsrestricted open shell calculation. If you use this method, alpha and beta MOs are not forced to be in the same spatial MOs. And this freedom leads to different energies alpha and beta MOs.


 6

The short answer is, yes, computational chemistry gets them wrong because we don't have a way to calculate orbital energies exactly. In the end, the root cause is the general problem with quantum many-fermion systems - we can't solve the Schrödinger equation exactly.


Let's have a look to different computational chemistry methods to understand why that's the case. In general, more sophisticated methods will reduce the energy gap between alpha and beta electrons, but as we will see, none will guarantee that it doesn't exist.


**Hartree-Fock methods**


Hartree-Fock methods are inherently approximate: when we solve a Hartree-Fock calculation, we are solving not a system of *N* electrons, but *N* coupled systems of electrons that only interact indirectly. Of the energy terms in the Hamiltonian, HF accounts for the kinetic energy of each electron, its interaction with all the nuclei, the coulombic repulsion between electrons, and the exchange energy with other electrons with the same spin (because they are fermions), *but* we ignore electronic correlation. In HF, electrons only interact with each other indirectly - by contributing to a background external potential in the resolution of each electron's particular Hamiltonian. This makes it necessary to solve the equations iteratively using a mean field approach - what we usually call "self-consistent field" or SCF. Correlation requires explicitly considering more than one electron at the same time, so it's not compatible with the core HF model.


But, even if it's approximate, shouldn't we still get the same, not exact, values for alpha and beta electrons? - sure we do, in systems where the error introduced by the HF assumptions is the same for alpha and beta electrons.


For instance, if you run an unrestricted HF calculation on singlet $\ce{O\_2}$, where you have the same number of alpha and beta electrons, you'll get exactly the same values for the alpha and beta orbitals (and these will be exactly the same as you'd get for a restricted HF calculation of the same system).


However, for triplet $\ce{O\_2}$, alpha and beta electrons are exposed to slightly different interactions - since correlation depends on the spin of the interacting electrons, ignoring it affects alpha and beta electrons differently.


**Post-Hartree-Fock methods**


Post-Hartree-Fock methods attempt to correct the HF limitations by first resolving a Hartree-Fock system and then adding a correction for electronic correlation. The problem, however is that the corrections aren't exact.


**Møller-Plesset perturbation methods** take the HF wavefunction and then add a numerical correction to the energy in a perturbative fashion - you apply the HF method to double-excited (MP2), triple-excited (MP3), and so on versions of the same ground state HF wavefunction and then introduce a correction to the energy of your orbitals (not to the wavefunctions themselves) based on that. It has several problems - you are only accounting for dynamic correlation (how electrons influence the movement of other electrons within your HF wavefunction), and, because excited estates are not variationally limited, there's no guarantee that you will not *overestimate* your correction and will end up with an energy below the exact energy. But even if those problems could be completely corrected, since part of the correlation is not accounted for, it will still give slightly different energies for alpha and beta orbitals in a system with $\mathrm{S > 0}$.


**Configuration interaction (CI)** methods apply the HF method not only to the ground state but also to a number of excited states, and then find the correct solution to the system variationally, as a linear combination of those configurations. This method is in principle exact - and it's variational, so we are guaranteed not to underestimate the energy - the problem is that it is a series expansion, in which we add a term for each possible configuration state function (CSF). Still, even a full CI taking into account all the possible CSFs won't account for *all* of the correlation - CI is only exact in the limit of a complete basis set. In practice, expanding the series is computationally very expensive, and it is truncated relatively early - full CI is only practical for the simplest of systems.


When it's applicable, however, difference in energy between alpha and beta electrons should be negligible.


**Coupled cluster**


Coupled cluster (CC) has a theoretical foundation that is distinct from HF and post-HF methods; essentially, it takes a wavefunction (usually a Slater determinant) and performs a number of excitations on it, producing a coupled cluster description of the system that, like full CI, is in principle exact. However, like CI, CC is a series expansion - you have terms for single, double, triple, etc. excitations, and we would require all possible excitations to be taken into account to fully consider electron correlation. As with CI, we usually truncate the series quite early - the usual "gold standard" of computational chemistry, against which all other methods have usually been benchmarked in the last couple of decades, is CCSD(T) - meaning that we consider the first and second term, and the third term only partially. This is usually considered to be as exact as it gets in computational chemistry, but again, it is prohibitive for all but the simplest of systems.


Like with full CI, CC differences between alpha and beta electrons will be negligible.


**Density functional theory**


I left DFT for the end because it's an entirely different beast, theoretically speaking, but in the end the problem comes down to the same issue: you don't account for correlation correctly, and that flaw impacts alpha and beta electrons differently for $\mathrm{S > 0}$.


Density functional theory, in principle, shouldn't have the same problems Hartree-Fock-related methods have, since we only with a single, global, three-dimensional density function $\mathrm{\rho(r)}$; and, as per the second Hohenberg-Kohn theorem, the ground state electron density is unique and variationally limited, so in principle we could keep applying the functional to different three-dimensional functions until we find the absolute minimum - and that's our exact ground state density.


There, are, however, two big problems.


First, in practical terms, we need a way to translate the electron density to wavefunctions so we can understand what's happening at the electron level. For instance, the very values for orbital energies that you mention for a DFT calculation are not inherently in the density function that DFT deals with; there's no way to distinguish the first from the second occupied orbitals in the density function. What we do is introduce a formalism to be able to translate the density function to orbital terms - that's the Kohn-Sham formalism, which boils down to treating the *global* particle density function, which distributes all *N* particles in space, as the sum of *N* single-particle functions, each corresponding to one electron, and which interact with each other only *indirectly* through a background external potential, which accounts for the coulomb, exchange and correlation interactions, so we update densities and potentials iteratively in another mean field, self-consistent approach. You'll notice that this is pretty much the same thing we did with the wavefunction in Hartree-Fock.


Two considerations, however. First, the Kohn-Sham formalism allows us to map the total density to orbital occupations and from there calculate a number of properties that we couldn't derive from the total density function, but it guarantees *only* that the sum of the single-electron densities is equal to the total density function for the ground state - not that the depiction of a *particular* orbital is accurate. So, values for specific molecular orbitals (such as orbital energy) should be taken with a pinch of salt. In fact, in practice, orbital properties from DFT are most accurate around the HOMO-LUMO gap, while values for core orbitals or virtual orbitals can diverge quite wildly from reference values. This is so because orbitals are constructed in a way that *best fits* density, and not one that *best fits* the wavefunction. You can find a lot of discussions about whether Kohn-Sham orbitals are physically meaningful and how they relate to the actual orbitals of the system.


As you can imagine, if you have important spin effects in your system, spin deviation in these properties are almost guaranteed.


But there's also an advantage over Hartree-Fock: in HF, we had to drop correlation because we couldn't incorporate it - it required treating more than one electron wavefunction explicitly. But, in DFT, wavefunctions are an afterthought - what we are minimising is the total density function, so there's nothing that prevents us from introducing correlation. We could even ditch the Kohn-Sham formalism and calculate correlation directly in terms of the total density function. So, aside from the problem of breaking down the total density function into a number of one-electron densities, we could in principle at least get an *exact* total alpha electron density and a total beta electron density, couldn't we?


Unfortunately, no, because we don't know how the functional for electrons looks like. And, if we move from wavefunctions to density, we also lose the ability to calculate exchange interactions exactly. The only system for which we can calculate these exactly is the free electron gas; so all the functionals we regularly use are proposals for the exchange or correlation terms, or both, that approximate the exchange-correlation interactions but that don't represent them exactly.


And as with all the other cases, if we are not capturing all the exchange and correlation effects, that will impact alpha and beta electrons differently in a system with $\mathrm{S > 0}$.


 2