Dealing with Spin

We need to be able to compute S^2 in PyCI, as determinant-based calculations are not naturally/automatically eigenfunctions of S^2.

We could use ModelHamiltonian to make the S^2 operator (it would be linked to issue #150 there) or we could just code it directly here, which might be more direct. The easiest way to do it, I think, is to use the formula in terms of the raising and lowering operators:

$$
\hat{S}^2 = \tfrac{1}{2} (\hat{S}{+} \hat{S}{-} + \hat{S}{-} \hat{S}{+} ) + S_z^2
$$

We could also implement a spin-purification/projection method to try to get a non-spin-pure eigenfunction a bit closer to a spin eigenstate.

Start by comparing the proposed ModelHamiltonian integration, linked issue #150 there, with implementing the S^2 operator directly in PyCI. Review the raising/lowering-operator formula in the issue and determine whether spin purification or projection is also in scope. Done means PyCI can compute S^2 and the agreed spin-purification behavior is defined and implemented.