From 1e14e1d20600a953c24cc3e5fea8a8285a3e3e66 Mon Sep 17 00:00:00 2001
From: Max Rossmannek <oss@zurich.ibm.com>
Date: Sat, 6 Aug 2022 15:00:37 +0200
Subject: [PATCH] refactor: match return type of `BaseProblem.second_q_ops()`
 to Terra

The algorithms in Terra always take a tuple of the main operator and
some auxiliary operators. This commit changes the `BaseProblem`
interface to already return the generated operators in the same format.
This simplifies the handling of the main operator, a scenario that
occurs very often and required manual intervention thus far.
---
 README.md                                     |  3 +-
 docs/tutorials/01_electronic_structure.ipynb  |  8 ++--
 docs/tutorials/02_vibrational_structure.ipynb | 10 ++---
 docs/tutorials/08_property_framework.ipynb    | 38 +++++++++----------
 .../excited_states_eigensolver.py             | 10 +----
 .../algorithms/excited_states_solvers/qeom.py |  6 +--
 .../ground_state_eigensolver.py               | 10 +----
 .../second_q/problems/base_problem.py         |  9 ++---
 .../problems/electronic_structure_problem.py  | 14 +++----
 .../problems/lattice_model_problem.py         |  9 ++---
 .../problems/vibrational_structure_problem.py | 15 +++-----
 .../test_advanced_ucc_variants.py             |  5 +--
 .../test_groundstate_eigensolver.py           | 12 ++----
 test/second_q/mappers/test_qubit_converter.py |  3 +-
 .../test_electronic_structure_problem.py      | 22 +++++------
 .../problems/test_lattice_model_problem.py    |  4 +-
 .../test_vibrational_structure_problem.py     | 10 ++---
 test/test_end2end_with_vqe.py                 |  6 +--
 18 files changed, 78 insertions(+), 116 deletions(-)

diff --git a/README.md b/README.md
index 2863482712..4ea39e85e9 100644
--- a/README.md
+++ b/README.md
@@ -81,8 +81,7 @@ driver = PySCFDriver(atom='H .0 .0 .0; H .0 .0 0.735',
 problem = ElectronicStructureProblem(driver)
 
 # generate the second-quantized operators
-second_q_ops = problem.second_q_ops()
-main_op = second_q_ops['ElectronicEnergy']
+main_op, _ = problem.second_q_ops()
 
 particle_number = problem.grouped_property_transformed.get_property("ParticleNumber")
 
diff --git a/docs/tutorials/01_electronic_structure.ipynb b/docs/tutorials/01_electronic_structure.ipynb
index 5d86d781e6..def6b6d56d 100644
--- a/docs/tutorials/01_electronic_structure.ipynb
+++ b/docs/tutorials/01_electronic_structure.ipynb
@@ -189,8 +189,8 @@
    ],
    "source": [
     "es_problem = ElectronicStructureProblem(driver)\n",
-    "second_q_op = es_problem.second_q_ops()\n",
-    "print(second_q_op[0])"
+    "main_op, aux_ops = es_problem.second_q_ops()\n",
+    "print(main_op)"
    ]
   },
   {
@@ -229,7 +229,7 @@
    ],
    "source": [
     "qubit_converter = QubitConverter(mapper=JordanWignerMapper())\n",
-    "qubit_op = qubit_converter.convert(second_q_op[0])\n",
+    "qubit_op = qubit_converter.convert(main_op)\n",
     "print(qubit_op)"
    ]
   },
@@ -259,7 +259,7 @@
    ],
    "source": [
     "qubit_converter = QubitConverter(mapper=ParityMapper(), two_qubit_reduction=True)\n",
-    "qubit_op = qubit_converter.convert(second_q_op[0], num_particles=es_problem.num_particles)\n",
+    "qubit_op = qubit_converter.convert(main_op, num_particles=es_problem.num_particles)\n",
     "print(qubit_op)"
    ]
   },
