feat: implement Fermionic.index_ordered

Closes #901

qiskit_nature/second_q/operators/fermionic_op.py

@@ -489,6 +489,55 @@ def _normal_ordered(self, terms: list[tuple[str, int]], coeff: complex) -> Fermi
         ordered_op += self._new_instance({new_label: coeff})
         return ordered_op
 
+    def index_ordered(self) -> FermionicOp:
+        """Convert to the equivalent operator with the terms of each label ordered by index.
+
+        Returns a new operator (the original operator is not modified).
+
+        .. note::
+
+            You can use this method to achieve the most aggressive simplification of an operator
+            without changing the operation order per index. :meth:`simplify` does *not* reorder the
+            terms and, thus, cannot deduce ``-_0 +_1`` and ``+_1 -_0 +_0 -_0`` to be
+            identical labels. Calling this method will reorder the latter label to
+            ``-_0 +_0 -_0 +_1``, after which :meth:`simplify` will be able to correctly collapse
+            these two labels into one.
+
+        Returns:
+            The index ordered operator.
+        """
+        data = defaultdict(complex)  # type: dict[str, complex]
+        for terms, coeff in self.terms():
+            label, coeff = self._index_ordered(terms, coeff)
+            data[label] += coeff
+
+        # after successful index ordering, we remove all zero coefficients
+        return self._new_instance(
+            {
+                label: coeff
+                for label, coeff in data.items()
+                if not np.isclose(coeff, 0.0, atol=self.atol)
+            }
+        )
+
+    def _index_ordered(self, terms: list[tuple[str, int]], coeff: complex) -> tuple[str, complex]:
+        if not terms:
+            return "", coeff
+
+        # perform insertion sorting
+        for i in range(1, len(terms)):
+            for j in range(i, 0, -1):
+                right = terms[j]
+                left = terms[j - 1]
+
+                if left[1] > right[1]:
+                    terms[j - 1] = right
+                    terms[j] = left
+                    coeff *= -1.0
+
+        new_label = " ".join(f"{term[0]}_{term[1]}" for term in terms)
+        return new_label, coeff
+
     def is_hermitian(self, *, atol: float | None = None) -> bool:
         """Checks whether the operator is hermitian.
 

test/second_q/operators/test_fermionic_op.py

@@ -169,6 +169,12 @@ def test_simplify(self):
             fer_op = FermionicOp({"-_0 +_1": 1}, num_spin_orbitals=2)
             self.assertEqual(fer_op.simplify().normal_ordered(), fer_op.normal_ordered().simplify())
 
+        with self.subTest("simplify + index order"):
+            orig = FermionicOp({"+_1 -_0 +_0 -_0": 1, "-_0 +_1": 2})
+            fer_op = orig.simplify().index_ordered()
+            targ = FermionicOp({"-_0 +_1": 1})
+            self.assertEqual(fer_op, targ)
+
     def test_hermiticity(self):
         """test is_hermitian"""
         with self.subTest("operator hermitian"):
@@ -323,6 +329,52 @@ def test_normal_ordered(self):
             targ = FermionicOp({"+_1 -_0": -2}, num_spin_orbitals=2)
             self.assertEqual(fer_op, targ)
 
+    def test_index_ordered(self):
+        """test index_ordered method"""
+        with self.subTest("Test for creation operator"):
+            orig = FermionicOp({"+_0": 1})
+            fer_op = orig.index_ordered()
+            self.assertEqual(fer_op, orig)
+
+        with self.subTest("Test for annihilation operator"):
+            orig = FermionicOp({"-_0": 1})
+            fer_op = orig.index_ordered()
+            self.assertEqual(fer_op, orig)
+
+        with self.subTest("Test for number operator"):
+            orig = FermionicOp({"+_0 -_0": 1})
+            fer_op = orig.index_ordered()
+            self.assertEqual(fer_op, orig)
+
+        with self.subTest("Test for empty operator"):
+            orig = FermionicOp({"-_0 +_0": 1})
+            fer_op = orig.index_ordered()
+            self.assertEqual(fer_op, orig)
+
+        with self.subTest("Test for multiple operators 1"):
+            orig = FermionicOp({"+_1 -_0": 1})
+            fer_op = orig.index_ordered()
+            targ = FermionicOp({"-_0 +_1": -1})
+            self.assertEqual(fer_op, targ)
+
+        with self.subTest("Test for multiple operators 2"):
+            orig = FermionicOp({"+_2 -_0 +_1 -_0": 1, "-_0 +_1": 2})
+            fer_op = orig.index_ordered()
+            targ = FermionicOp({"-_0 -_0 +_1 +_2": 1, "-_0 +_1": 2})
+            self.assertEqual(fer_op, targ)
+
+        with self.subTest("Test index ordering simplifies"):
+            orig = FermionicOp({"-_0 +_1": 1, "+_1 -_0": -1, "+_0": 0.0})
+            fer_op = orig.index_ordered()
+            targ = FermionicOp({"-_0 +_1": 2})
+            self.assertEqual(fer_op, targ)
+
+        with self.subTest("index order + simplify"):
+            orig = FermionicOp({"+_1 -_0 +_0 -_0": 1, "-_0 +_1": 2})
+            fer_op = orig.index_ordered().simplify()
+            targ = FermionicOp({"-_0 +_1": 1})
+            self.assertEqual(fer_op, targ)
+
     def test_induced_norm(self):
         """Test induced norm."""
         op = 3 * FermionicOp({"+_0": 1}, num_spin_orbitals=1) + 4j * FermionicOp(