# Several New Papers now Online

The end of the year is rapidly approaching, and I’d like to draw your attention to three papers of mine that have been recently published online.

Theoretical Characterization of Conduction-Band Electrons in Photodoped and Aluminum-Doped Zinc Oxide (AZO) Quantum Dots

This first paper was a fun collaboration between our theoretical work in the Li group and the experimental work done by Dan Gamelin at the UW. The paper looks at calculated UV-Vis spectra of photodoped and Aluminum-doped ZnO quantum dots. Both types of doping added “extra” electrons to the conduction band, and this resulted in some interesting properties that you can see in both theory and experiment.

What I was most surprised to see was that the “extra” electron simply acted like an electron moving in a spherical potential, just like a hydrogen atom! On the left we have the HOMO, which looks like an s-orbital. Then the two on the right are different LUMOs, one that looks like a p-orbital and another that looks like a d-orbital. We called the “super-orbitals”, since they look like atomic orbitals but exist in these tiny crystal structures. I think these are the first images of DFT-computed “superorbitals” ever published, though the “artificial atom” paradigm has been around for some time. It’s interesting to see complicated systems behaving like the simple quantum mechanical models we study when we first learn quantum mechanics!

J. J. Goings, A. Schimpf, J. W. May, R. Johns. D. R. Gamelin, X. Li, “ Theoretical Characterization of Conduction-Band Electrons in Photodoped and Aluminum-Doped Zinc Oxide (AZO) Quantum Dots,” J. Phys. Chem. C, 2014 , 118, 26584.

Assessment of Low-scaling Approximations to the Equation of Motion Coupled-Cluster Singles and Doubles Equations

This second paper details the implementation of several low-scaling methods for computing excited states (e.g. computing UV/Vis absorption). The idea was to take the highly accurate EOM-CCSD equations and simplify them using perturbation theory. That way, we could keep most of the accuracy of the method, while still having fast enough equations that large molecules could be studied. The resulting equations are closely related to other methods like CC2, CIS(D), and ADC(2), and I showed how they all relate to each other. We compared the performance to EOM-CCSD as well as experimental values. The results were promising, and the perturbative equations performed particularly well for Rydberg states. CC2, on the other hand, performs great for valence excitations.

J. J. Goings, M. Caricato, M. Frisch, X. Li, “ Assessment of Low-scaling Approximations to the Equation of Motion Coupled-Cluster Singles and Doubles Equations,” J. Chem. Phys., 2014 , 141, 164116.

Ab Initio Non-Relativistic Spin Dynamics

Finally, in this paper, we extended the generalized Hartree-Fock method to the time domain. In this proof-of-concept paper, we showed how a magnetic field can guide the spin dynamics of simple spin-frustrated systems. The key is reformulating the real-time time-dependent Hartree-Fock equations in the complex spinor basis. This allows the spin magnetizations on each atom to vary as a response to, say, an externally applied magnetic field. Here’s an example with a lithium trimer. Initially (left picture), all the spin magnetizations (that is, the spatial average of spin magnetic moment) point away from each other. Then, applying a static magnetic field (right picture) into the plane of the molecule causes each magnetization to precess at the Larmor frequency. The precession is shown in picoseconds by the colorization.

It’s really a beautiful idea in my opinion, and there is so much more to be done with it. For example, in our simple ab initio model, the spins only “talk” through Pauli repulsion, so they behave more or less independently. What would happen if we include spin-orbit coupling and other perturbations? That remains to be seen.

F. Ding, J. J. Goings, M. Frisch, X. Li, “ Ab Initio Non-Relativistic Spin Dynamics,” J. Chem. Phys., 2014 , 141, 214111.

# Broken Symmetries in Hartree-Fock

It’s been a busy summer, with traveling and giving talks and writing papers.

I’m back in Seattle again (yay!), and have been thinking a lot about how we get the simplest models of electrons in molecules, the Hartree-Fock equations.