diff --git a/docs/tutorials/02_vibrational_structure.ipynb b/docs/tutorials/02_vibrational_structure.ipynb
index 45a2c7d45e..3f2d53fcbf 100644
--- a/docs/tutorials/02_vibrational_structure.ipynb
+++ b/docs/tutorials/02_vibrational_structure.ipynb
@@ -201,7 +201,7 @@
     "from qiskit_nature.second_q.mappers import DirectMapper\n",
     "\n",
     "vibrational_problem = VibrationalStructureProblem(driver, num_modals=2, truncation_order=2)\n",
-    "second_q_ops = vibrational_problem.second_q_ops()"
+    "main_op, aux_ops = vibrational_problem.second_q_ops()"
    ]
   },
   {
@@ -270,7 +270,7 @@
     }
    ],
    "source": [
-    "print(second_q_ops[0])"
+    "print(main_op)"
    ]
   },
   {
@@ -342,7 +342,7 @@
    ],
    "source": [
     "qubit_converter = QubitConverter(mapper=DirectMapper())\n",
-    "qubit_op = qubit_converter.convert(second_q_ops[0])\n",
+    "qubit_op = qubit_converter.convert(main_op)\n",
     "\n",
     "print(qubit_op)"
    ]
@@ -611,11 +611,11 @@
    ],
    "source": [
     "vibrational_problem = VibrationalStructureProblem(driver, num_modals=3, truncation_order=2)\n",
-    "second_q_ops = vibrational_problem.second_q_ops()\n",
+    "main_op, aux_ops = vibrational_problem.second_q_ops()\n",
     "\n",
     "qubit_converter = QubitConverter(mapper=DirectMapper())\n",
     "\n",
-    "qubit_op = qubit_converter.convert(second_q_ops[0])\n",
+    "qubit_op = qubit_converter.convert(main_op)\n",
     "\n",
     "print(qubit_op)"
    ]
diff --git a/docs/tutorials/08_property_framework.ipynb b/docs/tutorials/08_property_framework.ipynb
index a52a1aa986..36f048f36b 100644
--- a/docs/tutorials/08_property_framework.ipynb
+++ b/docs/tutorials/08_property_framework.ipynb
@@ -454,7 +454,7 @@
     }
    ],
    "source": [
-    "hamiltonian = electronic_energy.second_q_ops()[0]  # here, output length is always 1\n",
+    "hamiltonian = electronic_energy.second_q_ops()[\"ElectronicEnergy\"]\n",
     "print(hamiltonian)"
    ]
   },
