Approximate Spin Projection for Broken-Symmetry Method and Its Application

Abstract

A broken-(spin) symmetry (BS) method is now widely used for systems that involve (quasi) degenerated frontier orbitals because of their lower cost of computation. The BS method splits up-spin and down-spin electrons into two different special orbitals, so that a singlet spin state of the degenerate system is expressed as a singlet biradical. In the BS solution, therefore, the spin symmetry is no longer retained. Due to such spin-symmetry breaking, the BS method often suffers from a serious problem called a spin contamination error, so that one must eliminate the error by some kind of projection method. An approximate spin projection (AP) method, which is one of the spin projection procedures, can eliminate the error from the BS solutions by assuming the Heisenberg model and can recover the spin symmetry. In this chapter, we illustrate a theoretical background of the BS and AP methods, followed by some examples of their applications, especially for calculations of the exchange interaction and for the geometry optimizations.

1. Introduction

For the past few decades, many reports about “polynuclear metal complexes” have been presented actively in the field of the coordination chemistry [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. Those systems usually have complicated electronic structures that are constructed by metal–metal (d-d) and metal–ligand (d-p) interactions. Those electronic structures caused by their unique molecular structures often bring many interesting and noble physical functionalities such as a magnetism [8, 9, 10, 11, 12, 13, 14, 15, 16, 17], a nonlinear optics [18], an electron conductivity [19], as well as their chemical functionalities, e.g., a catalyst and so on. For example, some three-dimensional (3D) metal complexes show interesting magnetic behaviors and are expected to be possible candidates for a single molecule magnet, a quantum dot, and so on [11, 12, 13, 14, 15, 16]. On the other hand, one-dimensional (1D) metal complexes are studied for the smallest electric wire, i.e., the nanowire [3, 4, 5, 6, 7, 17, 19]. In addition, it has been elucidated that the polynuclear metal complexes play an important role in the biosystems [20, 21, 22, 23, 24], e.g., Mn cluster [25, 26] in photosystem II and 4Fe-4S cluster [27, 28, 29, 30] in electron transfer proteins. In this way, the polynuclear metal complexes are widely noticed from a viewpoint of fundamental studies on their peculiar characters and of applications to materials. From those reasons, an elucidation of a relation among electronic structures, molecular structures, and physical properties is a quite important current subject.

Physical properties of molecules are sometimes discussed by using several parameters such as an exchange integrals (J_ab), on-site Coulomb repulsion, and transfer integrals of Heisenberg and Hubbard Hamiltonians, respectively, in material physics [31, 32, 33, 34, 35]. In recent years, on the other hand, direct predictions of such electronic structures, molecular structure, and physical properties of those metal complexes are fairly realized by the recent progress in computers and computational methods. In this sense, theoretical calculations are now one of the powerful tools for understanding of such systems. However, those systems are, in a sense, still challenging subjects because they are usually large and orbitally degenerated systems with localized electron spins (localized orbitals). The localized spins are caused by an electron correlation effect called a static (or a non-dynamical) correlation [36]. In addition, a dynamical correlation effect of core electrons also must be treated together with the static correlation in the case of the metal complexes. A treatment of both the static correlation and the dynamical correlation in large molecules is still a difficult task and a serious problem in this field. For those systems, a standard method for the static and dynamical correlation corrections is a complete active space (CAS) method [37, 38] or a multi-reference (MR) method [39] that considers all configuration interaction in active valence orbitals, together with the second-order perturbation correction, e.g., CASPT2 or MPMP2 methods. In addition to these methods, recently, other multi-configuration methods such as DDCI [40, 41, 42], CASDFT [43, 44, 45], MRCC [46, 47, 48], and DMRG-CT [49, 50, 51] methods are also proposed for the same purpose. These newer methods are developing and seem to be promising tools in terms of accuracy; however, real molecules such as polynuclear metal complexes are still too large to treat computationally with those methods at this state. An alternative way is a broken-symmetry (BS) method, which approximates the static correlation with a lower cost of computation [52, 53, 54, 55]. The BS method (or commonly known as an unrestricted (U) method) splits up and down spins (electrons) into two different spatial orbitals (it is sometimes called as different orbitals for different spins; DODS), so a singlet spin state of the orbitally degenerated system is expressed as a singlet biradical, namely, the BS singlet [55]. The BS method such as the unrestricted Hartree-Fock (UHF) and the unrestricted DFT (UDFT) methods are now widely used for the first principle calculations of such large degenerate systems. In this sense, the BS method seems to be the most possible quantum chemical approach for the polynuclear metal complexes, although it has a serious problem called the spin contamination error [56, 57, 58, 59, 60, 61, 62, 63, 64, 65]. Therefore one must eliminate the error by some kind of projection method. An approximate spin projection (AP) method, which is one of the spin projection procedures, can eliminate the error from the BS solutions and can recover the spin symmetry. In this chapter, we illustrate a theoretical background of the BS and AP methods, followed by some examples of their applications.

2. Theoretical background of AP method

In this section, the theoretical background of the BS and AP methods for the biradical systems is explained with the simplest two-spin model (e.g., a dissociated H₂) as illustrated in Figure 1(a).

Figure 1.

(a) Illustration of the two-spin states of the simplest two-spin model. (b) HOMO and LUMO of spin-adapted (SA) and BS methods. (c) Illustration of spin-symmetry recovery of BS method by AP method.

2.1. Broken-symmetry (BS) solution and approximate spin projection (AP) methods for the (two-spin) biradical state

In the BS method, the spin-polarized orbitals are obtained from HOMO-LUMO mixing [55, 56]. For example, HOMO orbitals for up-spin (${\psi }_{\text{HOMO}}$) and down-spin (${\overline{\psi }}_{\text{HOMO}}$) electrons of the simple H₂ molecule are expressed as follows (Figure 1(b)):

$${\psi }_{\text{HOMO}}^{\mathrm{BS}}=\text{cos}{\mathit{\theta \psi }}_{\text{HOMO}}+\text{sin}{\mathit{\theta \psi }}_{\text{LUMO}},$$E1

$${\overline{\psi }}_{\text{HOMO}}^{\mathrm{BS}}=\text{cos}{\mathit{\theta \psi }}_{\text{HOMO}}-\text{sin}{\mathit{\theta \psi }}_{\text{LUMO}},$$E2

where 0 ≤ θ ≤ 45° and ${\psi }_{\mathit{\text{HOMO}}}$ and ${\psi }_{\text{LUMO}}$ express HOMO and LUMO orbitals of spin-adapted (SA) (or spin-restricted (R)) calculations, respectively, as illustrated in Figure 1(b). And the wavefunction of the BS singlet (e.g., unrestricted Hartree-Fock (UHF)) becomes

$$\left|{\mathrm{\Psi }}_{\mathrm{BS}}^{\text{Singlet}}\right.⟩=\mathrm{co}{\mathrm{s}}^2\theta \left|{\psi }_{\text{HOMO}}{\overline{\psi }}_{\text{HOMO}}\right.⟩+\mathrm{si}{\mathrm{n}}^2\theta \left|{\psi }_{\text{LUMO}}{\overline{\psi }}_{\text{LUMO}}\right.⟩-\sqrt{2}\text{cos}\theta \text{sin}\theta \left|{\mathrm{\Psi }}^{\text{Triplet}}\right.⟩,$$E3

where ${\psi }_{\mathit{\text{HOMO}}}$ and ${\overline{\psi }}_{\mathit{\text{HOMO}}}$ express up- and down-spin electrons in orbital ${\psi }_{\mathit{\text{HOMO}}}$, respectively. If θ = 0, the BS wavefunction corresponds to the closed shell, i.e., SA wavefunctions, while if θ is not zero, one can have spin-polarized, i.e., BS wavefunctions. In the BS solution, ${\psi }_{\mathit{\text{HOMO}}}$ ≠ ${\overline{\psi }}_{\mathit{\text{HOMO}}}$ (Figure 1(b)), so that a spin symmetry is broken. In addition, it gives nonzero ${⟨{\widehat{\mathit{S}}}^2⟩}_{\mathrm{BS}}^{\text{Singlet}}$ value, and as described later, up- and down-spin densities appeared on the hydrogen atoms.

We often regard such spin densities as an existence of localized spins. An interaction between localized spins can be expressed by using Heisenberg Hamiltonian:

$$\widehat{\mathit{H}}=-2J_{\mathrm{ab}}{\widehat{\mathit{S}}}_{\mathbf{a}}\cdot {\widehat{\mathit{S}}}_{\mathbf{b}},$$E4

where ${\widehat{\mathit{S}}}_{\mathrm{a}}$ and ${\widehat{\mathit{S}}}_{\mathrm{b}}$ are spin operators for spin sites a and b, respectively, and J_ab is an effective exchange integral. Using a total spin operator of the system $\widehat{\mathit{S}}={\widehat{\mathit{S}}}_{\mathrm{a}}+{\widehat{\mathit{S}}}_{\mathrm{b}}$, Eq. (4) becomes

$$\widehat{\mathit{H}}=-2J_{\mathrm{ab}}\left(-{\widehat{\mathit{S}}}^{\mathit{2}}\mathit{+}{\widehat{\mathit{S}}}_{\mathbf{a}}^{\mathbf{2}}+{\widehat{\mathit{S}}}_{\mathbf{b}}^{\mathbf{2}}\right).$$E5

Operating Eq. (5) to Eq. (3), the singlet state energy in Heisenberg Hamiltonian ($E_{\mathrm{HH}}^{\text{Singlet}}$) is expressed as

$$E_{\mathrm{HH}}^{\text{Singlet}}=J_{\mathrm{ab}}\left(-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\text{Singlet}}+{⟨{\widehat{\mathit{S}}}_{\mathbf{a}}^{\mathit{2}}⟩}^{\text{Singlet}}+{⟨{\widehat{\mathit{S}}}_{\mathbf{b}}^{\mathit{2}}⟩}^{\text{Singlet}}\right).$$E6