While most theorists know about restricted Hartree Fock (RHF) and unrestricted Hartree-Fock (UHF), not as many know that these can be derived from simple (okay, sometimes not-so-simple) symmetry considerations.

RHF and UHF are actually subgroups of a large class of problems called the generalized Hartree Fock equations. GHF allows mixed spin descriptions (not just simple spin up or spin down), as well as complex wave functions.

These are actually quite important for describing frustrated spin systems, like this simple triangular system

The spins can’t all align favorably, so we call it “spin frustrated”. You can try this out on your own with, say, a triangle of hydrogens or lithiums or chromiums. UHF won’t give you the lowest energy solution. It will be GHF-unstable. Only GHF can give you the lowest energy solution, since it allows the MOs to take on both “spin up” and “spin down” character.

Weird, right?

If we insist that the GHF equations are invariant with respect to spin rotations along an axis (say, the z axis), we can get the UHF equations back.

If we insist that the GHF equations are invariant to time-reversal (doesn’t that sound cool?) as well as spin rotations along all axes (x,y,z), we get the real RHF equations.

Unfortunately, the literature on it is pretty sparse, and the content pretty dense. A paradox, if I’ve ever seen one. Thankfully, I found you can derive all the relationships using the tools of linear algebra.

The results are illuminating, and it’s fun to see how the different broken symmetry HF solutions relate to each other.

So here’s my take on re-deriving the different HF classes from generalized Hartree-Fock!

We want to classify broken symmetry wave function by investigating how they transform under the action of the invariance operators ${\hat{g}}$ constituting the symmetry group ${\mathcal{G}}$ of the spin-free electronic Hamiltonian ${\hat{H}}$ which is equivalent to

This simply means ${\hat{H}}$ is invariant to transformation by ${\hat{g}}$, whatever ${\hat{g}}$ may be. In general, when ${\vert \psi\rangle}$ is an eigenstate of ${\hat{H}}$, then ${\hat{g}\vert \psi\rangle}$ is also an eigenstate belonging to the same eigenvalue as ${\psi\rangle}$.
Now, exact eigenstates of ${\hat{H}}$ can be chosen to be simultaneous eigenstates of the various symmetries in ${\mathcal{G}}$. In contrast, for approximate variational wave functions (e.g. Hartree-Fock) the symmetry requirements represent additional constraints.

This is Löwdin’s ‘‘symmetry dilemma’’.

To state this dilemma again, we know the exact solution (lowest energy solution) will have certain symmetries, but if we include these symmetries in our approximate, variational Hamiltonian, we can only raise the energy and not lower it.

This means we can get closer to the exact solution by removing physical constraints! This is troublesome, but not a huge deal in practice.

It’s crucial to recognize that this problem arises only because we use an approximate independent-particle Hamiltonian.

There are several ways to break the single-determinental Hartree-Fock model into its various broken symmetry subgroups. However we do this, what we ultimately want to determine is the form of the operators ${\hat{A}}$ that are invariant with respect to similarity transformations by various subgroups ${\mathcal{H} \subset \mathcal{G}}$, i.e.

I find it easier to put ${\hat{A}}$ into a finite basis and treat this problem with the tools of linear algebra.

The basis for our independent particle model will be the spinor basis,

so any arbitrary spin function can be written

We will only worry about one body operators ${\hat{A}}$, which include the Fock operator as well as the unitary parameterization of the single determinant (c.f. Thouless representation).

Putting ${\hat{A}}$ into a spinor basis gives us

with

Or, in second quantization

The transformation properties of ${\hat{A}}$ under symmetry operations can be determined by examining the transformation of the ${U(2n)}$ generators ${\hat{E}^{i \sigma_1}_{j \sigma_2}}$, but it is much simpler to consider the transformations of the ${2\times2}$ (block) matrix ${\mathbf{A}}$, e.g.

Basically, we are looking for the constraints on ${\mathbf{A}}$ that make the above equation true, for any given symmetry operation.

