diff --git a/src/qutip_qip/circuit/circuit.py b/src/qutip_qip/circuit/circuit.py
index 61cca8586..9ef5af7e0 100644
--- a/src/qutip_qip/circuit/circuit.py
+++ b/src/qutip_qip/circuit/circuit.py
@@ -490,30 +490,19 @@ def run(
         final_state : Qobj
                 output state of the circuit run.
         """
-
         if state.isket:
-            sim = CircuitSimulator(
-                self,
-                state,
-                cbits,
-                U_list,
-                measure_results,
-                "state_vector_simulator",
-                precompute_unitary,
-            )
+            mode = "state_vector_simulator"
         elif state.isoper:
-            sim = CircuitSimulator(
-                self,
-                state,
-                cbits,
-                U_list,
-                measure_results,
-                "density_matrix_simulator",
-                precompute_unitary,
-            )
+            mode = "density_matrix_simulator"
         else:
             raise TypeError("State is not a ket or a density matrix.")
-        return sim.run(state, cbits).get_final_states(0)
+        sim = CircuitSimulator(
+            self,
+            U_list,
+            mode,
+            precompute_unitary,
+        )
+        return sim.run(state, cbits, measure_results).get_final_states(0)
 
     def run_statistics(
         self, state, U_list=None, cbits=None, precompute_unitary=False
@@ -530,11 +519,6 @@ def run_statistics(
                 initialization of the classical bits.
         U_list: list of Qobj, optional
             list of predefined unitaries corresponding to circuit.
-        measure_results : tuple of ints, optional
-            optional specification of each measurement result to enable
-            post-selection. If specified, the measurement results are
-            set to the tuple of bits (sequentially) instead of being
-            chosen at random.
         precompute_unitary: Boolean, optional
             Specify if computation is done by pre-computing and aggregating
             gate unitaries. Possibly a faster method in the case of
@@ -546,27 +530,18 @@ def run_statistics(
             Return a CircuitResult object containing
             output states and and their probabilities.
         """
-
         if state.isket:
-            sim = CircuitSimulator(
-                self,
-                state,
-                cbits,
-                U_list,
-                mode="state_vector_simulator",
-                precompute_unitary=precompute_unitary,
-            )
+            mode = "state_vector_simulator"
         elif state.isoper:
-            sim = CircuitSimulator(
-                self,
-                state,
-                cbits,
-                U_list,
-                mode="density_matrix_simulator",
-                precompute_unitary=precompute_unitary,
-            )
+            mode = "density_matrix_simulator"
         else:
             raise TypeError("State is not a ket or a density matrix.")
+        sim = CircuitSimulator(
+            self,
+            U_list,
+            mode,
+            precompute_unitary,
+        )
         return sim.run_statistics(state, cbits)
 
     def resolve_gates(self, basis=["CNOT", "RX", "RY", "RZ"]):
diff --git a/src/qutip_qip/circuit/circuitsimulator.py b/src/qutip_qip/circuit/circuitsimulator.py
index d0b2daacb..d756fcd8a 100644
--- a/src/qutip_qip/circuit/circuitsimulator.py
+++ b/src/qutip_qip/circuit/circuitsimulator.py
@@ -11,6 +11,7 @@
     gate_sequence_product,
 )
 from qutip import basis, ket2dm, Qobj, tensor
+import warnings
 
 
 __all__ = ["CircuitSimulator", "CircuitResult"]
@@ -253,12 +254,12 @@ class CircuitSimulator:
     def __init__(
         self,
         qc,
-        state=None,
-        cbits=None,
         U_list=None,
-        measure_results=None,
         mode="state_vector_simulator",
         precompute_unitary=False,
+        state=None,
+        cbits=None,
+        measure_results=None,
     ):
         """
         Simulate state evolution for Quantum Circuits.
@@ -268,21 +269,9 @@ def __init__(
         qc : :class:`.QubitCircuit`
             Quantum Circuit to be simulated.
 
-        state: ket or oper
-            ket or density matrix
-
-        cbits: list of int, optional
-            initial value of classical bits
-
         U_list: list of Qobj, optional
             list of predefined unitaries corresponding to circuit.
 
-        measure_results : tuple of ints, optional
-            optional specification of each measurement result to enable
-            post-selection. If specified, the measurement results are
-            set to the tuple of bits (sequentially) instead of being
-            chosen at random.
-
         mode: string, optional
             Specify if input state (and therefore computation) is in
             state-vector mode or in density matrix mode.
@@ -321,7 +310,13 @@ def __init__(
         else:
             self._process_ops()
 
-        self.initialize(state, cbits, measure_results)
+        if any(p is not None for p in (state, cbits, measure_results)):
+            warnings.warn(
+                "Initializing the quantum state, cbits and measure_results "
+                "when initializing the simulator is deprecated. "
+                "The inputs are ignored. "
+                "They should, instead, be provided when running the simulation."
+            )
 
     def _process_ops(self):
         """
diff --git a/tests/test_circuit.py b/tests/test_circuit.py
index 713a2cc0c..157e581d2 100644
--- a/tests/test_circuit.py
+++ b/tests/test_circuit.py
@@ -598,6 +598,23 @@ def test_measurement_circuit(self):
                 else:
                     assert simulator.cbits[0] != simulator.cbits[1]
 
+    def test_circuit_with_selected_measurement_result(self):
+        qc = QubitCircuit(N=1, num_cbits=1)
+        qc.add_gate("SNOT", targets=0)
+        qc.add_measurement("M0", targets=0, classical_store=0)
+
+        # We reset the random seed so that
+        # if we don's select the measurement result,
+        # the two circuit should return the same value.
+        np.random.seed(0)
+        final_state = qc.run(qp.basis(2, 0), cbits=[0], measure_results=[0])
+        fid = pytest.approx(qp.fidelity(final_state, basis(2, 0)))
+        assert fid == 1.0
+        np.random.seed(0)
+        final_state = qc.run(qp.basis(2, 0), cbits=[0], measure_results=[1])
+        fid = pytest.approx(qp.fidelity(final_state, basis(2, 1)))
+        assert fid == 1.0
+
     def test_gate_product(self):
 
         filename = "qft.qasm"