diff --git a/tangelo/algorithms/projective/__init__.py b/tangelo/algorithms/projective/__init__.py
index 6db96b66b..e83d338e3 100644
--- a/tangelo/algorithms/projective/__init__.py
+++ b/tangelo/algorithms/projective/__init__.py
@@ -13,3 +13,5 @@
 # limitations under the License.
 
 from .quantum_imaginary_time import QITESolver
+from .qpe import QPESolver
+from .iqpe import IterativeQPESolver
diff --git a/tangelo/algorithms/projective/iqpe.py b/tangelo/algorithms/projective/iqpe.py
new file mode 100644
index 000000000..03ae9f00e
--- /dev/null
+++ b/tangelo/algorithms/projective/iqpe.py
@@ -0,0 +1,297 @@
+# Copyright 2023 Good Chemistry Company.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Implements the iterative Quantum Phase Estimation (iQPE) algorithm to solve
+electronic structure calculations.
+
+Ref:
+    M. Dobsicek, G. Johansson, V. Shumeiko, and G. Wendin, Arbitrary Accuracy Iterative Quantum Phase
+    Estimation Algorithm using a Single Ancillary Qubit: A two-qubit benchmark, Phys. Rev. A 76, 030306(R)
+    (2007).
+"""
+
+from typing import Optional, Union, List
+from collections import Counter
+from enum import Enum
+
+import numpy as np
+
+from tangelo import SecondQuantizedMolecule
+from tangelo.linq import get_backend, Circuit, Gate, ClassicalControl, generate_applied_gates
+from tangelo.toolboxes.operators import QubitOperator
+from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
+from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_mapped_vector, vector_to_circuit, get_reference_circuit
+import tangelo.toolboxes.unitary_generator as ugen
+import tangelo.toolboxes.ansatz_generator as agen
+from tangelo.toolboxes.post_processing.histogram import Histogram
+
+
+class BuiltInUnitary(Enum):
+    """Enumeration of the ansatz circuits supported by iQPE."""
+    TrotterSuzuki = ugen.TrotterSuzukiUnitary
+    CircuitUnitary = ugen.CircuitUnitary
+
+
+class IterativeQPESolver:
+    r"""Solve the electronic structure problem for a molecular system by using
+    the iterative Quantum Phase Estimation (iQPE) algorithm.
+
+    This algorithm evaluates the energy of a molecular system by performing
+    a series of controlled time-evolutions
+
+    Users must first set the desired options of the iterative QPESolver object through the
+    __init__ method, and call the "build" method to build the underlying objects
+    (mean-field, hardware backend, unitary...). They are then able to call the
+    the simulate method. In particular, simulate
+    runs the iQPE algorithm, returning the optimal energy found by the most probable
+    measurement as a binary fraction.
+
+    Attributes:
+        molecule (SecondQuantizedMolecule) : The molecular system.
+        qubit_mapping (str) : one of the supported qubit mapping identifiers.
+        unitary (Unitary) : one of the supported unitary evolutions.
+        backend_options (dict): parameters to build the underlying compute backend (simulator, etc).
+        simulate_options (dict): Options for fine-control of the simulator backend, including desired measurement results, etc.
+        penalty_terms (dict): parameters for penalty terms to append to target
+            qubit Hamiltonian (see penalty_terms for more details).
+        unitary_options (dict): parameters for the given ansatz (see given ansatz
+            file for details).
+        up_then_down (bool): change basis ordering putting all spin up orbitals
+            first, followed by all spin down. Default, False has alternating
+                spin up/down ordering.
+        qubit_hamiltonian (QubitOperator): The Hamiltonian expressed as a sum of products of Pauli matrices.
+        verbose (bool): Flag for iQPE verbosity.
+        projective_circuit (Circuit): A terminal circuit that projects into the correct space, always added to
+            the end of the unitary circuit. Could be measurement gates for example
+        ref_state (array or Circuit): The reference configuration to use. Replaces HF state
+        size_qpe_register (int): The number of iterations of single qubit iQPE to use for the calculation.
+    """
+
+    def __init__(self, opt_dict):
+
+        default_backend_options = {"target": None, "n_shots": 1, "noise_model": None}
+        copt_dict = opt_dict.copy()
+
+        self.molecule: Optional[SecondQuantizedMolecule] = copt_dict.pop("molecule", None)
+        self.qubit_mapping: str = copt_dict.pop("qubit_mapping", "jw")
+        self.unitary: ugen.Unitary = copt_dict.pop("unitary", BuiltInUnitary.TrotterSuzuki)
+        self.backend_options: dict = copt_dict.pop("backend_options", default_backend_options)
+        self.penalty_terms: Optional[dict] = copt_dict.pop("penalty_terms", None)
+        self.simulate_options: dict = copt_dict.pop("simulate_options", dict())
+        self.unitary_options: dict = copt_dict.pop("unitary_options", dict())
+        self.up_then_down: bool = copt_dict.pop("up_then_down", False)
+        self.qubit_hamiltonian: QubitOperator = copt_dict.pop("qubit_hamiltonian", None)
+        self.verbose: bool = copt_dict.pop("verbose", False)
+        self.projective_circuit: Circuit = copt_dict.pop("projective_circuit", None)
+        self.ref_state: Optional[Union[list, Circuit]] = copt_dict.pop("ref_state", None)
+        self.n_qpe_qubits: int = copt_dict.pop("size_qpe_register", 1)
+
+        if len(copt_dict) > 0:
+            raise KeyError(f"The following keywords are not supported in {self.__class__.__name__}: \n {copt_dict.keys()}")
+
+        # If nothing is provided raise an Error:
+        if not (bool(self.molecule) or bool(self.qubit_hamiltonian) or isinstance(self.unitary, (ugen.Unitary, Circuit))):
+            raise ValueError(f"A Molecule or a QubitOperator or a Unitary object/Circuit must be provided in {self.__class__.__name__}")
+
+        # Raise error/warnings if input is not as expected. Only a single input
+        if (bool(self.molecule) and bool(self.qubit_hamiltonian)):
+            raise ValueError(f"Incompatible Options in {self.__class__.__name__}:"
+                             "Only one of the following can be provided by user: molecule OR qubit Hamiltonian.")
+        if isinstance(self.unitary, (Circuit, ugen.Unitary)) and bool(self.qubit_hamiltonian):
+            raise ValueError(f"Incompatible Options in {self.__class__.__name__}:"
+                             "Only one of the following can be provided by user: unitary OR qubit Hamiltonian.")
+        if isinstance(self.unitary, (Circuit, ugen.Unitary)) and bool(self.molecule):
+            raise Warning("The molecule is only being used to generate the reference state. The unitary is being used for the iQPE.")
+
+        # Initialize the reference state circuit.
+        if self.ref_state is not None:
+            if isinstance(self.ref_state, Circuit):
+                self.reference_circuit = self.ref_state
+            else:
+                self.reference_circuit = vector_to_circuit(get_mapped_vector(self.ref_state, self.qubit_mapping, self.up_then_down))
+        else:
+            if bool(self.molecule):
+                self.reference_circuit = get_reference_circuit(self.molecule.n_active_sos,
+                                                               self.molecule.n_active_electrons,
+                                                               self.qubit_mapping,
+                                                               self.up_then_down,
+                                                               self.molecule.spin)
+            else:
+                self.reference_circuit = Circuit()
+
+        default_backend_options.update(self.backend_options)
+        self.backend_options = default_backend_options
+        self.builtin_unitary = set(BuiltInUnitary)
+
+    def build(self):
+        """Build the underlying objects required to run the iQPE algorithm
+        afterwards.
+        """
+
+        if isinstance(self.unitary, Circuit):
+            self.unitary = ugen.CircuitUnitary(self.unitary, **self.unitary_options)
+
+        # Building QPE with a molecule as input.
+        if self.molecule:
+
+            # Compute qubit hamiltonian for the input molecular system
+            self.qubit_hamiltonian = fermion_to_qubit_mapping(fermion_operator=self.molecule.fermionic_hamiltonian,
+                                                              mapping=self.qubit_mapping,
+                                                              n_spinorbitals=self.molecule.n_active_sos,
+                                                              n_electrons=self.molecule.n_active_electrons,
+                                                              up_then_down=self.up_then_down,
+                                                              spin=self.molecule.active_spin)
+
+            if self.penalty_terms:
+                pen_ferm = agen.penalty_terms.combined_penalty(self.molecule.n_active_mos, self.penalty_terms)
+                pen_qubit = fermion_to_qubit_mapping(fermion_operator=pen_ferm,
+                                                     mapping=self.qubit_mapping,
+                                                     n_spinorbitals=self.molecule.n_active_sos,
+                                                     n_electrons=self.molecule.n_active_electrons,
+                                                     up_then_down=self.up_then_down,
+                                                     spin=self.molecule.active_spin)
+                self.qubit_hamiltonian += pen_qubit
+
+        if isinstance(self.unitary, BuiltInUnitary):
+            self.unitary = self.unitary.value(self.qubit_hamiltonian, **self.unitary_options)
+        elif not isinstance(self.unitary, ugen.Unitary):
+            raise TypeError("Invalid ansatz dataype. Expecting a custom Unitary (Unitary class).")
+
+        # Quantum circuit simulation backend options
+        self.backend = get_backend(**self.backend_options)
+
+        # Determine where to place QPE ancilla qubit index
+        self.n_state, self.n_ancilla = self.unitary.qubit_indices()
+        self.qft_qubit = max(list(self.n_state)+list(self.n_ancilla)) + 1
+
+        self.cfunc = IterativeQPEControl(self.n_qpe_qubits, self.qft_qubit, self.unitary)
+        self.circuit = Circuit(self.reference_circuit._gates+[Gate("CMEASURE", self.qft_qubit)],
+                               cmeasure_control=self.cfunc, n_qubits=self.qft_qubit+1)
+
+    def simulate(self):
+        """Run the iQPE circuit. Return the energy of the most probable bitstring
+
+        Attributes:
+            bitstring (str): The most probable bitstring.
+            histogram (Histogram): The full Histogram of measurements on the iQPE ancilla qubit representing energies.
+            qpe_freqs (dict): The dictionary of measurements on the iQPE ancilla qubit.
+            freqs (dict): The full dictionary of measurements on all qubits.
+
+        Returns:
+            float: The energy of the most probable bitstring measured during iQPE.
+        """
+
+        if not (self.unitary and self.backend):
+            raise RuntimeError(f"No unitary or hardware backend built. Have you called {self.__class__.__name__}.build ?")
+
+        self.freqs, _ = self.backend.simulate(self.circuit)
+        self.histogram = Histogram(self.freqs)
+        self.histogram.remove_qubit_indices(*(self.n_state+self.n_ancilla))
+        qpe_counts = Counter(self.cfunc.measurements[:self.backend.n_shots])
+        self.qpe_freqs = {key[::-1]: value/self.backend.n_shots for key, value in qpe_counts.items()}
+        self.bitstring = max(self.qpe_freqs.items(), key=lambda x: x[1])[0]
+
+        return self.energy_estimation(self.bitstring)
+
+    def get_resources(self):
+        """Estimate the resources required by iQPE, with the current unitary. This
+        assumes "build" has been run, as it requires the circuit and the
+        qubit Hamiltonian. Return information that pertains to the user, for the
+        purpose of running an experiment on a classical simulator or a quantum
+        device.
+        """
+
+        resources = dict()
+
+        # If the attribute of the applied_gates has been populated, use the exact resources else approximate the iQPE circuit
+        circuit = Circuit(self.circuit.applied_gates) if self.circuit.applied_gates else Circuit(generate_applied_gates(self.circuit))
+        resources["applied_circuit_width"] = circuit.width
+        resources["applied_circuit_depth"] = circuit.depth()
+        resources["applied_circuit_2qubit_gates"] = circuit.counts_n_qubit.get(2, 0)
+        return resources
+
+    def energy_estimation(self, bitstring):
+        """Estimate energy using the calculated frequency dictionary.
+
+        Args:
+            bitstring (str): The bitstring to evaluate the energy of in base 2.
+
+        Returns:
+             float: Energy of the given bitstring
+        """
+
+        return sum([0.5**(i+1) for i, b in enumerate(bitstring) if b == "1"])
+
+
+class IterativeQPEControl(ClassicalControl):
+    def __init__(self, n_bits: int, qft_qubit: int, u: ugen.Unitary):
+        """Iterative QPE with n_bits"""
+        self.n_bits: int = n_bits
+        self.bitplace: int = n_bits
+        self.phase: float = 0
+        self.measurements: List[str] = [""]
+        self.energies: List[float] = [0.]
+        self.n_runs: int = 0
+        self.qft_qubit: int = qft_qubit
+        self.unitary: ugen.Unitary = u
+        self.started: bool = False
+
+    def return_gates(self, measurement) -> List[Gate]:
+        """Return a list of gates based on the measurement outcome for the next step in iQPE.
+
+        Each measurement updates the current phase correction and returns a list of gates that
+        implements the next controlled time-evolution along with the phase correction to determine
+        the next bit value.
+
+        Args:
+            measurement (str): "1" or "0"
+
+        Returns:
+            List[Gate]: The next gates to apply to the circuit
+        """
+        # Ignore the first measurement as it is always 0 and meaningless.
+        if self.started:
+            self.measurements[self.n_runs] += measurement
+            self.energies[self.n_runs] += int(measurement)/2**self.bitplace
+        else:
+            self.started = True
+
+        if self.bitplace > 0:
+            # Update phase and determine reset gates
+            if measurement == "1":
+                self.phase += 1/2**(self.bitplace)
+                reset_to_zero = [Gate("X", self.qft_qubit)]
+            else:
+                reset_to_zero = []
+
+            # Decrease bitplace and apply next phase estimation
+            self.bitplace -= 1
+            phase_correction = [Gate("PHASE", self.qft_qubit, parameter=-np.pi*self.phase*2**(self.bitplace))]
+            gates = reset_to_zero + [Gate("H", self.qft_qubit)] + phase_correction
+            gates += self.unitary.build_circuit(2**(self.bitplace), self.qft_qubit)._gates + [Gate("H", self.qft_qubit)]
+            return gates + [Gate("CMEASURE", self.qft_qubit)]
+        else:
+            return []
+
+    def finalize(self):
+        """Called from simulator after all gates have been called.
+
+        Reinitialize all variables and store measurements and energies from the current iQPE run.
+        """
+        self.bitplace = self.n_bits
+        self.phase = 0
+        self.n_runs += 1
+        self.measurements += [""]
+        self.energies += [0.]
+        self.started = False
diff --git a/tangelo/algorithms/projective/qpe.py b/tangelo/algorithms/projective/qpe.py
index db78fc208..0a754a7f7 100644
--- a/tangelo/algorithms/projective/qpe.py
+++ b/tangelo/algorithms/projective/qpe.py
@@ -94,10 +94,17 @@ def __init__(self, opt_dict):
         if len(copt_dict) > 0:
             raise KeyError(f"The following keywords are not supported in {self.__class__.__name__}: \n {copt_dict.keys()}")
 
+        # If nothing is provided raise and Error:
+        if not (bool(self.molecule) or bool(self.qubit_hamiltonian) or isinstance(self.unitary, (unitary.Unitary, Circuit))):
+            raise ValueError(f"A Molecule or a QubitOperator or a Unitary object/Circuit must be provided in {self.__class__.__name__}")
+
         # Raise error/warnings if input is not as expected. Only a single input
-        # must be provided to avoid conflicts.
-        if not (bool(self.molecule) ^ bool(self.qubit_hamiltonian)):
-            raise ValueError(f"A molecule OR qubit Hamiltonian object must be provided when instantiating {self.__class__.__name__}.")
+        if (bool(self.molecule) and bool(self.qubit_hamiltonian)):
+            raise ValueError(f"Both a molecule and qubit Hamiltonian can not be provided when instantiating {self.__class__.__name__}.")
+        if isinstance(self.unitary, Circuit) and bool(self.qubit_hamiltonian):
+            raise ValueError(f"Both a qubit Hamiltonian and a circuit defining the unitary can not be provided in {self.__class__.__name__}.")
+        if isinstance(self.unitary, (Circuit, unitary.Unitary)) and bool(self.molecule):
+            raise Warning("The molecule is only being used to generate the reference state. The unitary is being used for the QPE.")
 
         if self.ref_state is not None:
             if isinstance(self.ref_state, Circuit):
@@ -150,7 +157,7 @@ def build(self):
         if isinstance(self.unitary, BuiltInUnitary):
             self.unitary = self.unitary.value(self.qubit_hamiltonian, **self.unitary_options)
         elif not isinstance(self.unitary, unitary.Unitary):
-            raise TypeError(f"Invalid ansatz dataype. Expecting a custom Unitary (Unitary class).")
+            raise TypeError("Invalid ansatz dataype. Expecting a custom Unitary (Unitary class).")
 
         # Quantum circuit simulation backend options
         self.backend = get_backend(**self.backend_options)
diff --git a/tangelo/algorithms/projective/tests/test_iqpe.py b/tangelo/algorithms/projective/tests/test_iqpe.py
new file mode 100644
index 000000000..e6fc4c92b
--- /dev/null
+++ b/tangelo/algorithms/projective/tests/test_iqpe.py
@@ -0,0 +1,138 @@
+# Copyright 2023 Good Chemsitry Company.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+
+import numpy as np
+from openfermion import get_sparse_operator
+
+from tangelo.helpers.utils import installed_backends
+from tangelo.algorithms.projective.iqpe import IterativeQPESolver
+from tangelo.toolboxes.ansatz_generator.ansatz_utils import trotterize
+from tangelo.toolboxes.operators import QubitOperator
+from tangelo.molecule_library import mol_H2_sto3g
+from tangelo.linq.helpers.circuits.statevector import StateVector
+
+
+class IterativeQPESolverTest(unittest.TestCase):
+
+    def test_instantiation(self):
+        """Try instantiating IterativeQPESolver with basic input."""
+
+        options = {"molecule": mol_H2_sto3g, "qubit_mapping": "jw"}
+        IterativeQPESolver(options)
+
+    def test_instantiation_incorrect_keyword(self):
+        """Instantiating with an incorrect keyword should return an error """
+
+        options = {"molecule": mol_H2_sto3g, "qubit_mapping": "jw", "dummy": True}
+        self.assertRaises(KeyError, IterativeQPESolver, options)
+
+    def test_instantiation_missing_molecule(self):
+        """Instantiating with no molecule should return an error."""
+
+        options = {"qubit_mapping": "jw"}
+        self.assertRaises(ValueError, IterativeQPESolver, options)
+
+    @unittest.skipIf("qulacs" not in installed_backends, "Test Skipped: Backend not available \n")
+    def test_simulate_h2(self):
+        """Run iQPE on H2 molecule, with scbk qubit mapping and exact simulator with the approximate initial state
+        """
+
+        iqpe_options = {"molecule": mol_H2_sto3g, "qubit_mapping": "scbk", "up_then_down": True, "size_qpe_register": 7,
+                       "backend_options": {"target": "qulacs", "n_shots": 20}, "unitary_options": {"time": 2*np.pi, "n_trotter_steps": 1,
+                                                                                    "n_steps_method": "repeat", "trotter_order": 4}}
+        iqpe_solver = IterativeQPESolver(iqpe_options)
+        iqpe_solver.build()
+
+        iqpe_solver.simulate()
+
+        # Use the highest probability circuit which is about 0.5. Will fail ~1 in every 2^20 times.
+        max_prob_key = max(iqpe_solver.circuit.success_probabilities, key=iqpe_solver.circuit.success_probabilities.get)
+        self.assertAlmostEqual(iqpe_solver.energy_estimation(max_prob_key[::-1]), 0.14, delta=1e-2)
+
+    @unittest.skipIf("qulacs" not in installed_backends, "Test Skipped: Backend not available \n")
+    def test_circuit_input(self):
+        """Run iQPE on a qubit Hamiltonian, providing only the Trotter circuit and the exact initial state.
+        """
+
+        # Generate qubit operator with state 9 having eigenvalue 0.25
+        qu_op = (QubitOperator("X0 X1", 0.125) + QubitOperator("Y1 Y2", 0.125) + QubitOperator("Z2 Z3", 0.125))
+
+        ham_mat = get_sparse_operator(qu_op.to_openfermion()).toarray()
+        e, wavefunction = np.linalg.eigh(ham_mat)
+
+        sv = StateVector(wavefunction[:, 9], order="lsq_first")
+        ref_circ = sv.initializing_circuit()
+
+        unit_circ = trotterize(qu_op, -2*np.pi, 2, 4, True)
+
+        # Test supplying circuit and applying iQPE controls to only gates marked as variational
+        iqpe_options = {"unitary": unit_circ, "size_qpe_register": 3, "ref_state": ref_circ,
+                       "backend_options": {"target": "qulacs", "n_shots": 1}, "unitary_options": {"control_method": "variational"}}
+        iqpe_solver = IterativeQPESolver(iqpe_options)
+        iqpe_solver.build()
+
+        energy = iqpe_solver.simulate()
+
+        self.assertAlmostEqual(energy, 0.125, delta=1e-3)
+
+        # Test supplying circuit with iQPE controls added to every gate.
+        iqpe_options = {"unitary": unit_circ, "size_qpe_register": 3, "ref_state": ref_circ,
+                       "backend_options": {"target": "qulacs", "n_shots": 1}, "unitary_options": {"control_method": "all"}}
+        iqpe_solver = IterativeQPESolver(iqpe_options)
+        iqpe_solver.build()
+
+        energy = iqpe_solver.simulate()
+
+        self.assertAlmostEqual(energy, 0.125, delta=1e-3)
+
+    def test_qubit_hamiltonian_input(self):
+        """Test with qubit hamiltonian input."""
+
+        # Generate qubit operator with state 9 having eigenvalue 0.25
+        qu_op = (QubitOperator("X0 X1", 0.125) + QubitOperator("Y1 Y2", 0.125) + QubitOperator("Z2 Z3", 0.125)
+                 + QubitOperator("", 0.125))
+
+        ham_mat = get_sparse_operator(qu_op.to_openfermion()).toarray()
+        _, wavefunction = np.linalg.eigh(ham_mat)
+
+        sv = StateVector(wavefunction[:, 9], order="lsq_first")
+        init_circ = sv.initializing_circuit()
+        init_circ.simplify()
+
+        iqpe = IterativeQPESolver({"qubit_hamiltonian": qu_op, "size_qpe_register": 2, "ref_state": init_circ,
+                                  "backend_options": {"noise_model": None, "target": "cirq"},
+                                  "unitary_options": {"time": -2*np.pi, "n_trotter_steps": 1,
+                                                      "n_steps_method": "repeat", "trotter_order": 4}})
+        iqpe.build()
+
+        # Test that get_resources returns expected results before running the simulation.
+        # And the number of two qubit gates is the same as after the simulation.
+        resources = iqpe.get_resources()
+        self.assertEqual(resources["applied_circuit_width"], 5)
+        two_qubit_gates_before = resources["applied_circuit_2qubit_gates"]
+
+        # Simulate and obtain energy estimation.
+        energy = iqpe.simulate()
+        self.assertAlmostEqual(energy, 0.25, delta=1.e-5)
+
+        # Test that get_resources returns expected results after running the simulation.
+        resources = iqpe.get_resources()
+        self.assertEqual(resources["applied_circuit_width"], 5)
+        self.assertEqual(resources["applied_circuit_2qubit_gates"], two_qubit_gates_before)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tangelo/algorithms/projective/tests/test_qpe.py b/tangelo/algorithms/projective/tests/test_qpe.py
index 5925ecfba..63d93ea8b 100644
--- a/tangelo/algorithms/projective/tests/test_qpe.py
+++ b/tangelo/algorithms/projective/tests/test_qpe.py
@@ -76,7 +76,7 @@ def test_simulate_h2_circuit(self):
                                 "n_electrons": mol_H2_sto3g.n_active_electrons})
 
         # Test supplying circuit and applying QPE controls to only gates marked as variational
-        qpe_options = {"qubit_hamiltonian": qu_op, "unitary": unit_circ, "size_qpe_register": 8, "ref_state": ref_circ,
+        qpe_options = {"unitary": unit_circ, "size_qpe_register": 8, "ref_state": ref_circ,
                        "backend_options": {"target": "qulacs"}, "unitary_options": {"control_method": "variational"}}
         qpe_solver = QPESolver(qpe_options)
         qpe_solver.build()
@@ -86,7 +86,7 @@ def test_simulate_h2_circuit(self):
         self.assertAlmostEqual(energy, -(-1.13727-qu_op.constant), delta=1e-3)
 
         # Test supplying circuit with QPE controls added to every gate.
-        qpe_options = {"qubit_hamiltonian": qu_op, "unitary": unit_circ, "size_qpe_register": 8, "ref_state": ref_circ,
+        qpe_options = {"unitary": unit_circ, "size_qpe_register": 8, "ref_state": ref_circ,
                        "backend_options": {"target": "qulacs"}, "unitary_options": {"control_method": "all"}}
         qpe_solver = QPESolver(qpe_options)
         qpe_solver.build()
diff --git a/tangelo/linq/__init__.py b/tangelo/linq/__init__.py
index a5f05f5a1..5a077d376 100644
--- a/tangelo/linq/__init__.py
+++ b/tangelo/linq/__init__.py
@@ -13,7 +13,7 @@
 # limitations under the License.
 
 from .gate import *
-from .circuit import Circuit, stack, remove_small_rotations, remove_redundant_gates, get_unitary_circuit_pieces, ClassicalControl
+from .circuit import Circuit, stack, remove_small_rotations, remove_redundant_gates, get_unitary_circuit_pieces, ClassicalControl, generate_applied_gates
 from .translator import *
 from .simulator import get_backend
 from .target.backend import get_expectation_value_from_frequencies_oneterm
diff --git a/tangelo/linq/circuit.py b/tangelo/linq/circuit.py
index 5c908fa43..d52235fc3 100644
--- a/tangelo/linq/circuit.py
+++ b/tangelo/linq/circuit.py
@@ -20,7 +20,8 @@
 
 import copy
 import abc
-from typing import List, Tuple, Iterator, Union, Set, Dict, Callable
+from typing import List, Tuple, Union, Set, Dict, Callable
+import warnings
 
 import numpy as np
 from cirq.contrib.svg import SVGCircuit
@@ -189,7 +190,7 @@ def draw(self):
 
     def copy(self):
         """Return a deepcopy of circuit"""
-        return Circuit(copy.deepcopy(self._gates), n_qubits=self._qubits_simulated, name=self.name)
+        return Circuit(copy.deepcopy(self._gates), n_qubits=self._qubits_simulated, name=self.name, cmeasure_control=copy.deepcopy(self._cmeasure_control))
 
     def add_gate(self, g):
         """Add a new gate to a circuit object and update other fields of the
