How to calculate S² value of a broken-symmetry wave function?
The spin operators for the $p^\textrm{th}$ MO are defined as follows:
\begin{align}
S_{x,p} &amp;= \frac{\hbar}{2} \sum_{\tau,\tau'} c_{p,\tau}^\dagger \hat{\sigma}^x_{\tau,\tau'} c_{p,\tau'}, \\
S_{y,p} &amp;= \frac{\hbar}{2} \sum_{\tau,\tau'} c_{p,\tau}^\dagger \hat{\sigma}^y_{\tau,\tau'} c_{p,\tau'}, \\
S_{z,p} &amp;= \frac{\hbar}{2} \sum_{\tau,\tau'} c_{p,\tau}^\dagger \hat{\sigma}^z_{\tau,\tau'} c_{p,\tau'}, \\
\end{align}
where $\sigma_{\tau,\tau'}^{x}$ denotes the matrix elements of the Pauli matrix $\sigma^x$ and so on. Any spin-observable can then be constructed from these basic blocks. For instance, $S^2$ operator takes the form
\begin{equation}
S^2 = \sum_{p \in MO} \frac{3}{4}\left( c_{p,\uparrow}^\dagger c_{p,\uparrow} +
c^\dagger_{p,\downarrow} c_{p,\downarrow} - 2c^\dagger_{p,\uparrow} c^\dagger_{p,
\downarrow} c_{p,\downarrow} c_{p,\uparrow} \right )
\end{equation}
All one need is to evaluate the expectation value of these operators over the wavefunction of interest.