The general invariance group of the spin-free electronic Hamiltonian involves the spin rotation group SU(2) and the time reversal group ${\mathcal{T}}$, i.e. SU(2) ${\otimes \mathcal{T}}$. SU(2) can be given in terms of spin operators,

All of this amounts to performing rotations in spin space. The form looks awfully similar to rotations in 3D space (in fact, the relationship between rotations in spin space and 3D space is much deeper than that). The time reversal group is given as

In general,

Where ${\hat{K}}$ is the complex conjugation operator.

There is nothing special about using ${\hat{S}_2}$, it’s just convention.

Time reversal doesn’t really affect time per se, but rather changes the direction of movement, be it linear momentum or angular momentum.

It’s an antiunitary operator (it must be, in fact) that consists of the spin-looking part (unitary) and then the complex conjugation operator (antiunitary), to be antiunitary overall. It affects electrons by flipping their spin and then taking the complex conjugate.

One final detail about the symmetry operators before we go on, and that is to observe that since SU(2) is a double cover, if we rotate in spin space by ${2\pi}$, we flip the sign of the wave function. If we rotate by another ${2\pi}$, we get the original state back. This is characteristic of fermions, of which electrons are a part.

Let us now consider the unitary transformations of the form ${U= exp(i\hat{B})}$ where ${\hat{B} = \hat{B}^{\dagger}}$, acting on some operator ${\hat{A}}$.This is valid for any unitary transformation that depends on a continuous parameter (which we have absorbed into ${\hat{B}}$). We insist

Now, by the Baker-Campbell-Hausdorff transformation, we can rewrite the right hand side as

Since by definition ${\hat{A} = e^{(-i\hat{B})} \hat{A} e^{(i\hat{B})}}$ It suffices to show that the constraints on ${\hat{A}}$ introduced by the symmetry operation are satisfied when ${\left[\hat{A},i\hat{B}\right] = 0}$.
This is an excellent result, because it means that we can just look at how the ${2\times2}$ matrix ${\mathbf{A}}$ transforms under the Pauli spin matrices which define the SU(2) spin rotation as well as the time-reversal operation (in the time-reversal case, we can use the BCH expansion of the unitary part ${\hat{S}_y}$ and then absorb the complex conjugation into the commutator expression).

We can show how RHF, UHF, and GHF all fall out of the different symmetry combinations (or lack thereof).

From here on out, I’ll drop the hats on my operators for clarity.

Let’s start by considering how many different types of symmetry subgroups we can have. A moments thought at the form of the spin rotation and time-reversal groups gives the following subgroups:

${\mathbf{I}}$. This means there is no symmetry.
${\Theta}$. This means we are just symmetric to time-reversal symmetry.
${K}$. Only symmetric to complex conjugation.
${S_z}$. Only symmetric to spin rotations about the z-axis.
${S_z}$,${\Theta}$. Symmetric to both spin rotation about z, as well as time-reversal symmetry.
${S_z}$, ${K}$. Symmetric to spin rotation about z, as well as complex conjugation.
${S^2}$. Symmetric to all spin rotations.
${S^2}$, ${\Theta}$. Symmetric to all spin rotations and time reversal symmetries.

That’s it! That’s all the cases you can get.

Let’s go through them, case by case, using the complex Fock matrix as an example, i.e.

Where we solve, for each symmetry (or group of symmetries),

where ${\mathbf{h}}$ is the generator of the symmetry.

${\mathbf{I}}$. No symmetry.

In this case, our transformation on ${\mathbf{F}}$ is rather simple. It looks like

So we get no constraints. We can mix spin, as well as take on complex values. This is the structure of the complex generalized Hartree Fock (GHF) Fock matrix.

${K}$. Complex conjugation symmetry.

If the only symmetry that holds is complex conjugation, our transformation looks like