@@ -278,7 +279,7 @@ def reindex_qubits(self, new_indices):
         """
 
         if len(new_indices) != len(self._qubit_indices):
-            raise ValueError(f"The number of indices does not match the length of self._qubit_indices")
+            raise ValueError("The number of indices does not match the length of self._qubit_indices")
 
         qubits_in_use = self._qubit_indices
         mapping = {i: j for i, j in zip(qubits_in_use, new_indices)}
@@ -417,7 +418,7 @@ def controlled_measurement_op(self, measure):
         if callable(self._cmeasure_control):
             return Circuit(self._cmeasure_control(measure), n_qubits=self.width)
         elif isinstance(self._cmeasure_control, ClassicalControl):
-            return Circuit(self._cmeasure_control.return_circuit(measure), n_qubits=self.width)
+            return Circuit(self._cmeasure_control.return_gates(measure), n_qubits=self.width)
         else:
             raise TypeError(f"cmeasure_control must either be a function or an instance of {ClassicalControl}")
 
@@ -652,24 +653,100 @@ def __init__(self):
         """instantiate classical control operations"""
 
     @abc.abstractmethod
-    def return_circuit(self, measurement) -> List[Gate]:
-        """Return a circuit based on the measurement outcome.
+    def return_gates(self, measurement) -> List[Gate]:
+        """Return the list of gates based on the measurement outcome.
 
         Args:
             measurement (str): "1" or "0"
             qubit (int): The qubit index
 
         Returns:
-            Circuit: The next circuit to apply
+            List[Gate]: The next gates to apply to the circuit
         """
 
     @abc.abstractmethod
