Fix edge cases around settings.use_pauli_sum_op

Closes #1076

qiskit_nature/second_q/algorithms/excited_states_solvers/qeom.py

@@ -386,8 +386,13 @@ def solve(
             untap_aux_ops_sumop,  # Auxiliary observables
         ) = self.get_qubit_operators(problem, aux_operators)
 
-        untap_main_op = untap_main_op_sumop.primitive
-        untap_aux_ops = {key: op.primitive for key, op in untap_aux_ops_sumop.items()}
+        untap_main_op = untap_main_op_sumop
+        if isinstance(untap_main_op, PauliSumOp):
+            untap_main_op = untap_main_op.primitive
+        untap_aux_ops = {
+            key: op.primitive if isinstance(op, PauliSumOp) else op
+            for key, op in untap_aux_ops_sumop.items()
+        }
 
         # 2. Run ground state calculation with fully tapered custom auxiliary operators
         # Note that the solve() method includes the `second_q' auxiliary operators
@@ -710,7 +715,10 @@ def _prepare_expansion_basis(
         else:
             untap_hopping_ops = hopping_operators
 
-        untap_hopping_ops_sparse = {key: op.primitive for key, op in untap_hopping_ops.items()}
+        untap_hopping_ops_sparse = {
+            key: op.primitive if isinstance(op, PauliSumOp) else op
+            for key, op in untap_hopping_ops.items()
+        }
 
         return untap_hopping_ops_sparse, type_of_commutativities, size
 

qiskit_nature/second_q/circuit/library/ansatzes/succd.py

@@ -22,6 +22,7 @@
 import logging
 
 from qiskit.circuit import QuantumCircuit
+from qiskit.opflow import PauliSumOp
 from qiskit_nature import QiskitNatureError
 from qiskit_nature.deprecation import deprecate_arguments
 from qiskit_nature.second_q.mappers import QubitConverter, QubitMapper
@@ -148,6 +149,10 @@ def _filter_operators(self, operators):
         valid_operators, valid_excitations = [], []
         for op, ex in zip(operators, self._excitations_dict.values()):
             if op is not None:
+                # TODO: remove wrapping into PauliSumOp after the EvolvedOperatorAnsatz supports
+                # SparsePauliOp instances, too: https://github.com/Qiskit/qiskit-terra/pull/9537
+                if not isinstance(op, PauliSumOp):
+                    op = PauliSumOp(op)
                 valid_operators.append(op)
                 valid_excitations.extend(ex)
 

qiskit_nature/second_q/circuit/library/ansatzes/ucc.py

@@ -21,7 +21,7 @@
 
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import EvolvedOperatorAnsatz
-from qiskit.opflow import PauliTrotterEvolution
+from qiskit.opflow import PauliSumOp
 
 from qiskit_nature import QiskitNatureError
 from qiskit_nature.deprecation import deprecate_arguments, deprecate_property, warn_deprecated_type
@@ -216,7 +216,7 @@ def __init__(
         self._generalized = generalized
         self._preserve_spin = preserve_spin
 
-        super().__init__(reps=reps, evolution=PauliTrotterEvolution(), initial_state=initial_state)
+        super().__init__(reps=reps, initial_state=initial_state)
 
         # To give read access to the actual excitation list that UCC is using.
         self._excitation_list: list[tuple[tuple[int, ...], tuple[int, ...]]] | None = None
@@ -359,6 +359,10 @@ def _filter_operators(self, operators):
         valid_operators, valid_excitations = [], []
         for op, ex in zip(operators, self._excitation_list):
             if op is not None:
+                # TODO: remove wrapping into PauliSumOp after the EvolvedOperatorAnsatz supports
+                # SparsePauliOp instances, too: https://github.com/Qiskit/qiskit-terra/pull/9537
+                if not isinstance(op, PauliSumOp):
+                    op = PauliSumOp(op)
                 valid_operators.append(op)
                 valid_excitations.append(ex)
 
@@ -565,8 +569,8 @@ def _build_fermionic_excitation_ops(self, excitations: Sequence) -> list[Fermion
                 label.append(f"-_{unocc}")
             op = FermionicOp({" ".join(label): 1}, num_spin_orbitals=num_spin_orbitals)
             op -= op.adjoint()
-            # we need to account for an additional imaginary phase in the exponent (see also
-            # `PauliTrotterEvolution.convert`)
+            # we need to account for an additional imaginary phase in the exponent accumulated from
+            # the first-order trotterization routine implemented in Qiskit Terra
             op *= 1j  # type: ignore
             operators.append(op)
 

qiskit_nature/second_q/circuit/library/ansatzes/uvcc.py

@@ -20,7 +20,7 @@
 
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import EvolvedOperatorAnsatz
-from qiskit.opflow import PauliTrotterEvolution
+from qiskit.opflow import PauliSumOp
 
 from qiskit_nature import QiskitNatureError
 from qiskit_nature.deprecation import deprecate_arguments, deprecate_property, warn_deprecated_type
@@ -129,7 +129,7 @@ def __init__(
         self._num_modals = num_modals
         self._excitations = excitations
 
-        super().__init__(reps=reps, evolution=PauliTrotterEvolution(), initial_state=initial_state)
+        super().__init__(reps=reps, initial_state=initial_state)
 
         # To give read access to the actual excitation list that UVCC is using.
         self._excitation_list: list[tuple[tuple[int, ...], tuple[int, ...]]] | None = None
@@ -247,17 +247,24 @@ def operators(self):  # pylint: disable=invalid-overridden-method
                 else:
                     operators = self.qubit_mapper.map(excitation_ops)
 
-                valid_operators, valid_excitations = [], []
-                for op, ex in zip(operators, self._excitation_list):
-                    if op is not None:
-                        valid_operators.append(op)
-                        valid_excitations.append(ex)
-
-                self._excitation_list = valid_excitations
-                self.operators = valid_operators
+                self._filter_operators(operators=operators)
 
         return super(UVCC, self.__class__).operators.__get__(self)
 
+    def _filter_operators(self, operators):
+        valid_operators, valid_excitations = [], []
+        for op, ex in zip(operators, self._excitation_list):
+            if op is not None:
+                # TODO: remove wrapping into PauliSumOp after the EvolvedOperatorAnsatz supports
+                # SparsePauliOp instances, too: https://github.com/Qiskit/qiskit-terra/pull/9537
+                if not isinstance(op, PauliSumOp):
+                    op = PauliSumOp(op)
+                valid_operators.append(op)
+                valid_excitations.append(ex)
+
+        self._excitation_list = valid_excitations
+        self.operators = valid_operators
+
     def _invalidate(self):
         self._excitation_ops = None
         super()._invalidate()
@@ -388,8 +395,8 @@ def _build_vibration_excitation_ops(self, excitations: Sequence) -> list[Vibrati
                 label.append(f"-_{VibrationalOp.build_dual_index(self.num_modals, unocc)}")
             op = VibrationalOp({" ".join(label): 1}, self.num_modals)
             op -= op.adjoint()
-            # we need to account for an additional imaginary phase in the exponent (see also
-            # `PauliTrotterEvolution.convert`)
+            # we need to account for an additional imaginary phase in the exponent accumulated from
+            # the first-order trotterization routine implemented in Qiskit Terra
             op *= 1j  # type: ignore
             operators.append(op)
 

qiskit_nature/second_q/mappers/interleaved_qubit_mapper.py

@@ -96,7 +96,9 @@ def __init__(self, mapper: FermionicMapper):
     def _map_single(
         self, second_q_op: FermionicOp, *, register_length: int | None = None
     ) -> SparsePauliOp | PauliSumOp:
-        blocked_op = self.mapper._map_single(second_q_op, register_length=register_length).primitive
+        blocked_op = self.mapper._map_single(second_q_op, register_length=register_length)
+        if isinstance(blocked_op, PauliSumOp):
+            blocked_op = blocked_op.primitive
 
         def blocked_to_interleaved(label: str) -> str:
             return label[::2] + label[1::2]

qiskit_nature/second_q/mappers/qubit_converter.py

@@ -498,16 +498,21 @@ def symmetry_reduce_clifford(
         if converted_ops is None or self._z2symmetries is None or self._z2symmetries.is_empty():
             return_ops = converted_ops
         else:
-            wrapped_type = type(converted_ops)
-            wrapped_converted_ops: _ListOrDict[PauliSumOp] = _ListOrDict(converted_ops)
+            wrapped_converted_ops, wrapped_type = _ListOrDict.wrap(converted_ops)
+
+            pauli_sum_ops: ListOrDictType[PauliSumOp] = _ListOrDict()
+            for name, qubit_op in iter(wrapped_converted_ops):
+                if not isinstance(qubit_op, PauliSumOp):
+                    qubit_op = PauliSumOp(qubit_op)
+                pauli_sum_ops[name] = qubit_op
 
             if check_commutes:
                 logger.debug("Checking operators commute with symmetry:")
                 symmetry_ops = []
                 for sq_pauli in self._z2symmetries._sq_paulis:
                     symmetry_ops.append(PauliSumOp.from_list([(sq_pauli.to_label(), 1.0)]))
                 commuted = {}
-                for name, qubit_op in iter(wrapped_converted_ops):
+                for name, qubit_op in iter(pauli_sum_ops):
                     commutes = QubitConverter._check_commutes(symmetry_ops, qubit_op)
                     commuted[name] = commutes
                     logger.debug("Qubit operator '%s' commuted with symmetry: %s", name, commutes)
@@ -517,15 +522,13 @@ def symmetry_reduce_clifford(
                 for name, commutes in commuted.items():
                     if commutes:
                         tapered_qubit_ops[name] = self._z2symmetries.taper_clifford(
-                            wrapped_converted_ops[name]
+                            pauli_sum_ops[name]
                         )
             else:
                 logger.debug("Tapering operators whether they commute with symmetry or not:")
                 tapered_qubit_ops = _ListOrDict()
-                for name, qubit_op in iter(wrapped_converted_ops):
-                    tapered_qubit_ops[name] = self._z2symmetries.taper_clifford(
-                        wrapped_converted_ops[name]
-                    )
+                for name, qubit_op in iter(pauli_sum_ops):
+                    tapered_qubit_ops[name] = self._z2symmetries.taper_clifford(pauli_sum_ops[name])
 
             # NOTE: _ListOrDict.unwrap takes care of the conversion to/from PauliSumOp based on
             # settings.use_pauli_sum_op
@@ -549,19 +552,21 @@ def convert_clifford(
         if qubit_ops is None or self._z2symmetries is None or self._z2symmetries.is_empty():
             converted_ops = qubit_ops
         else:
-            if isinstance(qubit_ops, PauliSumOp):
-                converted_ops = self._z2symmetries.convert_clifford(qubit_ops)
-            else:
-                wrapped_type = type(qubit_ops)
-                wrapped_second_q_ops: _ListOrDict[PauliSumOp] = _ListOrDict(qubit_ops)
+            wrapped_second_q_ops, wrapped_type = _ListOrDict.wrap(qubit_ops)
 
-                converted_ops = _ListOrDict()
-                for name, second_q_op in iter(wrapped_second_q_ops):
-                    converted_ops[name] = self._z2symmetries.convert_clifford(second_q_op)
+            pauli_sum_ops: ListOrDictType[PauliSumOp] = _ListOrDict()
+            for name, qubit_op in iter(wrapped_second_q_ops):
+                if not isinstance(qubit_op, PauliSumOp):
+                    qubit_op = PauliSumOp(qubit_op)
+                pauli_sum_ops[name] = qubit_op
 
-                # NOTE: _ListOrDict.unwrap takes care of the conversion to/from PauliSumOp based on
-                # settings.use_pauli_sum_op
-                converted_ops = converted_ops.unwrap(wrapped_type)
+            converted_ops = _ListOrDict()
+            for name, second_q_op in iter(pauli_sum_ops):
+                converted_ops[name] = self._z2symmetries.convert_clifford(second_q_op)
+
+            # NOTE: _ListOrDict.unwrap takes care of the conversion to/from PauliSumOp based on
+            # settings.use_pauli_sum_op
+            converted_ops = converted_ops.unwrap(wrapped_type)
 
         return converted_ops
 

qiskit_nature/second_q/problems/base_problem.py

@@ -20,6 +20,7 @@
 import numpy as np
 from qiskit.algorithms.eigensolvers import EigensolverResult
 from qiskit.algorithms.minimum_eigensolvers import MinimumEigensolverResult
+from qiskit.opflow import PauliSumOp
 from qiskit.opflow.primitive_ops import Z2Symmetries as OpflowZ2Symmetries
 from qiskit.quantum_info.analysis.z2_symmetries import Z2Symmetries
 
@@ -149,7 +150,9 @@ def get_tapered_mapper(self, mapper: QubitMapper) -> TaperedQubitMapper:
             )
 
         qubit_op, _ = self.second_q_ops()
-        mapped_op = mapper.map(qubit_op).primitive
+        mapped_op = mapper.map(qubit_op)
+        if isinstance(mapped_op, PauliSumOp):
+            mapped_op = mapped_op.primitive
         z2_symmetries = Z2Symmetries.find_z2_symmetries(mapped_op)
         # pylint: disable=assignment-from-none
         # Known issue for abstract class methods https://github.com/PyCQA/pylint/issues/2559

test/second_q/algorithms/excited_state_solvers/test_qeom_electronic_ops.py

@@ -15,6 +15,7 @@
 from test import QiskitNatureTestCase
 
 from qiskit.utils import algorithm_globals
+from qiskit.opflow import PauliSumOp
 
 from qiskit_nature.units import DistanceUnit
 from qiskit_nature.second_q.mappers import QubitConverter, JordanWignerMapper, TaperedQubitMapper
@@ -67,7 +68,9 @@ def test_build_hopping_operators(self):
                 hopping_operators.keys(), expected_hopping_operators_electronic.keys()
             ):
                 self.assertEqual(key, exp_key)
-                val = hopping_operators[key].primitive
+                val = hopping_operators[key]
+                if isinstance(val, PauliSumOp):
+                    val = val.primitive
                 exp_val = expected_hopping_operators_electronic[exp_key]
                 if not val.equiv(exp_val):
                     print(val)
@@ -96,7 +99,9 @@ def test_build_hopping_operators_mapper(self):
                 hopping_operators.keys(), expected_hopping_operators_electronic.keys()
             ):
                 self.assertEqual(key, exp_key)
-                val = hopping_operators[key].primitive
+                val = hopping_operators[key]
+                if isinstance(val, PauliSumOp):
+                    val = val.primitive
                 exp_val = expected_hopping_operators_electronic[exp_key]
                 if not val.equiv(exp_val):
                     print(val)
@@ -125,7 +130,9 @@ def test_build_hopping_operators_taperedmapper(self):
                 hopping_operators.keys(), expected_hopping_operators_electronic.keys()
             ):
                 self.assertEqual(key, exp_key)
-                val = hopping_operators[key].primitive
+                val = hopping_operators[key]
+                if isinstance(val, PauliSumOp):
+                    val = val.primitive
                 exp_val = expected_hopping_operators_electronic[exp_key]
                 if not val.equiv(exp_val):
                     print(val)

test/second_q/algorithms/excited_state_solvers/test_qeom_vibrational_ops.py

@@ -17,6 +17,7 @@
 import unittest
 
 from qiskit.utils import algorithm_globals
+from qiskit.opflow import PauliSumOp
 
 from qiskit_nature.second_q.algorithms.excited_states_solvers.qeom_vibrational_ops_builder import (
     build_vibrational_ops,
@@ -98,7 +99,9 @@ def test_build_hopping_operators(self):
                 hopping_operators.keys(), expected_hopping_operators_vibrational.keys()
             ):
                 self.assertEqual(key, exp_key)
-                val = hopping_operators[key].primitive
+                val = hopping_operators[key]
+                if isinstance(val, PauliSumOp):
+                    val = val.primitive
                 exp_val = expected_hopping_operators_vibrational[exp_key]
                 if not val.equiv(exp_val):
                     print(val)
@@ -126,7 +129,9 @@ def test_build_hopping_operators_mapper(self):
                 hopping_operators.keys(), expected_hopping_operators_vibrational.keys()
             ):
                 self.assertEqual(key, exp_key)
-                val = hopping_operators[key].primitive
+                val = hopping_operators[key]
+                if isinstance(val, PauliSumOp):
+                    val = val.primitive
                 exp_val = expected_hopping_operators_vibrational[exp_key]
                 if not val.equiv(exp_val):
                     print(val)
@@ -154,7 +159,9 @@ def test_build_hopping_operators_taperedmapper(self):
                 hopping_operators.keys(), expected_hopping_operators_vibrational.keys()
             ):
                 self.assertEqual(key, exp_key)
-                val = hopping_operators[key].primitive
+                val = hopping_operators[key]
+                if isinstance(val, PauliSumOp):
+                    val = val.primitive
                 exp_val = expected_hopping_operators_vibrational[exp_key]
                 if not val.equiv(exp_val):
                     print(val)

test/second_q/hamiltonians/test_quadratic_hamiltonian.py

@@ -17,6 +17,7 @@
 import numpy as np
 from ddt import data, ddt
 from qiskit.quantum_info import random_hermitian
+from qiskit.opflow import PauliSumOp
 
 from qiskit_nature.second_q.hamiltonians import QuadraticHamiltonian
 from qiskit_nature.second_q.mappers import JordanWignerMapper, QubitConverter
@@ -105,11 +106,10 @@ def test_diagonalizing_bogoliubov_transform(self):
         )
 
         # confirm eigenvalues match with Jordan-Wigner transformed Hamiltonian
-        hamiltonian_jw = (
-            QubitConverter(mapper=JordanWignerMapper())
-            .convert(quad_ham.second_q_op())
-            .primitive.to_matrix()
-        )
+        qubit_op = QubitConverter(mapper=JordanWignerMapper()).convert(quad_ham.second_q_op())
+        if isinstance(qubit_op, PauliSumOp):
+            qubit_op = qubit_op.primitive
+        hamiltonian_jw = qubit_op.to_matrix()
         eigs, _ = np.linalg.eigh(hamiltonian_jw)
         expected_eigs = np.array(
             [

test/second_q/mappers/test_bksf_mapper.py

@@ -21,6 +21,7 @@
 
 import numpy as np
 from qiskit.quantum_info import SparsePauliOp, PauliList
+from qiskit.opflow import PauliSumOp
 
 from qiskit_nature.second_q.operators import FermionicOp
 from qiskit_nature.second_q.mappers import BravyiKitaevSuperFastMapper
@@ -162,10 +163,12 @@ def test_h2(self):
         aux = settings.use_pauli_sum_op
         try:
             settings.use_pauli_sum_op = True
-            pauli_sum_op = BravyiKitaevSuperFastMapper().map(h2_fop)
+            qubit_op = BravyiKitaevSuperFastMapper().map(h2_fop)
+            if isinstance(qubit_op, PauliSumOp):
+                qubit_op = qubit_op.primitive
 
             op1 = _sort_simplify(expected_pauli_op)
-            op2 = _sort_simplify(pauli_sum_op.primitive)
+            op2 = _sort_simplify(qubit_op)
 
             with self.subTest("Map H2 frome sto3g basis, number of terms"):
                 self.assertEqual(len(op1), len(op2))
@@ -245,8 +248,10 @@ def test_LiH(self):
         aux = settings.use_pauli_sum_op
         try:
             settings.use_pauli_sum_op = True
-            pauli_sum_op = BravyiKitaevSuperFastMapper().map(FERMIONIC_HAMILTONIAN)
-            self.assertEqual(pauli_sum_op._primitive, QUBIT_HAMILTONIAN)
+            qubit_op = BravyiKitaevSuperFastMapper().map(FERMIONIC_HAMILTONIAN)
+            if isinstance(qubit_op, PauliSumOp):
+                qubit_op = qubit_op.primitive
+            self.assertEqual(qubit_op, QUBIT_HAMILTONIAN)
             aux = settings.use_pauli_sum_op
 
             settings.use_pauli_sum_op = False