@@ -957,7 +957,7 @@
    "outputs": [],
    "source": [
     "from itertools import product\n",
-    "from typing import List\n",
+    "from typing import Dict\n",
     "\n",
     "import h5py\n",
     "\n",
@@ -1001,8 +1001,8 @@
     "    def from_hdf5(cls, h5py_group: h5py.Group) -> \"ElectronicDensity\":\n",
     "        return ElectronicDensity(h5py_group.attrs[\"num_molecular_orbitals\"])\n",
     "\n",
-    "    def second_q_ops(self) -> List[FermionicOp]:\n",
-    "        ops = []\n",
+    "    def second_q_ops(self) -> Dict[str, FermionicOp]:\n",
+    "        ops = {}\n",
     "\n",
     "        # iterate all pairs of molecular orbitals\n",
     "        for mo_i, mo_j in product(range(self._num_molecular_orbitals), repeat=2):\n",
@@ -1017,7 +1017,7 @@
     "            one_body_ints = OneBodyElectronicIntegrals(\n",
     "                ElectronicBasis.MO, (number_op_matrix, number_op_matrix)\n",
     "            )\n",
-    "            ops.append(one_body_ints.to_second_q_op())\n",
+    "            ops[f\"{mo_i}_{mo_j}\"] = one_body_ints.to_second_q_op()\n",
     "\n",
     "        return ops\n",
     "\n",
@@ -1065,28 +1065,28 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "0 : Fermionic Operator\n",
+      "0_0 : Fermionic Operator\n",
       "register length=4, number terms=2\n",
-      "  (1+0j) * ( +_0 -_0 )\n",
-      "+ (1+0j) * ( +_2 -_2 )\n",
-      "1 : Fermionic Operator\n",
+      "  1.0 * ( +_0 -_0 )\n",
+      "+ 1.0 * ( +_2 -_2 )\n",
+      "0_1 : Fermionic Operator\n",
       "register length=4, number terms=2\n",
-      "  (1+0j) * ( +_0 -_1 )\n",
-      "+ (1+0j) * ( +_2 -_3 )\n",
-      "2 : Fermionic Operator\n",
+      "  1.0 * ( +_0 -_1 )\n",
+      "+ 1.0 * ( +_2 -_3 )\n",
+      "1_0 : Fermionic Operator\n",
       "register length=4, number terms=2\n",
-      "  (1+0j) * ( +_1 -_0 )\n",
-      "+ (1+0j) * ( +_3 -_2 )\n",
-      "3 : Fermionic Operator\n",
+      "  1.0 * ( +_1 -_0 )\n",
+      "+ 1.0 * ( +_3 -_2 )\n",
+      "1_1 : Fermionic Operator\n",
       "register length=4, number terms=2\n",
-      "  (1+0j) * ( +_1 -_1 )\n",
-      "+ (1+0j) * ( +_3 -_3 )\n"
+      "  1.0 * ( +_1 -_1 )\n",
+      "+ 1.0 * ( +_3 -_3 )\n"
      ]
     }
    ],
    "source": [
-    "for idx, op in enumerate(density.second_q_ops()):\n",
-    "    print(idx, \":\", op)"
+    "for key, op in density.second_q_ops().items():\n",
+    "    print(key, \":\", op)"
    ]
   },
   {
diff --git a/qiskit_nature/second_q/algorithms/excited_states_solvers/excited_states_eigensolver.py b/qiskit_nature/second_q/algorithms/excited_states_solvers/excited_states_eigensolver.py
index 2f0c0c3123..b7c058d453 100644
--- a/qiskit_nature/second_q/algorithms/excited_states_solvers/excited_states_eigensolver.py
+++ b/qiskit_nature/second_q/algorithms/excited_states_solvers/excited_states_eigensolver.py
@@ -66,15 +66,7 @@ def get_qubit_operators(
         """Gets the operator and auxiliary operators, and transforms the provided auxiliary operators"""
         # Note that ``aux_ops`` contains not only the transformed ``aux_operators`` passed by the
         # user but also additional ones from the transformation
-        second_q_ops = problem.second_q_ops()
-        name = problem.main_property_name
-        main_second_q_op = second_q_ops.pop(name, None)
-        if main_second_q_op is None:
-            raise ValueError(
-                f"The main `SecondQuantizedOp` associated with the {name} property cannot be "
-                "`None`."
-            )
-        aux_second_q_ops = second_q_ops
+        main_second_q_op, aux_second_q_ops = problem.second_q_ops()
 
         main_operator = self._qubit_converter.convert(
             main_second_q_op,
diff --git a/qiskit_nature/second_q/algorithms/excited_states_solvers/qeom.py b/qiskit_nature/second_q/algorithms/excited_states_solvers/qeom.py
index 4ee78635e6..fd9c7a7b5f 100644
--- a/qiskit_nature/second_q/algorithms/excited_states_solvers/qeom.py
+++ b/qiskit_nature/second_q/algorithms/excited_states_solvers/qeom.py
@@ -118,11 +118,7 @@ def solve(
         groundstate_result = self._gsc.solve(problem, aux_operators)
 
         # 2. Prepare the excitation operators
-        second_q_ops = problem.second_q_ops()
-        if isinstance(second_q_ops, list):
-            main_second_q_op = second_q_ops[0]
-        elif isinstance(second_q_ops, dict):
-            main_second_q_op = second_q_ops.pop(problem.main_property_name)
+        main_second_q_op, _ = problem.second_q_ops()
 
         self._untapered_qubit_op_main = self._gsc.qubit_converter.convert_only(
             main_second_q_op, problem.num_particles
diff --git a/qiskit_nature/second_q/algorithms/ground_state_solvers/ground_state_eigensolver.py b/qiskit_nature/second_q/algorithms/ground_state_solvers/ground_state_eigensolver.py
index 6ef4a413ac..27f0e9d3cb 100644
--- a/qiskit_nature/second_q/algorithms/ground_state_solvers/ground_state_eigensolver.py
+++ b/qiskit_nature/second_q/algorithms/ground_state_solvers/ground_state_eigensolver.py
@@ -102,15 +102,7 @@ def get_qubit_operators(
         """Gets the operator and auxiliary operators, and transforms the provided auxiliary operators"""
         # Note that ``aux_ops`` contains not only the transformed ``aux_operators`` passed by the
         # user but also additional ones from the transformation
-        second_q_ops = problem.second_q_ops()
-        name = problem.main_property_name
-        main_second_q_op = second_q_ops.pop(name, None)
-        if main_second_q_op is None:
-            raise ValueError(
-                f"The main `SecondQuantizedOp` associated with the {name} property cannot be "
-                "`None`."
-            )
-        aux_second_q_ops = second_q_ops
+        main_second_q_op, aux_second_q_ops = problem.second_q_ops()
         main_operator = self._qubit_converter.convert(
             main_second_q_op,
             num_particles=problem.num_particles,
diff --git a/qiskit_nature/second_q/problems/base_problem.py b/qiskit_nature/second_q/problems/base_problem.py
index a22897dd7e..2ed37aa45b 100644
--- a/qiskit_nature/second_q/problems/base_problem.py
+++ b/qiskit_nature/second_q/problems/base_problem.py
@@ -80,13 +80,12 @@ def num_particles(self) -> Optional[Tuple[int, int]]:
         return None
 
     @abstractmethod
-    def second_q_ops(self) -> dict[str, SecondQuantizedOp]:
-        """Returns the second quantized operators associated with this Property.
-
-        The actual return-type is determined by `qiskit_nature.settings.dict_aux_operators`.
+    def second_q_ops(self) -> tuple[SecondQuantizedOp, dict[str, SecondQuantizedOp]]:
+        """Returns the second quantized operators associated with this problem.
 
         Returns:
-            A `list` or `dict` of `SecondQuantizedOp` objects.
+            A tuple, with the first object being the main operator and the second being a dictionary
+            of auxiliary operators.
         """
         raise NotImplementedError()
 
diff --git a/qiskit_nature/second_q/problems/electronic_structure_problem.py b/qiskit_nature/second_q/problems/electronic_structure_problem.py
index cf12b7c469..b612fca1d8 100644
--- a/qiskit_nature/second_q/problems/electronic_structure_problem.py
+++ b/qiskit_nature/second_q/problems/electronic_structure_problem.py
@@ -111,17 +111,12 @@ def num_spin_orbitals(self) -> int:
             )
         return self._grouped_property_transformed.get_property("ParticleNumber").num_spin_orbitals
 
-    def second_q_ops(self) -> dict[str, SecondQuantizedOp]:
+    def second_q_ops(self) -> tuple[SecondQuantizedOp, dict[str, SecondQuantizedOp]]:
         """Returns the second quantized operators associated with this Property.
 
-        If the arguments are returned as a `list`, the operators are in the following order: the
-        Hamiltonian operator, total particle number operator, total angular momentum operator, total
-        magnetization operator, and (if available) x, y, z dipole operators.
-
-        The actual return-type is determined by `qiskit_nature.settings.dict_aux_operators`.
-
         Returns:
-            A `list` or `dict` of `SecondQuantizedOp` objects.
+            A tuple, with the first object being the main operator and the second being a dictionary
+            of auxiliary operators.
         """
         driver_result = self.driver.run()
 
@@ -129,8 +124,9 @@ def second_q_ops(self) -> dict[str, SecondQuantizedOp]:
         self._grouped_property_transformed = self._transform(self._grouped_property)
 
         second_quantized_ops = self._grouped_property_transformed.second_q_ops()
+        main_op = second_quantized_ops.pop(self._main_property_name)
 
-        return second_quantized_ops
+        return main_op, second_quantized_ops
 
     def hopping_qeom_ops(
         self,
diff --git a/qiskit_nature/second_q/problems/lattice_model_problem.py b/qiskit_nature/second_q/problems/lattice_model_problem.py
index f733f38692..84eb96ccf8 100644
--- a/qiskit_nature/second_q/problems/lattice_model_problem.py
+++ b/qiskit_nature/second_q/problems/lattice_model_problem.py
@@ -41,16 +41,15 @@ def __init__(self, lattice_model: LatticeModel) -> None:
         super().__init__(main_property_name="LatticeEnergy")
         self._lattice_model = lattice_model
 
-    def second_q_ops(self) -> dict[str, SecondQuantizedOp]:
+    def second_q_ops(self) -> tuple[SecondQuantizedOp, dict[str, SecondQuantizedOp]]:
         """Returns the second quantized operators created based on the lattice models.
 
         Returns:
-            A ``list`` or ``dict`` of
-            :class:`~qiskit_nature.second_q.operators.SecondQuantizedOp`
+            A tuple, with the first object being the main operator and the second being a dictionary
+            of auxiliary operators.
         """
         second_q_op = self._lattice_model.second_q_ops()
-        second_q_ops = {self._main_property_name: second_q_op}
-        return second_q_ops
+        return second_q_op, {}
 
     def interpret(
         self,
diff --git a/qiskit_nature/second_q/problems/vibrational_structure_problem.py b/qiskit_nature/second_q/problems/vibrational_structure_problem.py
index 7e3ac78604..db49968e27 100644
--- a/qiskit_nature/second_q/problems/vibrational_structure_problem.py
+++ b/qiskit_nature/second_q/problems/vibrational_structure_problem.py
@@ -59,16 +59,12 @@ def __init__(
         self.num_modals = num_modals
         self.truncation_order = truncation_order
 
-    def second_q_ops(self) -> dict[str, SecondQuantizedOp]:
-        """Returns the second quantized operators created based on the driver and transformations.
-
-        If the arguments are returned as a `list`, the operators are in the following order: the
-        Vibrational Hamiltonian operator, occupied modal operators for each mode.
-
-        The actual return-type is determined by `qiskit_nature.settings.dict_aux_operators`.
+    def second_q_ops(self) -> tuple[SecondQuantizedOp, dict[str, SecondQuantizedOp]]:
+        """Returns the second quantized operators associated with this problem.
 
         Returns:
-            A `list` or `dict` of `SecondQuantizedOp` objects.
+            A tuple, with the first object being the main operator and the second being a dictionary
+            of auxiliary operators.
         """
         driver_result = self.driver.run()
 
@@ -94,8 +90,9 @@ def second_q_ops(self) -> dict[str, SecondQuantizedOp]:
         self._grouped_property_transformed.basis = basis
 
         second_quantized_ops = self._grouped_property_transformed.second_q_ops()
+        main_op = second_quantized_ops.pop(self._main_property_name)
 
-        return second_quantized_ops
+        return main_op, second_quantized_ops
 
     def hopping_qeom_ops(
         self,
diff --git a/test/second_q/algorithms/ground_state_solvers/test_advanced_ucc_variants.py b/test/second_q/algorithms/ground_state_solvers/test_advanced_ucc_variants.py
index e318fd3d34..171664bc6c 100644
--- a/test/second_q/algorithms/ground_state_solvers/test_advanced_ucc_variants.py
+++ b/test/second_q/algorithms/ground_state_solvers/test_advanced_ucc_variants.py
@@ -53,10 +53,9 @@ def setUp(self):
 
         # because we create the initial state and ansatzes early, we need to ensure the qubit
         # converter already ran such that convert_match works as expected
+        main_op, _ = self.electronic_structure_problem.second_q_ops()
         _ = self.qubit_converter.convert(
-            self.electronic_structure_problem.second_q_ops()[
-                self.electronic_structure_problem.main_property_name
-            ],
+            main_op,
             self.num_particles,
         )
 
diff --git a/test/second_q/algorithms/ground_state_solvers/test_groundstate_eigensolver.py b/test/second_q/algorithms/ground_state_solvers/test_groundstate_eigensolver.py
index 45bd789841..bf367a9a7e 100644
--- a/test/second_q/algorithms/ground_state_solvers/test_groundstate_eigensolver.py
+++ b/test/second_q/algorithms/ground_state_solvers/test_groundstate_eigensolver.py
@@ -157,9 +157,7 @@ def _setup_evaluation_operators(self):
 
         # now we decide that we want to evaluate another operator
         # for testing simplicity, we just use some pre-constructed auxiliary operators
-        second_q_ops = self.electronic_structure_problem.second_q_ops()
-        # Remove main op to leave just aux ops
-        second_q_ops.pop(self.electronic_structure_problem.main_property_name)
+        _, second_q_ops = self.electronic_structure_problem.second_q_ops()
         aux_ops_dict = self.qubit_converter.convert_match(second_q_ops)
         return calc, res, aux_ops_dict
 
@@ -276,9 +274,7 @@ def test_eval_op_qasm(self):
         calc = GroundStateEigensolver(self.qubit_converter, solver)
         res_qasm = calc.solve(self.electronic_structure_problem)
 
-        hamiltonian = self.electronic_structure_problem.second_q_ops()[
-            self.electronic_structure_problem.main_property_name
-        ]
+        hamiltonian, _ = self.electronic_structure_problem.second_q_ops()
         qubit_op = self.qubit_converter.map(hamiltonian)
 
         ansatz = solver.get_solver(self.electronic_structure_problem, self.qubit_converter).ansatz
@@ -306,9 +302,7 @@ def test_eval_op_qasm_aer(self):
         calc = GroundStateEigensolver(self.qubit_converter, solver)
         res_qasm = calc.solve(self.electronic_structure_problem)
 
-        hamiltonian = self.electronic_structure_problem.second_q_ops()[
-            self.electronic_structure_problem.main_property_name
-        ]
+        hamiltonian, _ = self.electronic_structure_problem.second_q_ops()
         qubit_op = self.qubit_converter.map(hamiltonian)
 
         ansatz = solver.get_solver(self.electronic_structure_problem, self.qubit_converter).ansatz
diff --git a/test/second_q/mappers/test_qubit_converter.py b/test/second_q/mappers/test_qubit_converter.py
index cc2bd5660e..2fa9568e93 100644
--- a/test/second_q/mappers/test_qubit_converter.py
+++ b/test/second_q/mappers/test_qubit_converter.py
@@ -271,8 +271,9 @@ def test_molecular_problem_sector_locator_z2_symmetry(self):
 
         mapper = JordanWignerMapper()
         qubit_conv = QubitConverter(mapper, two_qubit_reduction=True, z2symmetry_reduction="auto")
+        main_op, _ = problem.second_q_ops()
         qubit_op = qubit_conv.convert(
-            problem.second_q_ops()[problem.main_property_name],
+            main_op,
             self.num_particles,
             sector_locator=problem.symmetry_sector_locator,
         )
diff --git a/test/second_q/problems/test_electronic_structure_problem.py b/test/second_q/problems/test_electronic_structure_problem.py
index 9d606656d5..ac195f65ad 100644
--- a/test/second_q/problems/test_electronic_structure_problem.py
+++ b/test/second_q/problems/test_electronic_structure_problem.py
@@ -43,7 +43,7 @@ class TestElectronicStructureProblem(QiskitNatureTestCase):
     def test_second_q_ops_without_transformers(self):
         """Tests that the list of second quantized operators is created if no transformers
         provided."""
-        expected_num_of_sec_quant_ops = 7
+        expected_num_of_sec_quant_ops = 6
         expected_fermionic_op_path = self.get_resource_path(
             "H2_631g_ferm_op_two_ints",
             "second_q/problems/resources",
@@ -55,9 +55,7 @@ def test_second_q_ops_without_transformers(self):
         )
         electronic_structure_problem = ElectronicStructureProblem(driver)
 
-        second_quantized_ops = electronic_structure_problem.second_q_ops()
-        electr_sec_quant_op = second_quantized_ops[electronic_structure_problem.main_property_name]
-        second_quantized_ops = list(second_quantized_ops.values())
+        electr_sec_quant_op, second_quantized_ops = electronic_structure_problem.second_q_ops()
 
         with self.subTest("Check that the correct properties aren't None"):
             # properties should never be None
@@ -67,7 +65,7 @@ def test_second_q_ops_without_transformers(self):
         with self.subTest("Check expected length of the list of second quantized operators."):
             assert len(second_quantized_ops) == expected_num_of_sec_quant_ops
         with self.subTest("Check types in the list of second quantized operators."):
-            for second_quantized_op in second_quantized_ops:
+            for second_quantized_op in second_quantized_ops.values():
                 assert isinstance(second_quantized_op, SecondQuantizedOp)
         with self.subTest("Check components of electronic second quantized operator."):
             assert all(
@@ -78,7 +76,7 @@ def test_second_q_ops_without_transformers(self):
     def test_second_q_ops_with_active_space(self):
         """Tests that the correct second quantized operator is created if an active space
         transformer is provided."""
-        expected_num_of_sec_quant_ops = 7
+        expected_num_of_sec_quant_ops = 6
         expected_fermionic_op_path = self.get_resource_path(
             "H2_631g_ferm_op_active_space",
             "second_q/problems/resources",
@@ -90,9 +88,7 @@ def test_second_q_ops_with_active_space(self):
         trafo = ActiveSpaceTransformer(num_electrons=2, num_molecular_orbitals=2)
 
         electronic_structure_problem = ElectronicStructureProblem(driver, [trafo])
-        second_quantized_ops = electronic_structure_problem.second_q_ops()
-        electr_sec_quant_op = second_quantized_ops[electronic_structure_problem.main_property_name]
-        second_quantized_ops = list(second_quantized_ops.values())
+        electr_sec_quant_op, second_quantized_ops = electronic_structure_problem.second_q_ops()
 
         with self.subTest("Check that the correct properties aren't None"):
             # properties should never be None
@@ -102,7 +98,7 @@ def test_second_q_ops_with_active_space(self):
         with self.subTest("Check expected length of the list of second quantized operators."):
             assert len(second_quantized_ops) == expected_num_of_sec_quant_ops
         with self.subTest("Check types in the list of second quantized operators."):
-            for second_quantized_op in second_quantized_ops:
+            for second_quantized_op in second_quantized_ops.values():
                 assert isinstance(second_quantized_op, SecondQuantizedOp)
         with self.subTest("Check components of electronic second quantized operator."):
             assert all(
@@ -125,8 +121,9 @@ def test_sector_locator_h2o(self):
         qubit_conv = QubitConverter(
             mapper=ParityMapper(), two_qubit_reduction=True, z2symmetry_reduction="auto"
         )
+        main_op, _ = es_problem.second_q_ops()
         qubit_conv.convert(
-            es_problem.second_q_ops()[es_problem.main_property_name],
+            main_op,
             num_particles=es_problem.num_particles,
             sector_locator=es_problem.symmetry_sector_locator,
         )
@@ -146,8 +143,9 @@ def test_sector_locator_homonuclear(self):
         qubit_conv = QubitConverter(
             mapper=ParityMapper(), two_qubit_reduction=True, z2symmetry_reduction="auto"
         )
+        main_op, _ = es_problem.second_q_ops()
         qubit_conv.convert(
-            es_problem.second_q_ops()[es_problem.main_property_name],
+            main_op,
             num_particles=es_problem.num_particles,
             sector_locator=es_problem.symmetry_sector_locator,
         )
diff --git a/test/second_q/problems/test_lattice_model_problem.py b/test/second_q/problems/test_lattice_model_problem.py
index 188071a247..896d50dc93 100644
--- a/test/second_q/problems/test_lattice_model_problem.py
+++ b/test/second_q/problems/test_lattice_model_problem.py
@@ -48,8 +48,10 @@ def test_second_q_ops(self):
         boundary_condition = BoundaryCondition.OPEN
         line_lattice = LineLattice(num_nodes=num_nodes, boundary_condition=boundary_condition)
         fhm = FermiHubbardModel(lattice=line_lattice, onsite_interaction=5.0)
+        expected_op = fhm.second_q_ops()
         lmp = LatticeModelProblem(fhm)
-        self._compare_second_q_op(fhm.second_q_ops(), lmp.second_q_ops()[lmp.main_property_name])
+        main_op, _ = lmp.second_q_ops()
+        self._compare_second_q_op(expected_op, main_op)
 
     def test_interpret(self):
         """Tests that the result is interpreted"""
diff --git a/test/second_q/problems/test_vibrational_structure_problem.py b/test/second_q/problems/test_vibrational_structure_problem.py
index 86ec3dddc6..13c54dbd01 100644
--- a/test/second_q/problems/test_vibrational_structure_problem.py
+++ b/test/second_q/problems/test_vibrational_structure_problem.py
@@ -53,15 +53,14 @@ def setUp(self) -> None:
     def test_second_q_ops_without_transformers(self):
         """Tests that the list of second quantized operators is created if no transformers
         provided."""
-        expected_num_of_sec_quant_ops = 5
+        expected_num_of_sec_quant_ops = 4
         expected_len_of_vibrational_op = 47
         num_modals = 2
         truncation_order = 3
         num_modes = self.props.num_modes
         num_modals = [num_modals] * num_modes
         vibrational_problem = VibrationalStructureProblem(self.driver, num_modals, truncation_order)
-        second_quantized_ops = vibrational_problem.second_q_ops()
-        vibrational_op = second_quantized_ops[vibrational_problem.main_property_name]
+        vibrational_op, second_quantized_ops = vibrational_problem.second_q_ops()
 
         with self.subTest("Check expected length of the list of second quantized operators."):
             assert len(second_quantized_ops) == expected_num_of_sec_quant_ops
@@ -73,15 +72,14 @@ def test_second_q_ops_without_transformers(self):
 
     def test_truncation_order(self):
         """Tests that the truncation_order is being respected."""
-        expected_num_of_sec_quant_ops = 5
+        expected_num_of_sec_quant_ops = 4
         expected_len_of_vibrational_op = 10
         num_modals = 2
         truncation_order = 1
         num_modes = self.props.num_modes
         num_modals = [num_modals] * num_modes
         vibrational_problem = VibrationalStructureProblem(self.driver, num_modals, truncation_order)
-        second_quantized_ops = vibrational_problem.second_q_ops()
-        vibrational_op = second_quantized_ops[vibrational_problem.main_property_name]
+        vibrational_op, second_quantized_ops = vibrational_problem.second_q_ops()
 
         with self.subTest("Check expected length of the list of second quantized operators."):
             assert len(second_quantized_ops) == expected_num_of_sec_quant_ops
diff --git a/test/test_end2end_with_vqe.py b/test/test_end2end_with_vqe.py
index 7e698c2c33..6f6de2a944 100644
--- a/test/test_end2end_with_vqe.py
+++ b/test/test_end2end_with_vqe.py
@@ -37,14 +37,14 @@ def setUp(self):
             hdf5_input=self.get_resource_path("test_driver_hdf5.hdf5", "second_q/drivers/hdf5d")
         )
         problem = ElectronicStructureProblem(driver)
-        second_q_ops = [problem.second_q_ops()[problem.main_property_name]]
+        main_op, aux_ops = problem.second_q_ops()
         converter = QubitConverter(mapper=ParityMapper(), two_qubit_reduction=True)
         num_particles = (
             problem.grouped_property_transformed.get_property("ParticleNumber").num_alpha,
             problem.grouped_property_transformed.get_property("ParticleNumber").num_beta,
         )
-        self.qubit_op = converter.convert(second_q_ops[0], num_particles)
-        self.aux_ops = converter.convert_match(second_q_ops[1:])
+        self.qubit_op = converter.convert(main_op, num_particles)
+        self.aux_ops = converter.convert_match(aux_ops)
         self.reference_energy = -1.857275027031588
 
     def test_end2end_h2(self):