feat: implement Fermionic.index_order

qiskit_nature/second_q/operators/fermionic_op.py

@@ -413,8 +413,11 @@ def transpose(self) -> FermionicOp:
 
         return self._new_instance(data)
 
-    def normal_ordered(self) -> FermionicOp:
-        """Convert to the equivalent operator with normal order.
+    def normal_order(self) -> FermionicOp:
+        """Convert to the equivalent operator in normal order.
+
+        The normal order for fermions is defined
+        [here](https://en.wikipedia.org/wiki/Normal_order#Fermions).
 
         Returns a new operator (the original operator is not modified).
 
@@ -434,7 +437,7 @@ def normal_ordered(self) -> FermionicOp:
         ordered_op = FermionicOp.zero()
 
         for terms, coeff in self.terms():
-            ordered_op += self._normal_ordered(terms, coeff)
+            ordered_op += self._normal_order(terms, coeff)
 
         # after successful normal ordering, we remove all zero coefficients
         return self._new_instance(
@@ -445,7 +448,7 @@ def normal_ordered(self) -> FermionicOp:
             }
         )
 
-    def _normal_ordered(self, terms: list[tuple[str, int]], coeff: complex) -> FermionicOp:
+    def _normal_order(self, terms: list[tuple[str, int]], coeff: complex) -> FermionicOp:
         if not terms:
             return self._new_instance({"": coeff})
 
@@ -468,7 +471,7 @@ def _normal_ordered(self, terms: list[tuple[str, int]], coeff: complex) -> Fermi
                         # a_i a_i^\dagger = 1 - a_i^\dagger a_i
                         new_terms = terms[: (j - 1)] + terms[(j + 1) :]
                         # we can do so by recursion on this method
-                        ordered_op += self._normal_ordered(new_terms, -1.0 * coeff)
+                        ordered_op += self._normal_order(new_terms, -1.0 * coeff)
 
                 elif right[0] == left[0]:
                     # when we have identical neighboring operators, differentiate two cases:
@@ -489,6 +492,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_order(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_order(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_order(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.
 
@@ -499,7 +551,7 @@ def is_hermitian(self, *, atol: float | None = None) -> bool:
             True if the operator is hermitian up to numerical tolerance, False otherwise.
         """
         atol = self.atol if atol is None else atol
-        diff = (self - self.adjoint()).normal_ordered().simplify(atol=atol)
+        diff = (self - self.adjoint()).normal_order().simplify(atol=atol)
         return all(np.isclose(coeff, 0.0, atol=atol) for coeff in diff.values())
 
     def simplify(self, *, atol: float | None = None) -> FermionicOp:

test/second_q/mappers/test_bksf_mapper.py

@@ -169,7 +169,7 @@ def test_h2(self):
             self.assertEqual(op1, op2)
 
         with self.subTest("Sparse FermionicOp input"):
-            h2_fop_sparse = h2_fop.normal_ordered()
+            h2_fop_sparse = h2_fop.normal_order()
             pauli_sum_op_from_sparse = BravyiKitaevSuperFastMapper().map(h2_fop_sparse)
             op2_from_sparse = _sort_simplify(pauli_sum_op_from_sparse.primitive)
             self.assertEqual(op1, op2_from_sparse)

test/second_q/operators/test_fermionic_op.py

@@ -165,9 +165,15 @@ def test_simplify(self):
             targ = FermionicOp.zero()
             self.assertEqual(simplified_op, targ)
 
-        with self.subTest("simplify commutes with normal_ordered"):
+        with self.subTest("simplify commutes with normal_order"):
             fer_op = FermionicOp({"-_0 +_1": 1}, num_spin_orbitals=2)
-            self.assertEqual(fer_op.simplify().normal_ordered(), fer_op.normal_ordered().simplify())
+            self.assertEqual(fer_op.simplify().normal_order(), fer_op.normal_order().simplify())
+
+        with self.subTest("simplify + index order"):
+            orig = FermionicOp({"+_1 -_0 +_0 -_0": 1, "-_0 +_1": 2})
+            fer_op = orig.simplify().index_order()
+            targ = FermionicOp({"-_0 +_1": 1})
+            self.assertEqual(fer_op, targ)
 
     def test_hermiticity(self):
         """test is_hermitian"""
@@ -282,47 +288,93 @@ def test_to_matrix(self):
                 binary = f"{idx:0{4}b}"
                 self.assertEqual(binary.count("1"), 2)
 
-    def test_normal_ordered(self):
-        """test normal_ordered method"""
+    def test_normal_order(self):
+        """test normal_order method"""
         with self.subTest("Test for creation operator"):
             orig = FermionicOp({"+_0": 1}, num_spin_orbitals=1)
-            fer_op = orig.normal_ordered()
+            fer_op = orig.normal_order()
             self.assertEqual(fer_op, orig)
 
         with self.subTest("Test for annihilation operator"):
             orig = FermionicOp({"-_0": 1}, num_spin_orbitals=1)
-            fer_op = orig.normal_ordered()
+            fer_op = orig.normal_order()
             self.assertEqual(fer_op, orig)
 
         with self.subTest("Test for number operator"):
             orig = FermionicOp({"+_0 -_0": 1}, num_spin_orbitals=1)
-            fer_op = orig.normal_ordered()
+            fer_op = orig.normal_order()
             self.assertEqual(fer_op, orig)
 
         with self.subTest("Test for empty operator"):
             orig = FermionicOp({"-_0 +_0": 1}, num_spin_orbitals=1)
-            fer_op = orig.normal_ordered()
+            fer_op = orig.normal_order()
             targ = FermionicOp({"": 1, "+_0 -_0": -1}, num_spin_orbitals=1)
             self.assertEqual(fer_op, targ)
 
         with self.subTest("Test for multiple operators 1"):
             orig = FermionicOp({"-_0 +_1": 1}, num_spin_orbitals=2)
-            fer_op = orig.normal_ordered()
+            fer_op = orig.normal_order()
             targ = FermionicOp({"+_1 -_0": -1}, num_spin_orbitals=2)
             self.assertEqual(fer_op, targ)
 
         with self.subTest("Test for multiple operators 2"):
             orig = FermionicOp({"-_0 +_0 +_1 -_2": 1}, num_spin_orbitals=3)
-            fer_op = orig.normal_ordered()
+            fer_op = orig.normal_order()
             targ = FermionicOp({"+_1 -_2": 1, "+_0 +_1 -_0 -_2": 1}, num_spin_orbitals=3)
             self.assertEqual(fer_op, targ)
 
         with self.subTest("Test normal ordering simplifies"):
             orig = FermionicOp({"-_0 +_1": 1, "+_1 -_0": -1, "+_0": 0.0}, num_spin_orbitals=2)
-            fer_op = orig.normal_ordered()
+            fer_op = orig.normal_order()
             targ = FermionicOp({"+_1 -_0": -2}, num_spin_orbitals=2)
             self.assertEqual(fer_op, targ)
 
+    def test_index_order(self):
+        """test index_order method"""
+        with self.subTest("Test for creation operator"):
+            orig = FermionicOp({"+_0": 1})
+            fer_op = orig.index_order()
+            self.assertEqual(fer_op, orig)
+
+        with self.subTest("Test for annihilation operator"):
+            orig = FermionicOp({"-_0": 1})
+            fer_op = orig.index_order()
+            self.assertEqual(fer_op, orig)
+
+        with self.subTest("Test for number operator"):
+            orig = FermionicOp({"+_0 -_0": 1})
+            fer_op = orig.index_order()
+            self.assertEqual(fer_op, orig)
+
+        with self.subTest("Test for empty operator"):
+            orig = FermionicOp({"-_0 +_0": 1})
+            fer_op = orig.index_order()
+            self.assertEqual(fer_op, orig)
+
+        with self.subTest("Test for multiple operators 1"):
+            orig = FermionicOp({"+_1 -_0": 1})
+            fer_op = orig.index_order()
+            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_order()
+            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_order()
+            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_order().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(