Improve wire ordering in tapering workflows

doc/releases/changelog-dev.md

@@ -160,6 +160,9 @@
 * Autograph can now be used with custom operations defined outside of the pennylane namespace.
   [(#6931)](https://github.com/PennyLaneAI/pennylane/pull/6931)
 
+* `qml.qchem.taper` now handles wire ordering for the tapered observables more robustly.
+  [(#)](https://github.com/PennyLaneAI/pennylane/pull/)
+
 <h3>Breaking changes 💔</h3>
 
 * `MultiControlledX` no longer accepts strings as control values.

pennylane/fermi/fermionic.py

@@ -721,7 +721,7 @@ def _to_string(fermi_op, of=False):
     pl_to_of_map = {"+": "^", "-": ""}
 
     if len(fermi_op) == 0:
-        return "I"
+        return "I" if not of else ""
 
     op_list = ["" for _ in range(len(fermi_op))]
     for loc, wire in fermi_op:

pennylane/qchem/tapering.py

@@ -293,44 +293,39 @@
     for ps in _split_pauli_sentence(ps_h, max_size=PAULI_SENTENCE_MEMORY_SPLITTING_SIZE):
         ts_ps += ps_u @ ps @ ps_u  # helps restrict the peak memory usage for u @ h @ u
 
-    wireset = ps_u.wires + ps_h.wires
+    wireset = ps_h.wires + ps_u.wires
     wiremap = dict(zip(list(wireset.toset()), range(len(wireset) + 1)))
     paulix_wires = [x.wires[0] for x in paulixops]
 
-    o = []
-    val = np.ones(len(ts_ps))
-
-    wires_tap = [i for i in ts_ps.wires if i not in paulix_wires]
-    wiremap_tap = dict(zip(wires_tap, range(len(wires_tap) + 1)))
-
-    for i, pw_coeff in enumerate(ts_ps.items()):
-        pw, _ = pw_coeff
+    wires_tap = [i for i in wiremap.keys() if i not in paulix_wires]
+    wires_ord = list(range(len(wires_tap)))
+    wiremap_tap = dict(zip(wires_tap, wires_ord))
 
+    obs, val = [], qml.math.ones(len(ts_ps))
+    for i, (pw, _) in enumerate(ts_ps.items()):
         for idx, w in enumerate(paulix_wires):
             if pw[w] == "X":
                 val[i] *= paulix_sector[idx]
 
-        o.append(
-            qml.pauli.string_to_pauli_word(
-                "".join([pw[wiremap[i]] for i in wires_tap]),
-                wire_map=wiremap_tap,
+        obs.append(
+            qml.pauli.PauliWord({wiremap_tap[wire]: pw[wire] for wire in wires_tap}).operation(
+                wire_order=wires_ord
             )
         )
 
-    c = qml.math.stack(qml.math.multiply(val * complex(1.0), list(ts_ps.values())))
+    interface = qml.math.get_deep_interface(list(ps_h.values()))
+    coeffs = qml.math.multiply(val, qml.math.array(list(ts_ps.values()), like=interface))
 
-    tapered_ham = qml.simplify(qml.dot(c, o))
-    # If simplified Hamiltonian is missing wires, then add wires manually for consistency
-    if set(wires_tap) != tapered_ham.wires.toset():
-        identity_op = functools.reduce(
-            lambda i, j: i @ j,
-            [
-                qml.Identity(wire)
-                for wire in Wires.unique_wires([tapered_ham.wires, Wires(wires_tap)])
-            ],
-        )
+    if interface == "jax" and qml.math.is_abstract(coeffs):
+        tapered_ham = qml.sum(*(qml.s_prod(coeff, op) for coeff, op in zip(coeffs, obs)))
+    else:
+        if qml.math.all(qml.math.abs(qml.math.imag(coeffs)) <= 1e-8):
+            coeffs = qml.math.real(coeffs)
+        tapered_ham = qml.simplify(0.0 * qml.Identity(wires=wires_ord) + qml.dot(coeffs, obs))
 
-        return tapered_ham + (0.0 * identity_op)
+    # If simplified Hamiltonian is missing wires, then add wires manually for consistency
+    if set(wires_ord) != tapered_ham.wires.toset():
+        return 0.0 * qml.Identity(wires=wires_ord) + tapered_ham
 
     return tapered_ham
 
@@ -481,7 +476,7 @@
     # taper the HF observable using the symmetries obtained from the molecular hamiltonian
     fermop_taper = _taper_pauli_sentence(ferm_ps, generators, paulixops, paulix_sector)
     fermop_ps = pauli_sentence(fermop_taper)
-    fermop_mat = _binary_matrix_from_pws(list(fermop_ps), len(fermop_taper.wires))
+    fermop_mat = _binary_matrix_from_pws(list(fermop_ps), len(fermop_ps.wires))
 
     # build a wireset to match wires with that of the tapered Hamiltonian
     gen_wires = Wires.all_wires([generator.wires for generator in generators])
@@ -612,7 +607,7 @@
 
 
 # pylint: disable=too-many-branches, too-many-arguments, inconsistent-return-statements, no-member
 def taper_operation(
     operation, generators, paulixops, paulix_sector, wire_order, op_wires=None, op_gen=None
 ):
     r"""Transform a gate operation with a Clifford operator and then taper qubits.