How do DFT software packages enforce spin multiplicity?

New to matter modeling

In KS-DFT software packages (I particularly have Gaussian in mind although I know they are ultimately not open-source), how is finding a ground state with a certain spin multiplicity enforced?

As I understand it, it's generally discouraged to try to use TD-DFT to compute T1 state energies. A better approach is to instead take the difference in SCF energies between optimized geometry calculations with spin multiplicities of 1 and 3 specified.

How do these spin multiplicities alter the calculational procedure?

• Numerical wavefunction solvers can do one of three things: One way is to simply assume that all electrons will be paired, so that the final result is assumed to be spin-less. This allows you to compute just half the electron orbitals. Another way is to assume that the spin has to be aligned with the z axis. You can have an imbalance of filling of the spin orbitals, but typically one would assume that the orbital and spin parts are essentially decoupled. Finally, a code can directly solve the Dirac equations with no extra assumptions. Commented May 2 at 17:59

For collinear spin, where the spin quantization axis is along $$z$$, one has \begin{align} S_z = \frac{{N_\alpha}-{N_\beta}}{2} \end{align}

Multiplicity $$(2S+1)$$ can then be specified given a particular number of up $$(\alpha)$$ and down $$(\beta)$$ electrons.

In unrestricted Hartree-Fock (UHF), the $$\alpha$$ and $$\beta$$ molecular orbitals $$(C)$$ aren't restricted to be the same spatially. This allows for greater variational flexibility in the self-consistent-field (SCF) procedure for finding open-shell solutions, since spin symmetry constraints are lifted.

In this case, we solve the Pople-Nesbet equations — again using different orbitals for different spins

\begin{align} F^{\alpha} C^{\alpha} = SC^{\alpha}{\epsilon}^{\alpha} \end{align} \begin{align} F^{\beta} C^{\beta} = SC^{\beta}{\epsilon}^{\beta} \end{align}

The UHF energy can be expressed as

\begin{align} \begin{split} E_{\mathrm{UHF}} = \sum\limits_{i}^{N_{occ}^{\alpha}}h_{ii} + \sum\limits_{\overline{i}}^{N_{occ}^{\beta}}h_{\overline{ii}} + \frac{1}{2}\sum\limits_{ij}^{N_{occ}^{\alpha}}\mathcal{J}_{ij}^{\alpha\alpha}+ \frac{1}{2}\sum\limits_{\overline{ij}}^{N_{occ}^{\beta}}\mathcal{J}^{\beta\beta}_{\overline{ij}}\\+ \frac{1}{2}\sum\limits_{i\overline{j}}^{N_{occ}^{\alpha\beta}}\mathcal{J}^{\alpha\beta}_{i\overline{j}}+ \frac{1}{2}\sum\limits_{\overline{i}j}^{N_{occ}^{\beta\alpha}}\mathcal{J}^{\beta\alpha}_{\overline{i}{j}}- \frac{1}{2}\sum\limits_{ij}^{N_{occ}^{\alpha}}\mathcal{K}_{ij}^{\alpha\alpha}- \frac{1}{2}\sum\limits_{\overline{ij}}^{N_{occ}^{\beta}}\mathcal{K}_{\overline{ij}}^{\beta\beta} \end{split} \end{align}

For a given geometry, the benefit of performing separate SCF calculations for each spin state is orbital optimization. However, significant spin polarization can lead to spin contamination, where your single determinant solution is no longer an eigenstate of the total spin-squared operator $$\mathcal S^2$$.

Deviations from the ideal spin state is observed in the following expectation value

\begin{align} \langle{\mathcal{S}^2}\rangle_{\mathrm{UHF}} = S_z(S_z + 1) + N_{\beta} - \sum\limits_{i\overline{j}}^{N_{occ}} {\left| \mathbb{S}_{i\overline{j}} \right|}^2 \end{align}

Note that when the $$\alpha\beta$$ molecular orbital overlap term $${\left| \mathbb{S}_{i\overline{j}} \right|}^2$$ does not cancel with $$N_{\beta}$$, one can encounter a situation where $$\langle{\mathcal S^2}\rangle_{\mathrm{UHF}} > S_z(S_z + 1)$$. The amount $$\langle{\mathcal{S}^2}\rangle_{\mathrm{UHF}}$$ can deviate from the exact value while still retaining a viable wavefunction is variable and depends on the how the reference determinant will be used. In general, one would like to have an $$\langle{\mathcal{S}^2}\rangle_{\mathrm{UHF}}$$ no greater than $$10-15\%$$ of the target spin.

Spin contamination can of course be ameliorated or avoided via spin recoupling by expanding the wavefunction in a basis of determinants as done in configuration interaction and related methods. Restricted open-shell SCF methods that conserve spin exist, however, specifics regarding their construction and use can be found elsewhere in the literature.

In standard implementations, the principles of the UHF method can be applied to unrestricted Kohn-Sham (KS-)DFT as well. Though, the concept of spin in KS-DFT is a bit ambiguous and there are plenty of discussions on its validity/practicality. Here are just a few papers on the topic:

https://doi.org/10.1080/00268976.2017.1393573

https://doi.org/10.1002/qua.21423

https://doi.org/10.1021/jp803064k

https://doi.org/10.1002/qua.560560414

For routine energy calculations or structure optimizations, you're most often fine with using your preferred density functional approximations as a starting point to model the chemical problem of interest.

We know that the exact solution of the Schrödinger equation is an eigenstate of $$\hat{\bf S}^2$$, since the total spin operator commutes with the molecular Hamiltonian, $$[\hat{\bf S}^2,\hat{H}]=0$$.

The multiplicity of an eigenstate of $$\hat{\bf S}^2$$ with eigenvalue $$S(S+1)$$ is $$2S+1$$, that is, there are $$2S+1$$ eigenfunctions of the Hamiltonian with the same energy, but with a different value of the $$z$$ component $$S_z$$.

In quantum chemical packages, one typically only chooses the value for the projection $$S_z$$; in density functional theory (DFT) and unrestricted Hartree-Fock this is achieved by occupying the corresponding number of spin-up and spin-down orbitals, and in CASSCF and FCI by restricting the list of configurations i.e. determinant strings with the analogous description.

The problem is that fixing $$S_z$$ in this way does not ensure that the calculation will converge onto the pursued state, since all states with $$S \ge |S_z|$$ will have roots with the given value of $$S_z$$. Depending on the employed orbitals (e.g. in some steps of a CASSCF orbital optimization), the pursued state with spin $$S$$ can appear as an excited state, instead.

For unrestricted Hartree-Fock, you of course only have a single determinant. Here, the issue is that it cannot be a spin eigenstate for $$S > 0$$, and the solution will often be symmetry broken. For DFT the situation is even more complicated: spin does not really have a clear definition in DFT; see e.g. this review by Jacob and Reiher.