Note that ${K}$ is its own inverse. It also only acts to either the left or the right. The asterisk indicates complex conjugation (not and adjoint!).

The constraint we get here is that the values of the Fock matrix have to be identical on complex conjugation. Since this can only happen if the values are real, we get the real GHF Fock equations.
${\Theta}$. Time reversal symmetry.

Now we start to get slightly more complicated. Using the Pauli matrix

To represent the unitary ${S_2}$ operation gives us

We see that this really only introduces two constraints, so we choose to eliminate ${\mathbf{F}_{\beta\alpha}}$ and ${\mathbf{F}_{\beta\beta}}$. This gives the final result of paired GHF, or

${S_z}$. Rotation about spin z-axis.

Here we use the Pauli matrix

And show that

which is only satisfied if

This gives us the complex version of UHF. We see invariance with respect to ${S_z}$ results in two separate spin blocks, with no restriction on whether they take real or complex values, or the dimension of either spin block.

${S_z}$,${\Theta}$. Rotation about spin z-axis and time reversal.

We are now at the point where we examine the effect of invariance to multiple symmetry operations.

It might concern you that considering multiple symmetry operations means that the order in which we perform symmetry operations matters. In general symmetry operations do not commute, but we will see that for our purposes order really doesn’t matter. Because we insist on invariance, the multiple symmetries never actually act on each other, therefore we don’t need to consider the commutator between them.

Or, put a different way, because each symmetry operation returns the system to its original state, we can consider each operation separately.The system contains no memory of the previous symmetry operation.

We can show this another way using the BCH expansion. Consider two symmetry operations on ${F}$ parameterized by ${A}$ and ${B}$:

Which is true if and only if = [F, iA] = 0 Which decouples ${A}$ and ${B}$. Considering multiple symmetry operations only gives us more constraints, and order doesn’t matter. Let’s see it in action then for time reversal and z-axial spin symmetry. Using the results for time reversal symmetry, we have

Which means the off diagonals must go to zero, giving our final result of paired UHF,

${S_z}$,${K}$. Rotation about spin z-axis and complex conjugation.

We do a similar thing as above for rotation about spin z-axis and complex conjugation. This one is particularly easy to show, starting from the results of complex conjugation symmetry.

Since symmetry with respect to ${K}$ forces all matrix elements to be real, we just get the real version of the results of symmetry with respect to ${S_z}$ — we get the real UHF equations!

${S^2}$,${S_z}$. Rotation about all spin axes.

We finally move onto invariance with respect to rotations about all spin axes. Again, this is a little weird because these operations aren’t commutative, but we have already shown that insisting on invariance leads to a decoupling of the symmetry operations.

(I should mention that if you are still unsure of this, feel free to brute force through all orders of symmetry operations. You’ll see that it makes no difference.)

Two things worth mentioning here: first, technically ${S_z}$ is already a part of the total spin rotation group, so it’s a little weird to separate them into ${S^2}$ and ${S_z}$, but we understand this as drawing the distinction that you can be symmetric to ${S_z}$ but not ${S^2}$.

If you are invariant to ${S^2}$, though, you will be invariant to ${S_z}$.

Second, while ${S^2}$ does technically have an operator representation, it is not a symmetry operation. Think about the form of the operator ${S^2}$: it’s essentially the identity, right? So when we say invariant with respect to ${S^2}$, what we really mean is that we are invariant to the whole spin rotation group, ${S_x,S_y,S_z}$.

To show invariance with respect to the spin group, it suffices just to consider any two spin rotations, since each spin operator can be generated by the commutator of the other two.

You can show this using the Jacobi identity. Say we are looking the invariance of ${F}$ with respect to generators ${A}$, ${B}$, and ${C}$, and ${[A,B] = iC}$ (like our spin matrices do). We want to show

The Jacobi identity tells us

Now, by definition the first two terms are zero, and we can evaluate the commutator of ${[A,B] = iC}$, which means if ${[F,A] = [F,B] = 0}$, then it must follow that ${[F,C] = 0}$ as well (the imaginary in front doesn’t make a difference; expand to see).

