fix: make FermionicOp.simplify preserve the term order

qiskit_nature/second_q/operators/fermionic_op.py

@@ -518,50 +518,53 @@ def simplify(self, *, atol: float | None = None) -> FermionicOp:
     def _simplify_label(self, label: str, coeff: complex) -> tuple[str, complex]:
         bits = _BitsContainer()
 
+        # Since Python 3.7, dictionaries are guaranteed to be insert-order preserving. We use this
+        # to our advantage, to implement an ordered set, which allows us to preserve the label order
+        # and only remove canceling terms.
+        new_label: dict[str, None] = {}
+
         for lbl in label.split():
             char, index = lbl.split("_")
             idx = int(index)
             char_b = char == "+"
 
             if idx not in bits:
-                bits[idx] = int(f"{char_b:b}{not char_b:b}{char_b:b}{char_b:b}", base=2)
                 # we store all relevant information for each register index in 4 bits:
                 #   1. True if a `+` has been applied on this index
                 #   2. True if a `-` has been applied on this index
                 #   3. True if a `+` was applied first, False if a `-` was applied first
                 #   4. True if the last added operation on this index was `+`, False if `-`
+                bits[idx] = int(f"{char_b:b}{not char_b:b}{char_b:b}{char_b:b}", base=2)
+                # and we insert the encountered label into our ordered set
+                new_label[lbl] = None
 
             elif bits.get_last(idx) == char_b:
                 # we bail, if we apply the same operator as the last one
                 return "", 0
 
             elif bits.get_plus(idx) and bits.get_minus(idx):
-                # if both, `+` and `-`, have already been applied, we cancel the opposite to the
-                # current one (i.e. `+` will cancel `-` and vice versa)
+                # If both, `+` and `-`, have already been applied, we cancel the opposite to the
+                # current one (i.e. `+` will cancel `-` and vice versa). We update the bit container
                 bits.set_plus_or_minus(idx, not char_b, False)
+                # and pop the reverse label from the ordered set
+                new_label.pop(f"{'-' if char_b else '+'}_{idx}")
                 # we also update the last bit to the current char
                 bits.set_last(idx, char_b)
+                # and finally, we need to update the coefficient to account for the swap operations
+                # which were necessary to place the canceling terms next to each other
+                if idx != self.num_spin_orbitals:
+                    num_exchange = sum(i in bits for i in range(idx + 1, self.num_spin_orbitals))
+                    coeff *= -1 if num_exchange % 2 else 1
 
             else:
                 # else, we simply set the bit of the currently applied char
                 bits.set_plus_or_minus(idx, char_b, True)
+                # and track it in our ordered set
+                new_label[lbl] = None
                 # we also update the last bit to the current char
                 bits.set_last(idx, char_b)
 
-            if idx != self.num_spin_orbitals:
-                num_exchange = 0
-                for i in range(idx + 1, self.num_spin_orbitals):
-                    if i in bits:
-                        num_exchange += (bits.get_plus(i) + bits.get_minus(i)) % 2
-                coeff *= -1 if num_exchange % 2 else 1
-
-        new_label = []
-        for idx in sorted(bits):
-            plus = f"+_{idx}" if bits.get_plus(idx) else None
-            minus = f"-_{idx}" if bits.get_minus(idx) else None
-            new_label.extend([plus, minus] if bits.get_order(idx) else [minus, plus])
-
-        return " ".join(lbl for lbl in new_label if lbl is not None), coeff
+        return " ".join(new_label), coeff
 
 
 class _BitsContainer(MutableMapping):

test/second_q/operators/test_fermionic_op.py

@@ -84,7 +84,7 @@ def test_compose(self):
             ) @ FermionicOp({"": 1, "-_0 +_1": 1}, num_spin_orbitals=2)
             fer_op = fer_op.simplify()
             targ = FermionicOp(
-                {"+_0 +_1 -_1": 1, "-_0 +_0 -_1": 1, "+_0 -_0 +_1": 1, "-_0 -_1 +_1": -1},
+                {"+_0 +_1 -_1": 1, "-_0 +_0 -_1": 1, "+_0 +_1 -_0": 1, "-_0 -_1 +_1": -1},
                 num_spin_orbitals=2,
             )
             self.assertEqual(fer_op, targ)
@@ -155,12 +155,20 @@ def test_simplify(self):
             simplified_op = fer_op.simplify()
             self.assertEqual(simplified_op, fer_op)
 
+            fer_op = FermionicOp({"-_1 +_0": 1 + 0j}, num_spin_orbitals=2)
+            simplified_op = fer_op.simplify()
+            self.assertEqual(simplified_op, fer_op)
+
         with self.subTest("simplify zero"):
             fer_op = self.op1 - self.op1
             simplified_op = fer_op.simplify()
             targ = FermionicOp.zero()
             self.assertEqual(simplified_op, targ)
 
+        with self.subTest("simplify commutes with normal_ordered"):
+            fer_op = FermionicOp({"-_0 +_1": 1}, num_spin_orbitals=2)
+            self.assertEqual(fer_op.simplify().normal_ordered(), fer_op.normal_ordered().simplify())
+
     def test_hermiticity(self):
         """test is_hermitian"""
         with self.subTest("operator hermitian"):