Fix the AngularMomentum (again) (#1292)

* fix: use the symmetric formula for S^2

Prior to this fix, the `AngularMomentum` operator was working on the
assumption that the number of alpha-spin particles would always exceed
the number of beta-spin particles (for non-singlet systems). However,
this is not guaranteed to be the case. This commit fixes this problem by
using a symmetric formula for S^2 which is the average of the two cases
(alpha- vs. beta-spin particles exceeding the other).

* feat: log a warning when the alpha-beta overlap is non-unitary

When dealing with active spaces obtained from unrestricted-spin
orbitals, the alpha-beta overlap matrix which is used to construct the
S^2 operator can become non-unitary (for example, when the active set of
alpha- and beta-spin orbitals do not span the same space). This can
result in an <S^2> value measured on the active space that is largely
different from (e.g.) the UHF <S^2> value. More importantly, when this
is the case, the difference between these two <S^2> values is "hidden"
in the inactive+inactive and inactive+active interactions. Especially
when spin contamination is present, this can lead to vastly unexpected
results and may indicate a poor choice of active space.

* Add reno

* test: add a regression test

* fix: update warning tolerance

* Linting and formatting

* docs: fix verbatim in RST

qiskit_nature/second_q/properties/angular_momentum.py

@@ -14,6 +14,7 @@
 
 from __future__ import annotations
 
+import logging
 from typing import Mapping
 
 import numpy as np
@@ -23,6 +24,8 @@
 
 from .s_operators import s_minus_operator, s_plus_operator, s_z_operator
 
+LOGGER = logging.getLogger(__name__)
+
 
 class AngularMomentum:
     r"""The AngularMomentum property.
@@ -31,7 +34,7 @@ class AngularMomentum:
 
     .. math::
 
-       S^2 = S^- S^+ + S^z (S^z + 1)
+       S^2 = (S^+ S^- + S^- S^+) / 2 + S^z S^z
 
     .. warning::
 
@@ -49,20 +52,53 @@ class AngularMomentum:
 
     Attributes:
         num_spatial_orbitals (int): the number of spatial orbitals.
-        overlap (np.ndarray | None): the overlap-matrix between the $\alpha$- and $\beta$-spin
-            orbitals. When this is `None`, the overlap-matrix is assumed to be identity.
     """
 
     def __init__(self, num_spatial_orbitals: int, overlap: np.ndarray | None = None) -> None:
         r"""
         Args:
             num_spatial_orbitals: the number of spatial orbitals in the system.
             overlap: the overlap-matrix between the $\alpha$- and $\beta$-spin orbitals. When this
-                is `None`, the overlap-matrix is assumed to be identity.
+                is ``None``, the overlap-matrix is assumed to be identity.
         """
         self.num_spatial_orbitals = num_spatial_orbitals
+        self._overlap: np.ndarray | None = None
         self.overlap = overlap
 
+    @property
+    def overlap(self) -> np.ndarray | None:
+        r"""The overlap-matrix between the $\alpha$- and $\beta$-spin orbitals.
+
+        When this is ``None``, the overlap-matrix is assumed to be identity.
+        """
+        return self._overlap
+
+    @overlap.setter
+    def overlap(self, overlap: np.ndarray | None) -> None:
+        self._overlap = overlap
+
+        if overlap is not None:
+            norb = self.num_spatial_orbitals
+            delta = np.eye(2 * norb)
+            delta[:norb, :norb] -= overlap.T @ overlap
+            delta[norb:, norb:] -= overlap @ overlap.T
+            summed = np.einsum("ij->", np.abs(delta))
+            if not np.isclose(summed, 0.0, atol=1e-6):
+                LOGGER.warning(
+                    "The provided alpha-beta overlap matrix is NOT unitary! This can happen when "
+                    "the alpha- and beta-spin orbitals do not span the same space. To provide an "
+                    "example of what this means, consider an active space chosen from unrestricted-"
+                    "spin orbitals. Computing <S^2> within this active space may not result in the "
+                    "same <S^2> value as obtained on the single-reference starting point. More "
+                    "importantly, this implies that the inactive subspace will account for the "
+                    "difference between these two <S^2> values, possibly resulting in significant "
+                    "spin contamination in both subspaces. You should verify whether this is "
+                    "intentional/acceptable or whether your choice of active space can be improved."
+                    " As a reference, here is the summed-absolute deviation of `S^T @ S` from the "
+                    "identity: %s",
+                    str(summed),
+                )
+
     def second_q_ops(self) -> Mapping[str, FermionicOp]:
         """Returns the second quantized angular momentum operator.
 
@@ -75,7 +111,8 @@ def second_q_ops(self) -> Mapping[str, FermionicOp]:
         overlap_ba = overlap_ab.T if overlap_ab is not None else None
         s_m = s_minus_operator(self.num_spatial_orbitals, overlap=overlap_ba)
 
-        op = s_m @ s_p + s_z @ (s_z + FermionicOp.one())
+        spm_smp = (s_p @ s_m + s_m @ s_p).normal_order()
+        op = 0.5 * spm_smp + s_z @ s_z
 
         return {self.__class__.__name__: op}
 

releasenotes/notes/fix-angular-momentum-2-cebabde075b61a4e.yaml

@@ -0,0 +1,10 @@
+---
+features:
+  - |
+    The :meth:`.AngularMomentum.overlap` property will now warn the user, when
+    this matrix is non-unitary.
+fixes:
+  - |
+    Fixes the :class:`.AngularMomentum` operator further to also support cases
+    where the number of beta-spin particles exceeds the number of alpha-spin
+    particles.

test/second_q/properties/test_angular_momentum.py

@@ -104,6 +104,33 @@ def test_with_overlap(self):
         result = estimate_observables(Estimator(), hf_state, qubit_op)
         self.assertAlmostEqual(result["AngularMomentum"][0], 0.29663167846210015)
 
+    def test_with_non_unitary_overlap(self):
+        """Tests the result with a non-unitary overlap.
+
+        This is a regression test against
+        https://github.com/qiskit-community/qiskit-nature/issues/1291.
+        """
+        norb = 4
+        nelec = (4, 2)
+        ovlpab = np.asarray(
+            [
+                [0.987256, -0.001123, 0.00006, -0.0],
+                [-0.001123, -0.987256, -0.0, -0.00006],
+                [0.000019, 0.000055, -0.3195, -0.931662],
+                [0.000056, -0.000019, -0.931662, 0.3195],
+            ],
+        )
+
+        ang_mom = AngularMomentum(norb, ovlpab).second_q_ops()
+
+        mapper = ParityMapper(nelec)
+        qubit_op = mapper.map(ang_mom)
+
+        hf_state = HartreeFock(norb, nelec, mapper)
+
+        result = estimate_observables(Estimator(), hf_state, qubit_op)
+        self.assertAlmostEqual(result["AngularMomentum"][0], 1.9700743392855005)
+
 
 if __name__ == "__main__":
     unittest.main()