That being said, we can evaluate invariance with respect to all spin axes by using the results of ${S_z}$ and applying the generator ${S_x}$ defined by Pauli matrix ${\sigma_2}$ to it, where

Applying this gives

Which means that ${\mathbf{F}_{\alpha\alpha} = \mathbf{F}_{\beta\beta}}$, or

Invariance in this symmetry group results in the complex RHF equations, where the alpha and beta spin blocks are equivalent. Thus orbitals are doubly occupied.

${S^2}$,${S_z}$,${\Theta}$,${K}$. All spin rotations and time reversal.

Given the results we just obtained above, and understanding that time reversal contains the ${S_y}$ operator, we only need to take the previous results and make them invariant to complex conjugation. This is very simple, and we see that

In other words, the real RHF equations, since invariance with respect to complex conjugation forces the elements to be real.

N.B. Most of this work was first done by Fukutome, though my notes are based off of Stuber and Paldus. The notes here agree with both authors, though the derivations here make a rather large departure from the literature. Thus, it is helpful to reference the following work:

Fukutome, Hideo. ‘‘Unrestricted Hartree-Fock theory and its applications to molecules and chemical reactions.” International Journal of Quantum Chemistry 20.5 (1981): 955-1065.

Stuber, J. L., and J. Paldus. ‘‘Symmetry breaking in the independent particle model.” Fundamental World of Quantum Chemistry, A Tribute Volume to the Memory of Per-Olov Löwdin 1 (2003): 67-139.

I’d strongly recommend looking at Stuber and Paldus’ work first. Their terminology is more consistent with mainstream electronic structure theory.

# Linear Response Coupled Cluster

Linear response coupled cluster (LR-CC) is a way of getting electronic excitation energies at a coupled cluster level of theory. It’s also a general way of getting many response properties, which basically answer the question of “what do molecules do when you poke them”?

Here is a simple derivation of LR-CC, which I like because it parallels the derivation of linear response time-dependent DFT and Hartree Fock. There are much more robust ways to get these results, for example, propagator approaches have been particularly successful when applied to coupled cluster response theory, but I think it can mask a lot of the physics if you aren’t careful.

Anyway, here it is.

The coupled cluster ground state equations are given by

Since this is, by definition, a stationary state the following is equivalent

Since

Now, perturbing the system with monochromatic light of frequency ${\omega}$ gives us the modified Hamiltonian (in the dipole approximation):

This induces a response to the wavefunction, which is parameterized by the cluster operator ${\hat{T}}$.

The extra terms are higher order responses in the cluster operator. I should also note that this means T and R are in the same excitation manifold: in other words, if T is a single excitation operator, so is R. Same for doubles, triples, etc. Anyway, we are only considering the linear response here. Furthermore, from here on out I won’t write the complex conjugate parts. Though they are still there (and I will indicate them by dots), ultimately we collect terms of only ${e^{-i \omega t}}$ to get our final equations. We could collect the conjugates as well, but would get the same equations in the end.

The time dependent Schrödinger equation is

Plugging in the above perturbed expressions yields

Since ${\hat{T^{(0)}}}$ and ${\hat{R}^{(1)}}$ commute, we can rewrite as

Expanding the exponential containing ${\hat{R}^{(1)}}$

Finally perform the derivation, making sure not to forget the time dependent phase

Collect terms with ${e^{-i\omega t}}$, cancel out phase factors

In linear response theory, we assume the perturbation is small and send ${\hat{f}^{(1)} \rightarrow 0}$.

Thus the eigenvalues are just the ground state coupled cluster energy plus an excitation energy.