Similarly, for triplet state

$$E_{\mathrm{HH}}^{\text{Triplet}}=J_{\mathrm{ab}}\left(-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\text{Triplet}}+{⟨{\widehat{\mathit{S}}}_{\mathbf{a}}^{\mathit{2}}⟩}^{\text{Triplet}}+{⟨{\widehat{\mathit{S}}}_{\mathbf{b}}^{\mathit{2}}⟩}^{\text{Triplet}}\right).$$E7

The energy difference between singlet ($E_{\mathrm{HH}}^{\text{Singlet}}$) and triplet ($E_{\mathrm{HH}}^{\text{Triplet}}$) states (S-T gap) within Heisenberg Hamiltonian should be equal to the S-T gap calculated by the difference in total energies of ab initio calculations (here we denote $E_{\mathrm{BS}}^{\text{Singlet}}$ and $E^{\text{Triplet}}$ for the BS singlet and triplet states, respectively). And if we can assume that spin densities of the BS singlet state on spin site i (i = a or b) are almost equal to ones of the triplet state, i.e.,${⟨{\widehat{\mathit{S}}}_i^2⟩}^{\text{Triplet}}\cong {⟨{\widehat{\mathit{S}}}_i^2⟩}^{\text{Singlet}}$, then J_ab can be derived as

