Fix TypeError on non-spin-symmetric hamiltonians in low_rank trotter decomposition (#1296)

Fixes #1105

`simulate_trotter(..., algorithm=LOW_RANK)` on spin-exchange
Hamiltonians used to fail with a vague `TypeError`. LOW_RANK with
`spin_basis=True` only supports interactions that are symmetric in the
α/β channel sense used for spatial downfolding; spin-exchange is not.

This PR validates that earlier and raises a clear `ValueError` (with a
hint to try `spin_basis=False`), updates docs/tests, and rebases on
current `main`.

---------

Co-authored-by: gemini-code-assist[bot] <176961590+gemini-code-assist[bot]@users.noreply.github.com>

Author: rosspeili

src/openfermion/circuits/low_rank.py

@@ -44,8 +44,13 @@ def get_chemist_two_body_coefficients(two_body_coefficients, spin_basis=True):
             giving the $g_{pqrs}$ tensor in chemist notation.
 
     Raises:
-        TypeError: Input must be two-body number conserving
-            FermionOperator or InteractionOperator.
+        ValueError: Input two-body tensor is not spin-symmetric. The LOW_RANK
+            decomposition requires a spin-symmetric interaction when
+            spin_basis=True. A spin-symmetric interaction satisfies
+            h[p,q,r,s] == h[p+1,q+1,r+1,s+1] for all even p, q, r, s
+            (i.e., the same coefficient for alpha and beta spin channels).
+            Consider passing spin_basis=False if your Hamiltonian is not
+            spin-symmetric.
     """
     # Initialize.
     n_orbitals = two_body_coefficients.shape[0]
@@ -57,10 +62,34 @@ def get_chemist_two_body_coefficients(two_body_coefficients, spin_basis=True):
         n_orbitals = n_orbitals // 2
         alpha_indices = list(range(0, n_orbitals * 2, 2))
         beta_indices = list(range(1, n_orbitals * 2, 2))
-        chemist_two_body_coefficients = chemist_two_body_coefficients[
+        # Extract the αα→ββ block, which is the only block used by the
+        # spatial-orbital low-rank decomposition.
+        alpha_alpha_beta_beta = chemist_two_body_coefficients[
             numpy.ix_(alpha_indices, alpha_indices, beta_indices, beta_indices)
         ]
 
+        # Validate spin-symmetry by checking that the extracted block is
+        # symmetric when reshaped to a matrix. For a spin-symmetric interaction
+        # the chemist tensor satisfies g[p,q,r,s] == g[r,s,p,q] (up to
+        # permutation symmetry), which makes the reshaped matrix symmetric.
+        # General spin-dependent interactions such as spin-exchange Hamiltonians
+        # produce an asymmetric matrix here and cannot be handled by this
+        # spatial-orbital downfolding approach.
+        flat = numpy.reshape(alpha_alpha_beta_beta, (n_orbitals**2, n_orbitals**2))
+        spin_asymmetry = numpy.amax(numpy.absolute(flat - flat.T))
+        if spin_asymmetry > EQ_TOLERANCE:
+            raise ValueError(
+                'The two-body tensor is not spin-symmetric. The LOW_RANK '
+                'decomposition requires a spin-symmetric interaction when '
+                'spin_basis=True (i.e., the same coefficients for alpha and '
+                'beta spin channels). Spin-dependent interactions such as '
+                'spin-exchange Hamiltonians violate this requirement. '
+                'Consider passing spin_basis=False if your Hamiltonian is '
+                'not spin-symmetric.'
+            )
+
+        chemist_two_body_coefficients = alpha_alpha_beta_beta
+
     # Determine a one body correction in the spin basis from spatial basis.
     one_body_correction = numpy.zeros((2 * n_orbitals, 2 * n_orbitals), complex)
     for p, q, r, s in itertools.product(range(n_orbitals), repeat=4):
@@ -106,7 +135,11 @@ def low_rank_two_body_decomposition(
             $\sum_{l=0}^{L-1} (\sum_{pq} |g_{lpq}|)^2 |\lambda_l| < x$
 
     Raises:
-        TypeError: Invalid two-body coefficient tensor specification.
+        ValueError: The two-body tensor failed symmetry or reality checks
+            required for the low-rank decomposition. When spin_basis=True,
+            the tensor must be spin-symmetric (see
+            get_chemist_two_body_coefficients()). When spin_basis=False,
+            the chemist-ordered tensor must be real and symmetric.
     """
     # Initialize N^2 by N^2 interaction array.
     one_body_correction, chemist_two_body_coefficients = get_chemist_two_body_coefficients(
@@ -117,10 +150,15 @@ def low_rank_two_body_decomposition(
     interaction_array = numpy.reshape(chemist_two_body_coefficients, (full_rank, full_rank))
 
     # Make sure interaction array is symmetric and real.
-    asymmetry = numpy.sum(numpy.absolute(interaction_array - interaction_array.transpose()))
-    imaginary_norm = numpy.sum(numpy.absolute(interaction_array.imag))
+    asymmetry = numpy.amax(numpy.absolute(interaction_array - interaction_array.transpose()))
+    imaginary_norm = numpy.amax(numpy.absolute(interaction_array.imag))
     if asymmetry > EQ_TOLERANCE or imaginary_norm > EQ_TOLERANCE:
-        raise TypeError('Invalid two-body coefficient tensor specification.')
+        raise ValueError(
+            'The two-body coefficient tensor failed the symmetry or reality '
+            'checks required by the low-rank decomposition. If spin_basis=True, '
+            'ensure the Hamiltonian is spin-symmetric. If spin_basis=False, '
+            'ensure the chemist-ordered two-body tensor is real and symmetric.'
+        )
 
     # Decompose with exact diagonalization.
     eigenvalues, eigenvectors = numpy.linalg.eigh(interaction_array)

src/openfermion/circuits/low_rank_test.py

@@ -163,8 +163,8 @@ def test_molecular_operator_consistency(self):
         )
         self.assertAlmostEqual(trunc_error, 0.0)
 
-        # Check for property errors
-        with self.assertRaises(TypeError):
+        # Check for property errors — imaginary tensor should raise ValueError
+        with self.assertRaises(ValueError):
             eigenvalues, one_body_squares, _, trunc_error = low_rank_two_body_decomposition(
                 two_body_coefficients + 0.01j, truncation_threshold=1.0, final_rank=1
             )
@@ -315,3 +315,39 @@ def test_one_body_squared_nonhermitian_raises_error(self):
         one_body_matrix = numpy.array([[0, 1], [0, 0]])
         with self.assertRaises(ValueError):
             prepare_one_body_squared_evolution(one_body_matrix, spin_basis=False)
+
+
+class SpinExchangeTest(unittest.TestCase):
+    """Tests for spin-symmetry validation in the low-rank decomposition."""
+
+    def _make_spin_exchange_tensor(self):
+        """Return the two_body_tensor for H = a^dag_0 a^dag_3 a_1 a_2 + h.c."""
+        from openfermion import get_interaction_operator
+
+        h_sf = FermionOperator('0^ 3^ 1 2', 1.0) + FermionOperator('2^ 1^ 3 0', 1.0)
+        return get_interaction_operator(h_sf, n_qubits=4).two_body_tensor
+
+    def test_get_chemist_two_body_coefficients_raises_for_spin_exchange(self):
+        """Non-spin-symmetric input raises ValueError with informative message."""
+        two_body_tensor = self._make_spin_exchange_tensor()
+        with self.assertRaises(ValueError) as ctx:
+            get_chemist_two_body_coefficients(two_body_tensor, spin_basis=True)
+        self.assertIn('spin-symmetric', str(ctx.exception))
+
+    def test_low_rank_decomposition_raises_for_spin_exchange(self):
+        """Non-spin-symmetric input raises ValueError with informative message."""
+        two_body_tensor = self._make_spin_exchange_tensor()
+        with self.assertRaises(ValueError) as ctx:
+            low_rank_two_body_decomposition(two_body_tensor, spin_basis=True)
+        self.assertIn('spin-symmetric', str(ctx.exception))
+
+    def test_spin_symmetric_hamiltonian_succeeds(self):
+        """Spin-symmetric Hamiltonians decompose without error."""
+        filename = os.path.join(DATA_DIRECTORY, 'H2_sto-3g_singlet_0.7414')
+        molecule = MolecularData(filename=filename)
+        two_body_coefficients = molecule.get_molecular_hamiltonian().two_body_tensor
+        _, _ = get_chemist_two_body_coefficients(two_body_coefficients, spin_basis=True)
+        eigenvalues, _, _, _ = low_rank_two_body_decomposition(
+            two_body_coefficients, spin_basis=True
+        )
+        self.assertGreater(len(eigenvalues), 0)

src/openfermion/circuits/trotter/algorithms/low_rank.py

@@ -64,6 +64,13 @@ class LowRankTrotterAlgorithm(TrotterAlgorithm):
     or it is chosen so that
     $\sum_{l=0}^{L-1} (\sum_{pq} |g_{lpq}|)^2 |\lambda_l| < x$
     where x is a truncation threshold specified by user.
+
+    Note:
+        When spin_basis=True (the default), the input
+        InteractionOperator must have a
+        spin-symmetric two-body tensor, i.e. identical interaction
+        coefficients for the alpha and beta spin channels. Hamiltonians
+        that break this symmetry will raise a ValueError.
     """
 
     supported_types = {ops.InteractionOperator}

src/openfermion/circuits/trotter/simulate_trotter_test.py

@@ -499,3 +499,19 @@ def test_trotter_misspecified_control_raises_error(algorithm_type, hamiltonian):
             next(algorithm.trotter_step(qubits, time))
         with pytest.raises(TypeError):
             next(algorithm.trotter_step(qubits, time, control_qubit=2))
+
+
+def test_simulate_trotter_spin_exchange_raises_value_error():
+    """LOW_RANK raises ValueError for a non-spin-symmetric Hamiltonian."""
+    h_sf = openfermion.FermionOperator('0^ 3^ 1 2', 1.0) + openfermion.FermionOperator(
+        '2^ 1^ 3 0', 1.0
+    )
+    interaction_hamiltonian = openfermion.get_interaction_operator(h_sf, n_qubits=4)
+    qubits = cirq.LineQubit.range(4)
+
+    with pytest.raises(ValueError, match='spin-symmetric'):
+        next(
+            simulate_trotter(
+                qubits=qubits, hamiltonian=interaction_hamiltonian, time=1.0, algorithm=LOW_RANK
+            )
+        )