As an aside, it is possible (in fact, this is how LR-CC is done in practice) to get rid of the ${E_{CC}^{(0)}}$ entirely, so that you solve for the excitation energies directly. How? The simplest way to think about it is that when you evaluate the left hand side, you find that it contains expressions for the ground state energy (times ${R}$), so if we leave out these terms when solving the equations, we can solve for the excitation energies directly. In other words we force the ground state energies on both sides to cancel by being careful how we set up the equation. And so we get an equation like so:

Here’s the more jargony answer: if you apply Wick’s theorem to the Hausdorff expansion of the similarity transformed coupled cluster Hamiltonian, you retain only terms that share a sum with each other. Diagrammatically, this means we keep only the connected diagrams. If you go ahead and derive all diagrams, connected and disconnected, you find that the disconnected diagrams correspond exactly to the terms you want to cancel on the right.

# Metal-Organic Frameworks Paper Now Online

I had my second first author paper published a few weeks back in the Journal of Physical Chemistry A, and wanted to share it here. It’s an investigation of how the tiny little hydrogen molecule binds to open metal sites in metal organic frameworks.

Metal organic frameworks are like molecular sponges. They are very porous and small gases seem drawn to them, so researchers want to use them to store hydrogen gas for fuel cells. Another interesting use is to use them for purifying gases! Imagine if we could store and condense carbon dioxide right out of the air! Lots of potential environmental applications — problem is, no one really knows how they work.

We studied the hydrogen gas binding to a series of different metals in these metal organic frameworks. Here is a picture of one of the systems we studied. The full metal organic framework is on the left. To simplify our method, we made ‘mimics’, meaning we studied just one part of a metal organic framework. That’s what you see in the center and on the right. Here is a close up of an interaction Hydrogen is kind of weird in that it has no charge and only two electrons. So why does bind to the metal centers? To understand why, we used a method called Symmetry Adapted Perturbation Theory (SAPT). This method was designed to study weak interactions between molecules. In addition to calculating the energy of the hydrogen sticking to the metal, it also decomposes the interaction into physical forces: like electrostatic, dispersion, induction, etc. It can be an extremely accurate method (we compared it to CCSD(T) energies and got pretty much the same results).

What we found might surprise you, as it did me! Hydrogen binds weakly (like ~ 3 kcal/mol weakly) and in many cases a good third of the interaction comes from dispersion forces. Now, if you remember from general chemistry, the more electrons you have, the greater the effects of dispersion become. Yet hydrogen — that little H2, with only two electrons — had an interaction largely mediated by dispersion (electrostatics was the other major component, which wasn’t too surprising).

The conclusion? Don’t ignore (or underestimate) your dispersion forces! This is very important to researchers who continue to study these metal organic frameworks, as they like to use density functional theory (DFT), which has some problems describing dispersion (if you aren’t careful).

You can read the whole paper here.

# Kinetic Balance

Kinetic balance was discovered when early attempts at relativistic SCF calculations failed to converge to a bound state. Often, the energy was too low because of the negative energy continuum. The crux of the problem was that the four components of the spinor in relativistic methods were each allowed to vary independently. In other words, scientists treated each component with its own, independent basis without regard to how the components depended on each other. This problem was eventually solved by paying attention to the non-relativistic limit of the kinetic energy, and noting that the small and large components of the four-spinor are not independent. There is, in fact, a coupling that we refer to as ‘‘kinetic balance’’. I’ll show you how it works.

First, as you may have guessed, we only need to consider the one-electron operator. This is the matrix form of the Dirac equation, and it contains (in addition to other terms) the contributions to the kinetic energy. Written as a pair of coupled equations, we have

Where ${\mathbf{V}^{LL} = \langle \chi^{L} \vert V \vert \chi^{L} \rangle}$, ${\mathbf{S}^{LL} = \langle \chi^{L} \vert \chi^{L} \rangle}$, and ${\mathbf{\Pi}^{LS} = \langle \chi^{L} \vert -i\hbar \mathbf{\sigma} \cdot \mathbf{\nabla} \vert \chi^{S} \rangle}$, and so on. These are the potential, overlap, and momentum terms of the Dirac equation in a basis ${\vert \chi \rangle}$.