$$J_{\mathrm{ab}}=\frac{E_{\mathrm{HH}}^{\text{Singlet}}-E_{\mathrm{HH}}^{\text{Triplet}}}{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\text{Triplet}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\text{Singlet}}}=\frac{E_{\mathrm{BS}}^{\text{Singlet}}-E^{\text{Triplet}}}{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\text{Triplet}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\text{Singlet}}}.$$E8

If the method is exact and the spin contamination error is not found in both singlet and triplet states (i.e., ${⟨{\widehat{\mathit{S}}}^2⟩}_{\text{Exact}}^{\text{Singlet}}=0$ and ${⟨{\widehat{\mathit{S}}}^2⟩}_{\text{Exact}}^{\text{Triplet}}=2$), the S-T gap between those states can be expressed as

$$E_{\text{Exact}}^{\text{Singlet}}-E_{\text{Exact}}^{\text{Triplet}}=2J_{\mathrm{ab}}.$$E9

The spin contamination in the triplet state is usually negligible (i.e., ${⟨{\widehat{\mathit{S}}}^2⟩}_{\text{Exact}}^{\text{Triplet}}\cong {⟨{\widehat{\mathit{S}}}^2⟩}^{\text{Triplet}}\cong 2$), and one must consider the error only in the BS singlet state, so the S-T gap becomes

$$E_{\mathrm{BS}}^{\text{Singlet}}-E^{\text{Triplet}}=2J_{\mathrm{ab}}-J_{\mathrm{ab}}{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\text{Singlet}}.$$E10

A second term in a right side of Eq. (10) indicates the spin contamination error in the S-T gap, and consequently, a second term in a denominator of Eq. (8) eliminates the spin contamination in the BS singlet solution. In this way, Eq. (8) gives approximately spin-projected (AP) J_ab values. Eq. (8) can be easily expanded into any spin dimers, namely, the lowest spin (LS) state and the highest spin (HS) state, e.g., singlet-quintet for S_a = S_b = 2/2 pairs, singlet-sextet for S_a = S_b = 3/2 pairs, and so on, as follows:

$$J_{\mathrm{ab}}=\frac{E_{\mathrm{BS}}^{\mathrm{LS}}-E^{\mathrm{HS}}}{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}}.$$E11

Eq. (11) is the so-called Yamaguchi equation to calculate J_ab values with the AP procedure, which is simply denoted by J_ab here. The calculated J_ab value can explain an interaction between two spins. If a sign of calculated J_ab value is positive, the HS, i.e., ferromagnetic coupling state, is stable, while if it is negative, the LS, i.e., antiferromagnetic coupling state is stable. Therefore, one can discuss the magnetic interactions in a given system.

2.2. Approximate spin projection for BS energy and energy derivatives

Because J_ab calculated by Eq. (11) is a value that the spin contamination error is approximately eliminated, it should be equal to J_ab value calculated by the approximately spin-projected LS energy ($E_{\mathrm{AP}}^{\mathrm{LS}}$) as

$$J_{\mathrm{ab}}=\frac{E_{\mathrm{BS}}^{\mathrm{LS}}-E^{\mathrm{HS}}}{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}}=\frac{E_{\mathrm{AP}}^{\mathrm{LS}}-E^{\mathrm{HS}}}{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\text{exact}}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\text{ecact}}^{\mathrm{LS}}}.$$E12

Here, we assume ${⟨{\widehat{\mathit{S}}}^2⟩}_{\text{Exact}}^{\mathrm{HS}}\cong {⟨{\widehat{\mathit{S}}}^2⟩}^{\mathrm{HS}}$; then one can obtain a spin-projected energy of the singlet state without the spin contamination error as follows [62, 63, 64, 65]:

$$E_{\mathrm{AP}}^{\mathrm{LS}}=\alpha E_{\mathrm{BS}}^{\mathrm{LS}}-\beta E^{\mathrm{HS}},$$E13

where

$$\alpha =\frac{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\text{exact}}^{\mathrm{LS}}}{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}}$$and $$\beta =\alpha -1$$E14

Then, we explain about derivatives of this spin-projected energy ($E_{\mathrm{AP}}^{\mathrm{LS}}$). In order to carry out the geometry optimization using the AP method, an energy gradient of $E_{\mathrm{AP}}^{\mathrm{LS}}$ is necessary. $E_{\mathrm{AP}}^{\mathrm{LS}}$ can be expanded by using Taylor expansion:

$$E_{\mathrm{AP}}^{\mathrm{LS}}\left({\mathit{R}}_{\mathit{AP}}^{\mathit{LS}}\right)=E_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)+{\mathit{X}}^T{\mathit{G}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)+\frac{1}{2}{\mathit{X}}^T{\mathit{F}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)\mathit{X},$$E15

where ${\mathit{G}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)$ and ${\mathit{F}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)$ are the first and second derivatives (i.e., gradient and Hessian) of $E_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)$, respectively [62, 63, 64, 65]; ${\mathit{R}}_{\mathrm{AP}}^{\mathrm{LS}}$ and R are a stationary point of $E_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)$ and a present position, respectively; and X is a position vector ($\mathit{X}={\mathit{R}}_{\mathrm{AP}}^{\mathrm{LS}}-\mathit{R}$). The stationary point ${\mathit{R}}_{\mathrm{AP}}^{\mathrm{LS}}$ is a position where ${\mathit{G}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)=0$; therefore one can obtain ${\mathit{R}}_{\mathrm{AP}}^{\mathrm{LS}}$ if ${\mathit{G}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)$ can be calculated. By differentiating $E_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)$ in Eq. (13), we obtain

$${\mathit{G}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)=\frac{\mathit{\partial }{E}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)}{\mathit{\partial }\mathit{R}}=\left\{\alpha \left(\mathit{R}\right){\mathit{G}}_{\mathrm{BS}}^{\mathrm{LS}}\left(\mathit{R}\right)-\beta \left(\mathit{R}\right){\mathit{G}}^{\mathrm{HS}}\left(\mathit{R}\right)\right\}+\frac{\mathit{\partial \alpha }\left(\mathit{R}\right)}{\mathit{\partial }\mathit{R}}\left\{E_{\mathrm{BS}}^{\mathrm{LS}}\left(\mathit{R}\right)-E^{\mathrm{HS}}\left(\mathit{R}\right)\right\},$$E16

where ${\mathit{G}}_{\mathrm{BS}}^{\mathrm{LS}}$ and ${\mathit{G}}^{\mathrm{HS}}$ are the first energy derivatives (energy gradients) of the BS and the HS states, respectively. As mentioned above, the spin contamination in the HS state is negligible, so that ${⟨{\widehat{\mathit{S}}}^2⟩}^{\mathrm{HS}}$ is usually a constant. Then $\partial \alpha \left(\mathit{R}\right)/\partial \mathit{R}$ can be written as

$$\frac{\mathit{\partial \alpha }\left(\mathit{R}\right)}{\mathit{\partial }\mathit{R}}=\frac{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\text{exact}}^{\mathrm{LS}}}{{\left({⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}\right)}^2}\frac{\mathit{\partial }{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}}{\mathit{\partial }\mathit{R}}.$$E17

By using Eqs. (16) and (17), the AP optimization can be carried out. In addition, one can also calculate the spin-projected Hessian (AP Hessian; ${\mathit{F}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)$ in Eq. (15)) as follows:

${\mathit{F}}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)=\frac{{\mathit{\partial }}^2{E}_{\mathrm{AP}}^{\mathrm{LS}}\left(\mathit{R}\right)}{{\mathit{\partial }}^2\mathit{R}}=\left\{\alpha \left(\mathit{R}\right){\mathit{F}}_{\mathrm{BS}}^{\mathrm{LS}}\left(\mathit{R}\right)-\beta \left(\mathit{R}\right){\mathit{F}}^{\mathrm{HS}}\left(\mathit{R}\right)\right\}$,

$$+2\frac{\mathit{\partial \alpha }\left(\mathit{R}\right)}{\mathit{\partial }\mathit{R}}\left\{{\mathit{G}}_{\mathrm{BS}}^{\mathrm{LS}}\left(\mathit{R}\right)-{\mathit{G}}^{\mathrm{HS}}\left(\mathit{R}\right)\right\}+\frac{{\mathit{\partial }}^2\alpha \left(\mathit{R}\right)}{{\mathit{\partial }}^2\mathit{R}}\left\{E_{\mathrm{BS}}^{\mathrm{LS}}\left(\mathit{R}\right)-E^{\mathrm{HS}}\left(\mathit{R}\right)\right\},$$E18

where ${\mathit{F}}_{\mathrm{BS}}^{\mathrm{LS}}$ and ${\mathit{F}}^{\mathrm{HS}}$ are the Hessians calculated by the BS and the HS states, respectively. And a second derivative of α can be expressed by