-    def finalize(self, frequencies):
+    def finalize(self):
         """Called from simulator after all gates have been called.
 
         Can be used to reinitialize variables for the next run
-        and store results
-
-        Args:
-            frequencies (dict): The measured frequencies for the final state measurement.
+        and store results.
         """
+
+
+def generate_applied_gates(source_circuit: Circuit, desired_meas_result=None) -> List[Gate]:
+    """Generate the applied gates of a Circuit without explicitly simulating.
+
+    Useful to determine the resources required for a circuit with measurement controlled operations given certain
+    measurement outcomes.
+
+    Note: Measurement outcomes with zero probability can not be screened.
+
+    Args:
+        source_circuit (Circuit): A circuit in the abstract format to be simulated with
+            the classical control function called.
+        desired_meas_result (str): The binary string of the desired measurement.
+            Must have the same length as the number of CMEASURE+MEASURE gates in source_circuit
+    """
+
+    circuit = source_circuit.copy()
+    n_cmeas = circuit.counts.get("CMEASURE", 0)
+    if n_cmeas == 0:
+        warnings.warn("The supplied circuit does not contain CMEASURE gates."
+                      "This function will not modify the applied_gates attribute.")
+        return
+
+    applied_gates = []
+    dmeas = None if not desired_meas_result else list(desired_meas_result)
+    measurements = ""
+
+    # Break circuit into pieces that do not include CMEASURE or MEASURE gates
+    unitary_circuits, qubits, cmeasure_flags = get_unitary_circuit_pieces(circuit)
+    # Generate list of circuits that are extended by previous CMEASURE operations
+    precirc = [Circuit()]*len(unitary_circuits)
+
+    while len(unitary_circuits) > 1:
+        c = precirc[0]+unitary_circuits[0]
+        applied_gates += c._gates
+
+        # Perform measurement.
+        measure = dmeas[0] if desired_meas_result else "1"
+        measurements += measure
+        if desired_meas_result:
+            del dmeas[0]
+
+        # If a CMEASURE has occurred
+        if cmeasure_flags[0] is not None:
+            applied_gates += [Gate("CMEASURE", qubits[0], parameter=measure)]
+            if isinstance(cmeasure_flags[0], str):
+                newcirc = circuit.controlled_measurement_op(measure)
+            elif isinstance(cmeasure_flags[0], dict):
+                newcirc = Circuit(cmeasure_flags[0][measure], n_qubits=circuit.width)
+            new_unitary_circuits, new_qubits, new_cmeasure_flags = get_unitary_circuit_pieces(newcirc)
+
+        # No classical control
+        else:
+            applied_gates += [Gate("MEASURE", qubits[0], parameter=measure)]
+            new_unitary_circuits = [Circuit(n_qubits=circuit.width)]
+            new_qubits = []
+            new_cmeasure_flags = []
+
+        # Remove circuits, measurements and corresponding qubits that have been applied.
+        del unitary_circuits[0]
+        del qubits[0]
+        del cmeasure_flags[0]
+        del precirc[0]
+        precirc[0] = new_unitary_circuits[-1] + precirc[0]
+
+        # If new_unitary_circuits includes MEASURE or CMEASURE Gates, the number of unitary_circuits grows.
+        if len(new_unitary_circuits) > 1:
+            unitary_circuits = new_unitary_circuits[:-1] + unitary_circuits
+            qubits = new_qubits + qubits
+            cmeasure_flags = new_cmeasure_flags + cmeasure_flags
+            precirc = [Circuit()]*len(qubits) + precirc
+
+    # No more MEASURE or CMEASURE gates are present, run final unitary circuit segment and set attributes
+    final_circuit = precirc[0] + unitary_circuits[-1]
+    # Call the finalize method of ClassicalControl, used to reset variables, perform computation etc.
+    circuit.finalize_cmeasure_control()
+
+    return applied_gates + final_circuit._gates
diff --git a/tangelo/linq/target/backend.py b/tangelo/linq/target/backend.py
index b983d3c18..2349ee886 100644
--- a/tangelo/linq/target/backend.py
+++ b/tangelo/linq/target/backend.py
@@ -204,6 +204,8 @@ def simulate_circuit(self):
                 if available.
             initial_statevector (list/array) : A valid statevector in the format
                 supported by the target backend.
