fix: the behavior of SpinOp.to_matrix for multiple spin sites (#1189)

* fix: the behavior of SpinOp.to_matrix for multiple spin sites

Closes #1187

* fix: the spin mappers after fixing the SpinOp.to_matrix method

This change is necessary because the `num_spins` object does not get set
to a lower bound during the `_validate_keys` method, since this gets
bypassed when setting `copy=False`.

qiskit_nature/second_q/operators/spin_op.py

@@ -338,7 +338,7 @@ def x(cls, spin: float | Fraction = Fraction(1, 2)) -> SpinOp:
         Returns:
             The X spin operator for ``spin``.
         """
-        return cls({"X_0": 1.0}, spin=spin, copy=False)
+        return cls({"X_0": 1.0}, spin=spin, num_spins=1, copy=False)
 
     @classmethod
     def y(cls, spin: float | Fraction = Fraction(1, 2)) -> SpinOp:
@@ -347,7 +347,7 @@ def y(cls, spin: float | Fraction = Fraction(1, 2)) -> SpinOp:
         Returns:
             The Y spin operator for ``spin``.
         """
-        return cls({"Y_0": 1.0}, spin=spin, copy=False)
+        return cls({"Y_0": 1.0}, spin=spin, num_spins=1, copy=False)
 
     @classmethod
     def z(cls, spin: float | Fraction = Fraction(1, 2)) -> SpinOp:
@@ -356,7 +356,7 @@ def z(cls, spin: float | Fraction = Fraction(1, 2)) -> SpinOp:
         Returns:
             The Z spin operator for ``spin``.
         """
-        return cls({"Z_0": 1.0}, spin=spin, copy=False)
+        return cls({"Z_0": 1.0}, spin=spin, num_spins=1, copy=False)
 
     @classmethod
     def one(cls, spin: float | Fraction = Fraction(1, 2)) -> SpinOp:
@@ -602,7 +602,11 @@ def to_matrix(self) -> np.ndarray:
         # reorder and expand
         simplified_op = self.index_order().simplify()
 
-        final_matrix = np.zeros((dim, dim), dtype=np.complex128)
+        # the size of the final matrix needs to adjust for the total number of spins which this
+        # operator acts on
+        final_dim = dim**self.register_length
+        final_matrix = np.zeros((final_dim, final_dim), dtype=np.complex128)
+
         for label, coeff in simplified_op.items():
             matrix_per_idx = {}
             # after .simplify() all exponents will be 1,
@@ -615,7 +619,9 @@ def to_matrix(self) -> np.ndarray:
                 matrix_per_idx[idx] = matrix_per_idx[idx] @ char_map.get(char, i_mat)
 
             # fill out empty indices with identity
-            dense_matrix_per_idx = [matrix_per_idx.get(i, i_mat) for i in range(len(self))]
+            dense_matrix_per_idx = [
+                matrix_per_idx.get(i, i_mat) for i in range(self.register_length)
+            ]
             # add weighted kronecker product to final matrix
             final_matrix += coeff * tensorall(np.asarray(dense_matrix_per_idx))
 

releasenotes/notes/fix-spin-op-to-matrix-e84594248fc49fd7.yaml

@@ -0,0 +1,5 @@
+---
+fixes:
+  - |
+    Fixes the behavior of :meth:`.SpinOp.to_matrix` for operators acting on more
+    than a single spin.

test/second_q/operators/test_spin_op.py

@@ -192,6 +192,48 @@ def test_to_matrix(self, label):
         actual = actual.to_matrix()
         np.testing.assert_array_almost_equal(actual, self.spin_1_matrix[label])
 
+    def test_to_matrix_multi_site(self):
+        """Test the to_matrix method when multiple sites are involved.
+
+        This is a regression test against https://github.com/Qiskit/qiskit-nature/issues/1187.
+        """
+        with self.subTest("spin 1/2; 2 sites"):
+            op = SpinOp({"X_0 X_1": 1, "Y_0 Y_1": 1, "Z_0 Z_1": 1})
+            actual = op.to_matrix()
+            expected = 0.25 * np.diag([1, -1, -1, 1])
+            expected[1, 2] = 0.5
+            expected[2, 1] = 0.5
+            np.testing.assert_array_almost_equal(actual, expected)
+
+        with self.subTest("spin 1; 2 sites"):
+            op = SpinOp({"X_0 X_1": 1, "Y_0 Y_1": 1, "Z_0 Z_1": 1}, spin=1)
+            actual = op.to_matrix()
+            expected = np.zeros((9, 9))
+            for i, j, v in zip(
+                [0, 1, 2, 2, 3, 4, 4, 5, 6, 6, 7, 8],
+                [0, 3, 2, 4, 1, 2, 6, 7, 4, 6, 5, 8],
+                [1, 1, -1, 1, 1, 1, 1, 1, 1, -1, 1, 1],
+            ):
+                expected[i, j] = v
+            np.testing.assert_array_almost_equal(actual, expected)
+
+        with self.subTest("spin 1/2; 3 sites"):
+            op = SpinOp({"X_0 X_1 X_2": 1, "Y_0 Y_2": 1, "Z_1": 1})
+            actual = op.to_matrix()
+            expected = np.zeros((8, 8))
+            for i, j, v in zip(
+                [0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7],
+                [0, 5, 7, 1, 4, 6, 2, 5, 7, 3, 4, 6, 1, 3, 4, 0, 2, 5, 1, 3, 6, 0, 2, 7],
+                # fmt: off
+                [
+                    0.5, -0.25, 0.125, 0.5, 0.25, 0.125, -0.5, 0.125, -0.25, -0.5, 0.125, 0.25,
+                     0.25, 0.125, 0.5, -0.25, 0.125, 0.5, 0.125, 0.25, -0.5, 0.125, -0.25, -0.5,
+                ],
+                # fmt: on
+            ):
+                expected[i, j] = v
+            np.testing.assert_array_almost_equal(actual, expected)
+
     def test_index_order(self):
         """Test index_order method"""
         with self.subTest("Test for single operators"):