$$\frac{{\mathit{\partial }}^2\alpha \left(\mathit{R}\right)}{\mathit{\partial }{\mathit{R}}^2}=\frac{2\left({⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\text{exact}}^{\mathrm{LS}}\right)}{{\left({⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}\right)}^3}{\left(\frac{\mathit{\partial }{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}}{\mathit{\partial }\mathit{R}}\right)}^2+\frac{{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\text{exact}}^{\mathrm{LS}}}{{\left({⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}^{\mathrm{HS}}-{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}\right)}^2}\frac{\mathit{\partial }{⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{BS}}^{\mathrm{LS}}}{\mathit{\partial }\mathit{R}}.$$E19

By using Eqs. (18) and (19), the spin-projected vibrational frequencies are also calculated. The AP optimization can be carried out based on Eq. (16) with $\partial {⟨{\widehat{\mathit{S}}}^2⟩}_{\mathrm{BS}}^{\mathrm{LS}}/\partial \mathit{R}$ obtained by numerical fitting or analytical ways.

2.3. Relationship between the BS and projected wavefunctions

As well as a calculated energy and its derivatives, the BS wavefunction itself has also vital information. Here let us go back to Eq. (3). From the equation, an overlap between up-spin (so-called alpha) and down-spin (so-called beta) orbitals (T) becomes

$$T=⟨{\psi }_{\text{HOMO}}^{\mathrm{BS}}|{\overline{\psi }}_{\text{HOMO}}^{\mathrm{BS}}⟩=\mathrm{co}{\mathrm{s}}^2\theta -\mathrm{si}{\mathrm{n}}^2\theta =\text{cos}2\theta .$$E20

And because occupation number (n) of natural orbital (NO) for the corresponding orbital is expressed as $n=2{\text{cos}}^2\theta$, we get the relation:

$$T=\text{cos}2\theta =n-1$$E21

On the other hand, we can define projected wavefunction (PUHF) by eliminating triplet species from BS singlet wavefunction from Eq. (3) as follows:

$$\begin{array}{c}\left|{\mathrm{\Psi }}_{\text{PUHF}}^{\text{Singlet}}\right.⟩=\sqrt{\frac{2}{1+{\left(\text{cos}2\theta \right)}^2}}\left(\frac{1+\text{cos}2\theta }{2}\left|{\psi }_{\text{HOMO}}{\overline{\psi }}_{\text{HOMO}}\right.⟩-\frac{1-\text{cos}2\theta }{2}\left|{\psi }_{\text{LUMO}}{\overline{\psi }}_{\text{LUMO}}\right.⟩\right)\\ =\sqrt{\frac{2}{1+T^2}}\left(\frac{1+T}{2}\left|{\psi }_{\text{HOMO}}{\overline{\psi }}_{\text{HOMO}}\right.⟩-\frac{1-T}{2}\left|{\psi }_{\text{LUMO}}{\overline{\psi }}_{\text{LUMO}}\right.⟩\right).\end{array}$$E22

If we focus on the second term, which is related to double (two-electron) excitation, its weight (W_D) can be obtained from Eqs. (21) and (22) as follows:

$$W_D={\left\{\sqrt{\frac{2}{1+T^2}}\frac{1-T}{2}\right\}}^2=\frac{1}{2}\left\{1-\frac{2T}{1+T^2}\right\}$$E23

This is the weight of double excitation calculated by the BS wavefunction. By applying Eq. (21)–Eq. (23), the W_D is related to the occupation number of the corresponding NO as follows:

$$y=2W_D=\frac{n^2-4n+4}{n^2-2n+2}.$$E24

This y value is called an instability value of a chemical bond (or diradical character). In the case of the spin-restricted (or spin-adapted (SA)) calculations, the y value is zero. However if a couple of electrons tends to be separated and to be localized on each hydrogen atom, in other words the chemical bond becomes unstable with the strong static correlation effect, the y value becomes larger and finally becomes 1.0. So, the y value can be applied for the analyses of di- or polyradical species, and it is often useful to discuss the stability (or instability) of chemical bonds. The idea is also described by an effective bond order (b), which is defined by the difference in occupation numbers of occupied NO (n) and unoccupied NO (n^*):

$$b=\frac{n-n^{\ast }}{2}$$E25

Different from the y value, the b value becomes smaller when the chemical bond becomes unstable. If we define the effective bond order with the spin projection b(AP), it is related to the y value:

$$b\left(\mathrm{AP}\right)=1–y$$E26

Those indices show how the BS and AP wavefunctions are connected. In addition, one can utilize the indices to estimate the contribution of double excitation for very large systems that CAS and MR methods cannot be applied.

Finally, a relationship between the BS wavefunction and $⟨{\widehat{\mathit{S}}}^2⟩$ values are briefly explained. The $⟨{\widehat{\mathit{S}}}^2⟩$ values of the BS singlet states do not show the exact value by the spin contamination error. $⟨{\widehat{\mathit{S}}}^2⟩$ value of the SA calculation is.

