Dealing with Spin

[PaulWAyers]

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.

[msricher]

S^2 can be computed in a similar way to how we compute RDMs, given the equation you provided. Loop over determinants i, loop over determinant excitations j from i, and add the resulting <i|j> ($c_ic_j$) values.

[PaulWAyers]

It might also be sensible to compute it by post-processing the 2DM.

[msricher]

Yes, that's doable too. The question is mainly do you care about the memory taken up by the 2RDM. If not, then it's much simpler to just post-process the 2RDM, and that can be done from Python.

…

On Wed, May 7, 2025, 09:24 Paul W. Ayers ***@***.***> wrote: *PaulWAyers* left a comment (theochem/PyCI#106) <#106 (comment)> It might also be sensible to compute it by post-processing the 2DM. — Reply to this email directly, view it on GitHub <#106 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADLODOW4PT2EWGUNPNET47325ICPBAVCNFSM6AAAAAB4RSC4L6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDQNJYGU4DKOBYGQ> . You are receiving this because you commented.Message ID: ***@***.***>

[PaulWAyers]

I think it's sensible to do it with the 2DM. Admittedly painful in terms of memory for big systems, but there are ways around that if push came to shove.

[PaulWAyers]

Thinking about this yesterday, I think it makes sense to it the way @msricher mentioned above. The reduced density matrices are "observables", and S^2, S^z, and other things we can compute from the solved wavefunction are also observables. So let's do it in PyCI, not as a separate post-processing of the RDM.

This is a really nice starter project for someone, and also something that @msricher @DanielCalero1 or @marco-2023 could just knock out in short order, as it resembles other things already in PyCI.

My main concern with the direct implementation is the way the implementation will change as one goes from restricted (obvious) to generalized (less obvious). I think one has to write the spin-operators in a restricted basis, then transform to the generalized basis, giving new (more complicated) operators there.