diff --git a/python/ffsim/protocols/trace_protocol.py b/python/ffsim/protocols/trace_protocol.py
index f53be287e..4b1a62ccd 100644
--- a/python/ffsim/protocols/trace_protocol.py
+++ b/python/ffsim/protocols/trace_protocol.py
@@ -53,6 +53,16 @@ def trace(obj: Any, norb: int, nelec: tuple[int, int]) -> complex:
     )
 
 
+def _term_phase(op: tuple[tuple[bool, bool, int], ...]) -> int:
+    phase = 0
+    while len(op) > 0:
+        action, spin, orb = op[0]
+        conj_idx = op.index((not action, spin, orb))
+        phase += conj_idx - 1
+        op = op[1:conj_idx] + op[conj_idx + 1 :]
+    return (-1) ** phase
+
+
 def _trace_term(
     op: tuple[tuple[bool, bool, int], ...], norb: int, nelec: tuple[int, int]
 ) -> complex:
@@ -110,9 +120,12 @@ def _trace_term(
     # the trace is nontrival and is a product of
     # binom(#orbs not in the support of op, #elec on these orbs)
     # for each spin species
+    # and a phase that depends on the ordering between the actions
 
-    return math.comb(norb - norb_alpha, n_alpha - nelec_alpha) * math.comb(
-        norb - norb_beta, n_beta - nelec_beta
+    return (
+        _term_phase(op)
+        * math.comb(norb - norb_alpha, n_alpha - nelec_alpha)
+        * math.comb(norb - norb_beta, n_beta - nelec_beta)
     )
 
 
diff --git a/tests/python/operators/fermion_operator_test.py b/tests/python/operators/fermion_operator_test.py
index 1bec6fdf3..8f4550516 100644
--- a/tests/python/operators/fermion_operator_test.py
+++ b/tests/python/operators/fermion_operator_test.py
@@ -645,10 +645,20 @@ def test_mapping_methods():
     }
 
 def test_trace():
+    rng = np.random.default_rng(12345)
+    # compare for random_diagonal_coulomb_hamiltonian
     norb = 50
     nelec = (3, 3)
-    rng = np.random.default_rng(12345)
     ham = ffsim.random.random_diagonal_coulomb_hamiltonian(norb, real=True, seed=rng)
     t1 = ffsim.trace(ffsim.fermion_operator(ham), norb=norb, nelec=nelec)
     t2 = ffsim.trace(ham, norb=norb, nelec=nelec)
     np.testing.assert_allclose(t1, t2)
+
+    # compare for random_molecular_hamiltonian
+    norb = 20
+    nelec = (3, 3)
+    ham = ffsim.random.random_molecular_hamiltonian(norb, seed=rng, dtype=float)
+    t1 = ffsim.trace(ffsim.fermion_operator(ham), norb=norb, nelec=nelec)
+    t2 = ffsim.trace(ham, norb=norb, nelec=nelec)
+
+