$${⟨{\widehat{\mathit{S}}}^{\mathit{2}}⟩}_{\mathrm{SA}}=S\left(S+1\right),\text{where}\mathrm{S}={\mathrm{S}}_{\mathrm{a}}+{\mathrm{S}}_{\mathrm{b}}$$E27

However, in the case of the BS singlet state of H₂ molecule, it becomes

$${⟨{\widehat{S}}^2⟩}_{\mathit{BS}}={⟨{\widehat{S}}^2⟩}_{\mathit{\text{exact}}}+N^{\mathit{\text{down}}}-\sum _{\mathit{ij}}{T}_{\mathit{ij}}\cong 1-T$$E28

where N^down and T are number of down electrons and the overlap between spin-polarized up-spin and down-spin orbitals in Eq. (21). Therefore $⟨{\widehat{\mathit{S}}}^2⟩$ is also closely related to a degree of spin polarization. For the BS singlet state of the hydrogen molecule model, by substituting Eq. (21) into Eq. (28), we can obtain

$${⟨{\widehat{S}}^2⟩}_{\mathit{BS}}\cong 2-n$$E29

Here we explain another aspect of the spin projection method. As depicted in Figure 1(c), the BS wavefunction indicates only one spin-polarized configuration, e.g., BS1 in the figure. However, in order to obtain a pure singlet wavefunction, which satisfies the spin symmetry, the opposite spin-polarized state (BS2) must be included. The projection method can give a linear combination of the both BS states, and therefore it can give an appropriate quantum state for the singlet state.

3. Application of BS and AP methods to several biradical systems

3.1. Hydrogen molecule: comparison among SA, BS, and AP methods by simple biradical system

In this section, we briefly illustrate how the BS and AP methods approximate a dissociation of a hydrogen molecule. Figure 2(a) shows potential energy curves of Hartree-Fock and full CI methods. In the case of the spin-adapted (SA) HF, i.e., the spin-restricted (R) HF method, the curve does not converge to the dissociation limit. On the other hand, the BS HF, i.e., spin-unrestricted (U) HF calculation, successfully reproduces the dissociation limit of full CI method. This result indicates that the static correlation is included in the BS procedure. Around 1.2 Å, there is a bifurcation point between RHF and UHF methods. Within the closed shell (i.e., SA) region, where r_H-H < 1.2 Å, the UHF solution does not appear, and the singlet state is described by RHF (single slater determinant). In this region, the energy gap between full CI and RHF that is known as correlation energy indicates a necessity of the dynamical correlation correction as discussed later.

Figure 2.

(a) Calculated potential energy surface of H₂ molecule by spin-restricted (R), spin-unrestricted (U), and approximate spin-projected HF methods as well as full CI method. (b) Calculated $⟨{\widehat{\mathit{S}}}^2⟩$, occupation number (n), and 2–n values of H₂ molecule by UHF calculation. (c) a weight of double (two-electron) excitation (W_D) by double CI (CID) calculation and y/2 values in Eq. 24. (d) Calculated effective exchange integral (J) values of H₂ molecule with several H-H distances. For all calculations, 6-31G** basis set was used.

In order to elucidate how the double-excitation state is included in the BS solution, the occupation numbers of the highest occupied natural orbital (HONO) are plotted along the H-H distance in Figure 2(b). The figure indicates that the occupation number is 2.0 in the closed shell region, while it suddenly decreases at the bifurcation point. And it finally closes to 1.0 at the dissociation limit. In Figure 2(c), calculated y/2 values from the occupation numbers are compared with the weight of the double excitation (W_D) of CI double (CID) method. The figure indicates that the BS method approximates the bond dissociation by taking the double excitation into account. As frequently mentioned above, the BS wavefunction is not pure singlet state by the contamination of the triplet wavefunction. In Figure 2(b), $⟨{\widehat{\mathit{S}}}^2⟩$ values of the BS states are plotted. It suddenly increases at the bifurcation point and finally closes to the 1.0, which corresponds to occupation number n at the dissociation limit. And as mentioned above, $⟨{\widehat{\mathit{S}}}^2⟩$ and 2-n values are closely related.

Next, we illustrate results of calculated effective exchange integral (J_ab) values of the hydrogen molecule by Eq. (11). The calculated J values are shown in Figure 2(d). In a longer-distance region (r_H-H > 2.0 Å), the AP-UHF method reproduces the full CI result, indicating that the inclusion of double excitation state and elimination of the triplet state work well within the BS and AP framework. On the other hand, in a shorter-region (r_H-H < 1.2 Å), a hybrid DFT (B3LYP) method reproduces the full CI curve. In the region, the dynamical correlation that the RHF method cannot include is a dominant. Therefore the dynamical correlation must be compensated by other approaches, such as MP, CC, and DFT methods. The hybrid DFT methods are effective way in terms of the computational costs; however, one must be careful in a ratio of the HF exchange. It is reported that a larger HF exchange ratio is preferable in the intermediate region as well as the dissociation limit [69, 70].