+            desired_meas_result (str): The binary string of the desired measurement.
+                Must have the same length as the number of MEASURE gates in source_circuit
             save_mid_circuit_meas (bool): Save mid-circuit measurement results to
                 self.mid_circuit_meas_freqs. All measurements will be saved to
                 self.all_frequencies, with keys of length (n_meas + n_qubits).
@@ -256,6 +258,8 @@ def simulate(self, source_circuit, return_statevector=False, initial_statevector
         """
         n_meas = source_circuit.counts.get("MEASURE", 0)
         n_cmeas = source_circuit.counts.get("CMEASURE", 0)
+        if n_cmeas > 0:
+            save_mid_circuit_meas = True
 
         if desired_meas_result is not None:
             if not isinstance(desired_meas_result, str) or (len(desired_meas_result) != n_meas and n_cmeas == 0):
diff --git a/tangelo/linq/target/target_cirq.py b/tangelo/linq/target/target_cirq.py
index 97b24ee92..7450b5603 100644
--- a/tangelo/linq/target/target_cirq.py
+++ b/tangelo/linq/target/target_cirq.py
@@ -76,7 +76,9 @@ def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, in
             raise NotImplementedError(f"{self.__class__.__name__} does not currently support measurement-controlled"
                                "gates with a noise model.")
 
-        # Only DensityMatrixSimulator handles noise well, can use Simulator, but it is slower
+        # Only DensityMatrixSimulator handles noise well, can use Simulator, but it is slower.
+        # When not saving mid-circuit measurements, the measurement operation can be converted to a noise channel
+        # such that the circuit can be simulated once and the density matrix sampled at the end.
         if self._noise_model or (source_circuit.is_mixed_state and not save_mid_circuit_meas):
             cirq_simulator = self.cirq.DensityMatrixSimulator(dtype=np.complex128)
         else:
diff --git a/tangelo/linq/target/target_stim.py b/tangelo/linq/target/target_stim.py
index 5d9509e0f..b7118ab8c 100644
--- a/tangelo/linq/target/target_stim.py
+++ b/tangelo/linq/target/target_stim.py
@@ -34,7 +34,8 @@ def __init__(self, n_shots=None, noise_model=None):
         super().__init__(n_shots=n_shots, noise_model=noise_model)
         self.stim = stim
 
-    def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None, desired_meas_result=None):
+    def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None,
+                         desired_meas_result=None, save_mid_circuit_meas=False):
         """Perform state preparation corresponding to the input circuit on the
         target backend, return the frequencies of the different observables, and
         either the statevector or None depending on if return_statevector is set to True.
@@ -47,6 +48,7 @@ def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, in
             return_statevector (bool): option to return the statevector
             initial_statevector (list/array) : Not currently implemented, will raise an error
             desired_meas_result (str) : Not currently implemented, will raise an error
+            save_mid_circuit_meas (bool): Not currently implemented, will raise an error
 
         Returns:
             dict: A dictionary mapping multi-qubit states to their corresponding
@@ -58,6 +60,8 @@ def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, in
             raise NotImplementedError("initial_statevector not yet implemented with stim ")
         if desired_meas_result is not None:
             raise NotImplementedError("desired_meas_result not yet implemented with stim ")
+        if save_mid_circuit_meas:
+            raise NotImplementedError("save_mid_circuit_meas not yet implemented with stim ")
         if "CMEASURE" in source_circuit.counts:
             raise NotImplementedError(f"{self.__class__.__name__} does not currently support CMEASURE operations.")
 
diff --git a/tangelo/linq/target/target_sympy.py b/tangelo/linq/target/target_sympy.py
index edbfe9bc7..81f91b259 100644
--- a/tangelo/linq/target/target_sympy.py
+++ b/tangelo/linq/target/target_sympy.py
@@ -25,7 +25,8 @@ class SympySimulator(Backend):
     def __init__(self, n_shots=None, noise_model=None):
         super().__init__(n_shots, noise_model)
 
-    def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None):
+    def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, initial_statevector=None,
+                         desired_meas_result=None, save_mid_circuit_meas=False):
         """This simulator manipulates symbolic expressions, i.e. gates can have
         unspecified parameters (strings interpreted as variables). As with the
         other simulators, it performs state preparation corresponding to the
@@ -40,6 +41,8 @@ def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, in
                 if available.
             initial_statevector (array/matrix or sympy.physics.quantum.Qubit): A
                 valid statevector in the format supported by the target backend.
+            desired_meas_result (str) : Not currently implemented, will raise an error
+            save_mid_circuit_meas (bool): Not currently implemented, will raise an error
 
         Returns:
             dict: A dictionary mapping multi-qubit states to their corresponding
@@ -55,6 +58,10 @@ def simulate_circuit(self, source_circuit: Circuit, return_statevector=False, in
 
         if "CMEASURE" in source_circuit.counts:
             raise NotImplementedError(f"{self.__class__.__name__} does not currently support CMEASURE operations.")
+        if save_mid_circuit_meas:
+            raise NotImplementedError(f"{self.__class__.__name__} does not currently support mid-circuit measurements.")
+        if desired_meas_result is not None:
+            raise NotImplementedError(f"{self.__class__.__name__} does not currently support mid-circuit measurements.")
 
         translated_circuit = translate_c(source_circuit, "sympy")
 
diff --git a/tangelo/linq/tests/test_simulator.py b/tangelo/linq/tests/test_simulator.py
index 00565036c..b32f842de 100644
--- a/tangelo/linq/tests/test_simulator.py
+++ b/tangelo/linq/tests/test_simulator.py
@@ -26,7 +26,7 @@
 from openfermion import load_operator, get_sparse_operator
 
 from tangelo.toolboxes.operators import QubitOperator
-from tangelo.linq import Gate, Circuit, get_backend, ClassicalControl
+from tangelo.linq import Gate, Circuit, get_backend, ClassicalControl, generate_applied_gates
 from tangelo.linq.translator import translate_circuit as translate_c
 from tangelo.linq.gate import PARAMETERIZED_GATES
 from tangelo.linq.target.backend import Backend, get_expectation_value_from_frequencies_oneterm
@@ -702,6 +702,7 @@ def cfunc(measurement):
             assert_freq_dict_almost_equal(f, {"11": 1}, 1.e-7)
             assert_freq_dict_almost_equal(circuit.success_probabilities, {"001": 0.125}, 1.e-7)
             self.assertTrue(applied_circuit == Circuit(circuit.applied_gates))
+            self.assertTrue(applied_circuit == Circuit(generate_applied_gates(circuit, desired_meas_result="001")))
 
             # Test that the correct ordering of probabilities is obtained when running without specifying the
             # desired measurement result.
@@ -734,7 +735,7 @@ def __init__(self, n_bits: int):
                     self.energies: List[float] = [0.]
                     self.n_runs: int = 0
 
-                def return_circuit(self, measurement) -> List[Gate]:
+                def return_gates(self, measurement) -> List[Gate]:
                     self.measurements[self.n_runs] += measurement
                     self.energies[self.n_runs] += int(measurement)/2**self.bitplace