Now, if we look at the second of the two paired equations, we note that for potentials of chemical interest (e.g. molecular or atomic), ${\mathbf{V}^{SS}}$ is negative definite. We also know that when we solve the equation, we are looking for and energy above the negative energy continuum, which is to say we want ${E \> -2mc^2}$. Since the overlap is positive definite, putting all of these constraints together mean that we have a nonsingular matrix (and therefore invertible!) matrix in the second of our coupled equations. We can rewrite the second equation as

Substituting this expression back into the first (top) equation yields

This form is very useful for analysis. Now, we are going to use the matrix relation

with ${\mathbf{A} = 2mc^2\mathbf{S}^{SS}}$ and ${\mathbf{B} = \mathbf{V}^{SS} - E\mathbf{S}^{SS}}$. This leads to the rather long expression

Why this ridiculous form? Look closely at each side and their dependence on the speed of light, ${c}$. The left hand side has the ${c^0}$ terms, and the right hand side has the ${c^{-2}}$ terms. Since the non relativistic limit is found when the speed of light is infinite (${c \rightarrow \infty}$), the whole right hand side goes to zero. This gives us

Now, if this is indeed the true non-relativistic limit, then we find that our kinetic energy term ${\mathbf{T}^{LL}}$ is given by

Or, more explicitly,

Where that inner part, ${\vert \chi_{\kappa}^{S} \langle \left[S^{SS}\right]_{\kappa \lambda}^{-1} \langle \chi_{\lambda}^{S} \vert }$, is an inner projection onto the small component basis space. Less formally, the small component is the ‘‘mathematical glue’’ that connects the two momentum operators. If the small component spans the same space as the momentum operators in the large component space, ${\{\sigma\cdot\nabla\chi_{\mu}^{L}\}}$, then that inner projection just becomes the identity. This means that the expression becomes

Which is the kinetic energy term in the non-relativistic formulation! (N.B. We used the relation ${(\mathbf{\sigma}\cdot\mathbf{\nabla})(\mathbf{\sigma}\cdot\mathbf{\nabla}) = \mathbf{\nabla}^2}$).

So, when we set up our relativistic calculations, as long as we have the constraint that

Then we will find that we recover the correct non-relativistic limit of our equations. The basis is called ‘‘kinetically balanced’’, and we won’t collapse to energies lower than ${-2mc^2}$.

A few stray observations before we finish. First, if we enforce the relation between small and large component basis functions, then we find that ${\mathbf{\Pi}^{SL} = \mathbf{S}^{SS}}$ and ${\mathbf{\Pi}^{LS} = 2m\mathbf{T}^{LL} = 2m\mathbf{T}}$. Second, this constraint actually maximizes kinetic energy, and any approximation that does not satisfy the kinetic balance condition will lower the energy. This was weird to me coming from the non relativistic Hartree-Fock background, where variationally, if you remove basis functions you raise the energy. The thing is that when doing relativistic calculations, you aren’t bounded from below like their non-relativistic counterparts. While you can get variational stability, you are actually doing an ‘‘excited state’’ calculation (I am using the term ‘‘excited state’’ very loosely). Kind of odd, but the negative energy continuum does exist, and was a big factor in the predicting the existence of antimatter.

Finally, modern methods of relativistic electronic structure theory make use of the kinetic balance between large and small component basis functions to eliminate the small component completely. These are called ‘‘Dirac-exact’’ methods. One such example is NESC, or Normalized Elimination of the Small Component. In addition to reproducing the Dirac equation exactly, they have numerous computational benefits, as well as easily allowing for most (if not all) non-relativistic correlated methods to be applied directly. Thus after doing an NESC calculation, you get relativistic ‘‘orbitals’’ which can immediately be used in, say, a coupled cluster calculation with no computational modification.

Comments, questions, or corrections? Let me know!