3.2. Dichromium (II) complex: effectiveness of hybrid DFT method for calculation of J value

Next, the BS and AP methods are applied for Cr₂(O₂CCH₃)₄(OH₂)₂ (1) complex [1] as illustrated in Figure 3(a). This complex involves a quadruple Cr(II)-Cr(II) bond (σ, π _//, π_⊥, and δ orbitals). Due to the strong static correlation effect, it requires the multi-reference approach. Within the BS procedure, as a consequence, the electronic structure of the complex is expressed by the spin localization on each Cr(II) ions. First, let us examine the nature of the metal–metal bond between Cr(II) ions. For the purpose, natural orbitals and their occupation numbers are obtained from the BS wavefunctions using an experimental geometry.

Figure 3.

(a) Illustration of Cr₂(O₂CCH₃)₄(OH₂)₂ (1) complex. (b) Calculated natural orbitals of complex 1 by UB3LYP/basis set I (basis set I: Cr, MIDI+p; others, 6-31G*).

As depicted in Figure 3(b), there are eight magnetic orbitals, i.e., bonding and antibonding σ, π _//, π_⊥, and δ orbitals that concern about the direct bond between Cr(II) ions. The NO analysis clarifies the nature of the Cr-Cr bond. If d-orbitals of two Cr(II) ions have sufficient overlap to form the stable covalent bond, the occupation numbers of each occupied orbital will be almost 2.0 (i.e., T is close to 1.0). As summarized in Table 1, however, those bonds show much smaller values. The occupation numbers of all of occupied σ, π, and δ orbitals are close to 1.0, indicating that electronic structure of the complex 1 is described by a spin-polarized spin structure like the biradical singlet state.

Orbital	Occupation number (n)	Overlap (T)
δ	1.148	0.148
πave2	1.242	0.242
σ	1.625	0.625

Table 1.

n and T values of complex 1 calculated by UB3LYP/basis set I1.

¹

Cr, MIDI+p, and others, 6-31G*

²

Averaged value of π_⊥ and π_//

By substituting the obtained energies and $⟨{\widehat{\mathit{S}}}^2⟩$ values into Eq. (11), J_ab values of the complex 1 are calculated as summarized in Table 2. In comparison with the experimental value, HF method underestimates the effective exchange interaction, while B3LYP method overestimates it. This result is quite similar to a tendency of the J_ab curve of H₂ molecule at the intermediation region in Figure 2(d). In that region, BH and HLYP method, which involves 50% HF exchange, gives better value in comparison with B3LYP. The results also suggest an importance of the effect on the ratio of the HF/DFT exchange for estimation of the effective exchange interaction [71, 72].

Method	Jab values
B3LYP	−734
BH and HLYP	−520
HF	−264
Expt	−490

Table 2.

Calculated J_ab values1 of complex 1 by several functional sets2.

¹

In cm⁻¹

²

Basis set I was used.

3.3. Singlet methylene molecule: Spin contamination error in optimized geometry by BS method and its elimination by AP method

Finally, we examine the spin contamination error in the optimized structure. Here we focus on a singlet methylene (CH₂). As illustrated in Figure 4(a), the methylene molecule has two valence orbitals (ψ₁ and ψ₂) and two spins in those orbitals. Those two orbitals are orthogonal and energetically quasi-degenerate each other. The ground state of the molecule is ³B₁ (triplet) state, and ¹A₁ (singlet) state is the first excited state. Components of the wavefunction of ¹A₁ state obtained by BS method as illustrated in Figure 4(b) have been graphically explained [36]. The spin-restricted method such as RHF considers only single component (the first term of Figure 4(b)) although the BS wavefunction involves three components as illustrated in Figure 4(b). The existence of the triplet component is the origin of the spin contamination error in this system.

Figure 4.

Illustrations of (a) a methylene molecule and (b) components of BS wavefunctions.

Both ¹A₁ and ³B₁ methylene molecules have bent structures, but the experimental data indicates a large structural difference between them. For example, as summarized in Table 3, experimental HCH angles (θ_HCH) of ¹A₁ and ³B₁ states are 102.4° and 134.0°, respectively [66, 67]. There have also been many reports of the SA results as summarized in Ref. [68]. On the other hand, the BS method is a convenient substitute for CI and CAS method, so here we examined the optimized geometry of the ¹A₁ methylene by SA and BS methods. In order to elucidate a dependency of the spin contamination error on the calculation methods, HF, configuration interaction method with all double substitutions (CID), coupled-cluster method with double substitutions (CCD), several levels of Møller-Plesset energy correction methods (MP2, MP3, and MP4(SDQ)), and a hybrid DFT (B3LYP) method are also examined. In the case of ¹A₁ state, all SA results are in good agreement with the experimental values; however, it is reported that energy gap between the singlet and triplet (S-T gap) value is too much underestimated [65]. On the other hand, all BS results overestimate the HCH angle. The difference in HCH angle between the BS values and experimental one is about 10–20°. The HCH angles of UCI and UCC methods are especially larger than MP and DFT methods, indicating that the post-HF methods even require some correction for such systems if the BS procedure is utilized. Therefore it is difficult to use the BS solution for ¹A₁ state without some corrections. On the other hand, by applying the AP method to the BS solution, the error is drastically improved, and the optimized structural parameters became in good agreement with experimental ones. The difference in the optimized θ_HCH values between the BS and the AP method, i.e., the spin contamination error in the optimized geometry, is about 10–20°. Those results strongly indicate that the spin contamination sometimes becomes a serious problem in the structural optimization of spin-polarized systems and the AP method can work well for its elimination. On the other hand, the optimized structure with the AP-UHF method almost corresponds to CASSCF(2,2) result. This means that the AP method approximates two-electron excitation in the (2,2) active space well. The θ_HCH values become smaller by including higher electron correlation with the larger CAS space such as CASSCF(6,6) or with the dynamical correlation correction such as MRMP2(2,2) and MRMP2(6,6). The result of the spin-projected MP4 (AP MP4(SDQ)) successfully reproduced the MRMP2(6,6) result, indicating that the AP method plus dynamical correlation correction is a promising approach.

Method	rCHa					θHCHb
	SA	BS	AP	(^3B1)		SA	BS	AP	(^3B1)
HF	1.097	1.083	1.098	1.071		103.1	115.5	102.9	130.7
CID	1.114	1.091	1.112	1.081		101.6	119.7	101.9	131.8
CCD	1.116	1.087	1.113	1.082		101.7	125.1	102.4	132.0
MP2	1.109	1.091	1.109	1.077		102.0	114.7	100.9	131.6
MP3	1.109	1.094	1.112	1.080		102.0	114.9	101.0	131.8
MP4(SDQ)	1.117	1.096	1.114	1.081		101.2	115.0	101.0	131.9
B3LYP	1.120	1.100	1.113	1.082		100.3	112.9	103.2	133.1
CASSCF(2,2)	1.097					102.9
CASSCF(6,6)	1.124					100.9
MRMP2(2,2)	1.109					102.0
MRMP2(6,6)	1.122					101.1
Expt.d	1.107			1.077		102.4			134.0

Table 3.

Optimized C-H bond lengths (r_CH)a and H-C-H angle (θ_HCH)b by SA, BS, and AP approaches with several methodsc.

^a

In Å

^b

In degree

^c

6-31G* basis set was used

^d

In Refs. [66, 67] for singlet and triplet states, respectively

By calculating Hessian, one can also obtain frequencies of the normal modes. In Table 4, the calculated frequencies of the normal mode singlet methylene are summarized. The significant difference between the BS and AP methods can be found in a bending mode. The BS result underestimates the binding mode frequency by the contamination of the triplet state. On the other hand, the AP result gives close to the experimental result of ¹A₁ species. In this way, the AP method is also effective for the normal mode analysis as well as the geometry optimization.

Method	θHCH°	Mode
		Symmetry	Bent	Antisymmetry
BS	114.1	3008	1069	3152
AP	104.5	2959	1252	3054
Expt.c (^1A1)	102.4	2806	1353	2865
(^3B1)	134.0	2992	963	3190

Table 4.

Calculated vibrational frequenciesa of singlet methylene by SA, BS, and AP approaches with several methodsb.

^a

In cm⁻¹

^b

B3LYP/6–31++G(2d,2p) was used

^c

In Refs. [66, 67] for singlet and triplet states, respectively.

4. Summary

In this chapter, we explain how the BS method breaks the spin symmetry and AP method recover it. In addition, we also demonstrate how those methods work the biradical systems. The theoretical studies of the large biradical and polyradical systems such as polynuclear metal complexes have been fairly realized by the BS HDFT methods in this decade. The BS method is quite powerful for the large degenerate systems, but one must be careful about the spin contamination error. Therefore the AP method would be important for those studies. For example, it is suggested that the spin contamination error misleads a reaction path that involves biradical transition states (TS) or intermediate state (IM) [73]. In addition, in the case of the more larger systems, e.g., metalloproteins, some kind of semiempirical approach combined with the AP hybrid DFT method by ONIOM method will be effective [74]. By using the method, the mechanisms of the long-distance electron transfers and so on will be elucidated. In such cases, one also must be careful about the parameter of the semiempirical approach to fit the spin-polarized systems. Recently, some improvements for PM6 method have been proposed [75, 76]. Because the PM6 calculation can be utilized for the outer region in ONIOM approach, therefore the AP method is also the effective method for the larger systems. In addition, the BS wavefunction can be applied for other molecular properties by combining with other theoretical procedures. For example, it was reported that the electron conductivity of spin-polarized systems could be simulated by using the BS wavefunction together with elastic Green’s function method [77], and some applications for one-dimensional complexes have reported [78, 79]. The results indicate that the BS wavefunctions can be applied for calculations of the physical properties of the strong electron correlation systems as well as their electronic structures. The spin-projected wavefunctions seem to be effective for such simulations of the physical properties. From those points of view, the BS and AP methods have a great potential to clarify chemical and physical phenomena that are still open questions.