diff --git a/hamlib/qiskit/hamlib_simulation_benchmark.py b/hamlib/qiskit/hamlib_simulation_benchmark.py
index 3b418133..8f4d06e8 100644
--- a/hamlib/qiskit/hamlib_simulation_benchmark.py
+++ b/hamlib/qiskit/hamlib_simulation_benchmark.py
@@ -26,6 +26,7 @@
 
 import hamlib_simulation_kernel
 from hamlib_simulation_kernel import HamiltonianSimulation, kernel_draw, get_valid_qubits
+from hamlib_simulation_kernel import initial_state, create_circuit   # would like to remove these
 from hamlib_utils import create_full_filenames, construct_dataset_name
 from hamiltonian_simulation_exact import HamiltonianSimulationExact, HamiltonianSimulation_Noiseless
 
@@ -83,7 +84,10 @@ def generate_pattern(starting_bit):
 ############### Result Data Analysis
 
 #def analyze_and_print_result(qc: QuantumCircuit, result, num_qubits: int,
-def analyze_and_print_result(qc, result, num_qubits: int,
+def analyze_and_print_result(
+            qc,
+            result,
+            num_qubits: int,
             type: str,
             num_shots: int,
             hamiltonian: str,
@@ -112,8 +116,6 @@ def analyze_and_print_result(qc, result, num_qubits: int,
 
     hamiltonian = hamiltonian.strip().lower()
 
-    from hamlib_simulation_kernel import initial_state, create_circuit
-
     # calculate correct distribution on the fly
     
     # for method 1, compute expected dist using ideal quantum simulation of the circuit provided
@@ -133,23 +135,34 @@ def analyze_and_print_result(qc, result, num_qubits: int,
             print(f"... begin exact computation for id={type} ...")
             
         ts = time.time()
-        
+
+        # DEVNOTE: ideally, we can remove these next two lines by performing this code in the run() loop        
         # create quantum circuit with initial state
-        qc_initial = initial_state(n_spins=num_qubits, initial_state=init_state)
+        qc_initial = initial_state(n_spins=num_qubits, init_state=init_state)
         
         # get Hamiltonian operator by creating entire circuit (DEVNOTE: need to not require whole circuit)
         _, ham_op, _ = create_circuit(n_spins=num_qubits, init_state=init_state)
         
-        # compute the exact evolution
+        # compute the expected  distribution after exact evolution
         correct_dist = HamiltonianSimulationExact(qc_initial, n_spins=num_qubits,
                 hamiltonian_op=ham_op,
                 time=1.0)
                 
         if verbose:
-            print(f"... exact computation time = {round((time.time() - ts), 3)} sec") 
+            print(f"... exact computation time = {round((time.time() - ts), 3)} sec")
             
+    # for method 3, compute expected distribution from the initial state
     elif method == 3:
-        correct_dist = key_from_initial_state(num_qubits, num_shots, init_state, random_pauli_flag)
+    
+        # check simple distribution if not inserting random Paulis
+        if not random_pauli_flag:
+            correct_dist = key_from_initial_state(num_qubits, num_shots, init_state, random_pauli_flag)
+            
+        # if using random paulis, distribution is based on all ones in initial state
+        # random_pauli_flag is (for now) set to always return bitstring of 1's
+        else:
+            correct_dist = {"1" * num_qubits: num_shots}
+
         
     else:
         raise ValueError("Method is not 1 or 2 or 3, or hamiltonian is not valid.")
@@ -313,11 +326,14 @@ def execution_handler(qc, result, num_qubits, type, num_shots):
             ts = time.time()
             
             # create the HamLibSimulation kernel and the associated Hamiltonian operator
-            qc, ham_op = HamiltonianSimulation(num_qubits,
+            qc, ham_op = HamiltonianSimulation(
+                    num_qubits,
                     hamiltonian=hamiltonian, 
                     K=k, t=t,
                     init_state=init_state,
-                    method = method)
+                    method = method,
+                    random_pauli_flag=random_pauli_flag
+                    )
                     
             metrics.store_metric(num_qubits, circuit_id, 'create_time', time.time() - ts)
 
diff --git a/hamlib/qiskit/hamlib_simulation_kernel.py b/hamlib/qiskit/hamlib_simulation_kernel.py
index 610c6669..858357cb 100644
--- a/hamlib/qiskit/hamlib_simulation_kernel.py
+++ b/hamlib/qiskit/hamlib_simulation_kernel.py
@@ -38,7 +38,13 @@
 INV_ = None
 
 
-from hamlib_utils import process_hamiltonian_file, needs_normalization, normalize_data_format, parse_hamiltonian_to_sparsepauliop, determine_qubit_count
+from hamlib_utils import (
+    process_hamiltonian_file,
+    needs_normalization,
+    normalize_data_format,
+    parse_hamiltonian_to_sparsepauliop,
+    determine_qubit_count,
+)
 
 def set_default_parameter_values(filename):
     """
@@ -215,7 +221,14 @@ def create_trotter_steps(num_trotter_steps, evo, operator, circuit):
     circuit.barrier()
     return circuit
 
-def create_circuit(n_spins: int, time: float = 1, num_trotter_steps: int = 5, method = 1, init_state=None):
+def create_circuit(
+    n_spins: int,
+    time: float = 1,
+    num_trotter_steps: int = 5,
+    method: int = 1,
+    init_state: str = None,
+    random_pauli_flag: bool = False,
+):
     """
     Create a quantum circuit based on the Hamiltonian data from an HDF5 file.
 
@@ -255,9 +268,6 @@ def create_circuit(n_spins: int, time: float = 1, num_trotter_steps: int = 5, me
         if verbose:
             print(f"... Evolution operator = {ham_op}")
 
-        ham_op = ham_op  # Use the SparsePauliOp object directly
-        # print (ham_op)
-
         # Build the evolution gate
         # label = "e\u2071\u1D34\u1D57"    # superscripted, but doesn't look good
         evo_label = "e^-iHt"
@@ -267,7 +277,7 @@ def create_circuit(n_spins: int, time: float = 1, num_trotter_steps: int = 5, me
         circuit = QuantumCircuit(ham_op.num_qubits)
         circuit_without_initial_state = QuantumCircuit(ham_op.num_qubits)
         
-        # first insert the initial_state
+        # first create and append the initial_state
         # init_state = "checkerboard"
         i_state = initial_state(num_qubits, init_state)
         circuit.append(i_state, range(ham_op.num_qubits))
@@ -283,40 +293,54 @@ def create_circuit(n_spins: int, time: float = 1, num_trotter_steps: int = 5, me
         # Append K Trotter steps of inverse, if method 3
         inv = None
         if method == 3: 
-            inv = evo.inverse()
-            inv.name = "e^iHt"
-            circuit = create_trotter_steps(num_trotter_steps, inv, ham_op, circuit)
-            if n_spins <= 6:
-                INV_ = inv
-            
-        circuit.measure_all()
+    
+            # if not adding random Paulis, just create simple inverse Trotter steps
+            if not random_pauli_flag:
+                inv = evo.inverse()
+                inv.name = "e^iHt"
+                circuit = create_trotter_steps(num_trotter_steps, inv, ham_op, circuit)
+                if n_spins <= 6:
+                    INV_ = inv
+     
+            # if adding Paulis, do that here, with code from pyGSTi
+            else:
+                from pygsti_mirror import convert_to_mirror_circuit
+                circuit, _ = convert_to_mirror_circuit(circuit_without_initial_state)
+                
+        # convert_to_mirror_circuit adds its own measurement gates    
+        if not random_pauli_flag:
+            circuit.measure_all()
+        
         return circuit, ham_op, evo
+        
     else:
         # print(f"Dataset not available for n_spins = {n_spins}.")
         return None, None, None
 
 
-############### Circuit Definition
+############### Initial Circuit Definition
 
-def initial_state(n_spins: int, initial_state: str = "checker") -> QuantumCircuit:
+def initial_state(n_spins: int, init_state: str = "checker") -> QuantumCircuit:
     """
     Initialize the quantum state.
     
     Args:
         n_spins (int): Number of spins (qubits).
-        initial_state (str): The chosen initial state. By default applies the checkerboard state, but can also be set to "ghz", the GHZ state.
+        init_state (str): The chosen initial state. By default applies the checkerboard state, but can also be set to "ghz", the GHZ state.
 
     Returns:
         QuantumCircuit: The initialized quantum circuit.
     """
     qc = QuantumCircuit(n_spins)
 
-    if initial_state.strip().lower() == "checkerboard" or initial_state.strip().lower() == "neele":
+    init_state = init_state.strip().lower()
+    
+    if init_state == "checkerboard" or init_state == "neele":
         # Checkerboard state, or "Neele" state
         qc.name = "Neele"
         for k in range(0, n_spins, 2):
             qc.x([k])
-    elif initial_state.strip().lower() == "ghz":
+    elif init_state.strip().lower() == "ghz":
         # GHZ state: 1/sqrt(2) (|00...> + |11...>)
         qc.name = "GHZ"
         qc.h(0)
@@ -325,12 +349,16 @@ def initial_state(n_spins: int, initial_state: str = "checker") -> QuantumCircui
 
     return qc
 
+############### Hamiltonian Circuit Definition
 
-def HamiltonianSimulation(n_spins: int,
+def HamiltonianSimulation(
+            n_spins: int,
             hamiltonian: str, 
             K: int = 5, t: float = 1.0,
-            init_state=None,
-            method: int = 1) -> QuantumCircuit:
+            init_state = None,
+            method: int = 1,
+            random_pauli_flag = False
+        ) -> QuantumCircuit:
     """
     Construct a Qiskit circuit for Hamiltonian simulation.
 
@@ -340,6 +368,7 @@ def HamiltonianSimulation(n_spins: int,
         K (int): The Trotterization order.
         t (float): Duration of simulation.
         method (int): Type of comparison for fidelity
+        random_pauli_flag (bool): Insert random Pauli gates if method 3
 
     Returns:
         QuantumCircuit: The constructed Qiskit circuit.
@@ -351,13 +380,18 @@ def HamiltonianSimulation(n_spins: int,
     qr = QuantumRegister(n_spins)
     cr = ClassicalRegister(n_spins)
     qc = QuantumCircuit(qr, cr, name=f"hamsim-{num_qubits}-{secret_int}")
-    
-    # size of one Trotter step
-    tau = t / K
 
     hamiltonian = hamiltonian.strip().lower()
     
-    qc, ham_op, evo = create_circuit(n_spins = n_spins, method = method, init_state=init_state)
+    # create the quantum circuit for this Hamiltonian, along with the operator and trotter evolution circuit
+    qc, ham_op, evo = create_circuit(
+        n_spins=n_spins,
+        time=t,
+        method=method,
+        init_state=init_state,
+        num_trotter_steps=K,
+        random_pauli_flag=random_pauli_flag,
+        )
 
     # Save smaller circuit example for display
     global QC_, HAM_, EVO_, INV_
diff --git a/hamlib/qiskit/pygsti_mirror.py b/hamlib/qiskit/pygsti_mirror.py
new file mode 100644
index 00000000..43ee9ded
--- /dev/null
+++ b/hamlib/qiskit/pygsti_mirror.py
@@ -0,0 +1,347 @@
+import pygsti
+import tqdm
+import numpy as np
+import scipy as sp
+from qiskit import transpile, QuantumCircuit
+
+# Below code created by Timothy Proctor, July 16 2024, using code developed by various members of Sandia's QPL.
+
+def sample_mirror_circuits(c, num_mcs=100, randomized_state_preparation=True, rand_state=None):
+    if rand_state is None:
+        rand_state = np.random.RandomState()
+        
+    mirror_circuits = []
+    mirror_circuits_bitstrings = []
+    for j in range(num_mcs):
+        mc, bs = central_pauli_mirror_circuit(c, randomized_state_preparation=randomized_state_preparation,
+                                              rand_state=rand_state)
+        mirror_circuits.append(mc)
+        mirror_circuits_bitstrings.append(bs)
+
+    return mirror_circuits, mirror_circuits_bitstrings
+
+def sample_reference_circuits(qubits, num_ref=100, randomized_state_preparation=True, rand_state=None):
+    if rand_state is None:
+        rand_state = np.random.RandomState()
+              
+    ref_circuits = []
+    ref_circuits_bitstrings = []
+    for j in range(num_ref):
+        empty_circuit = pygsti.circuits.Circuit('', line_labels=qubits)
+        ref_c, bs = central_pauli_mirror_circuit(empty_circuit, 
+                                                 randomized_state_preparation=randomized_state_preparation,
+                                                 rand_state=rand_state)
+        ref_circuits.append(ref_c)
+        ref_circuits_bitstrings.append(bs)
+
+    return ref_circuits, ref_circuits_bitstrings
+
+def central_pauli_mirror_circuit(circ, randomized_state_preparation=True, rand_state=None):
+    if rand_state is None:
+        rand_state = np.random.RandomState()
+
+    qubits = circ.line_labels
+    if randomized_state_preparation:
+        prep_circ = pygsti.circuits.Circuit([haar_random_u3_layer(qubits, rand_state)], line_labels=qubits)
+        circuit_to_mirror = prep_circ + circ
+    else: 
+        circuit_to_mirror = circ
+
+    n = circuit_to_mirror.width
+    d = circuit_to_mirror.depth
+
+    # SONNY'S EDIT TO TIM CODE: MAKE THE CENTRAL PAULI (for now) NOT RANDOM
+    # Note that we could rework the code to pass in a desired bitstring here later, rather than randomly generating one here.
+    #central_pauli = 2 * rand_state.randint(0, 2, 2*n)
+    central_pauli = 2 * np.ones(2*n, dtype = np.int8 )  
+    central_pauli_layer = pauli_vector_to_u3_layer(central_pauli, qubits)
+    q = central_pauli.copy()
+
+    quasi_inverse_circ = pygsti.circuits.Circuit(line_labels=circ.line_labels, editable=True)
+
+    for i in range(d):
+        
+        layer = circuit_to_mirror.layer_label(d - i - 1).components
+        quasi_inverse_layer = [gate_inverse(gate_label) for gate_label in layer]
+        
+        # Update the u3 gates.
+        if len(layer) == 0 or layer[0].name == 'Gu3':
+            # update_u3_parameters now twirls/adds label for empty qubits, so don't prepad for speed
+            #padded_layer = pad_layer(quasi_inverse_layer, qubits)
+            #quasi_inverse_layer = update_u3_parameters(padded_layer, q, q, qubits)
+            quasi_inverse_layer = update_u3_parameters(quasi_inverse_layer, q, q, qubits)
+            
+        # Update q based on the CNOTs in the layer.
+        else:
+            for g in layer:
+                if g.name == 'Gcnot':
+                    (control, target) = g.qubits
+                    q[qubits.index(control)] = (q[qubits.index(control)] + q[qubits.index(target)]) % 4
+                    q[n + qubits.index(target)] = (q[n + qubits.index(control)] + q[n + qubits.index(target)]) % 4
+                else:
+                    raise ValueError("Circuit can only contain Gcnot and Gu3 gates in separate layers!")
+
+        quasi_inverse_circ.insert_layer_inplace(quasi_inverse_layer, i)
+
+    mc = circuit_to_mirror + pygsti.circuits.Circuit([central_pauli_layer], line_labels=circ.line_labels) + quasi_inverse_circ
+    mc.done_editing()
+
+    bs = ''.join([str(b // 2) for b in q[n:]])
+
+    return mc, bs
+
+def pauli_vector_to_u3_layer(p, qubits):
+    
+    n = len(qubits)
+    layer = []
+    for i, q in enumerate(qubits):
+
+        if p[i] == 0 and p[i+n] == 0:  # I
+            theta = 0.0
+            phi = 0.0
+            lamb = 0.0
+        if p[i] == 2 and p[i+n] == 0:  # Z
+            theta = 0.0
+            phi = np.pi / 2
+            lamb = np.pi / 2
+        if p[i] == 0 and p[i+n] == 2:  # X
+            theta = np.pi
+            phi = 0.0
+            lamb = np.pi
+        if p[i] == 2 and p[i+n] == 2:  # Y
+            theta = np.pi
+            phi = np.pi / 2
+            lamb = np.pi / 2
+
+        layer.append(pygsti.baseobjs.Label('Gu3', q, args=(theta, phi, lamb)))
+
+    return pygsti.baseobjs.Label(layer)
+
+def haar_random_u3_layer(qubits, rand_state=None):
+    
+    return pygsti.baseobjs.Label([haar_random_u3(q, rand_state) for q in qubits])
+
+def haar_random_u3(q, rand_state=None):
+    if rand_state is None:
+        rand_state = np.random.RandomState()
+
+    a, b = 2 * np.pi * rand_state.rand(2)
+    theta = mod_2pi(2 * np.arcsin(np.sqrt(rand_state.rand(1)))[0])
+    phi = mod_2pi(a - b + np.pi)
+    lamb = mod_2pi(-1 * (a + b + np.pi))
+    return pygsti.baseobjs.Label('Gu3', q, args=(theta, phi, lamb))
+
+def mod_2pi(theta):
+    while (theta > np.pi or theta <= -1 * np.pi):
+        if theta > np.pi:
+            theta = theta - 2 * np.pi
+        elif theta <= -1 * np.pi:
+            theta = theta + 2 * np.pi
+    return theta
+
+def inverse_u3(args):
+    theta_inv = mod_2pi(-float(args[0]))
+    phi_inv = mod_2pi(-float(args[2]))
+    lambda_inv = mod_2pi(-float(args[1]))
+    return (theta_inv, phi_inv, lambda_inv)
+
+def gate_inverse(label):
+    if label.name == 'Gcnot':
+        return label
+    elif label.name == 'Gu3':
+        return pygsti.baseobjs.label.LabelTupWithArgs.init('Gu3', label.qubits, args=inverse_u3(label.args))
+
+def pad_layer(layer, qubits):
+
+    padded_layer = list(layer)
+    used_qubits = []
+    for g in layer:
+        for q in g.qubits:
+            used_qubits.append(q)
+
+    for q in qubits:
+        if q not in used_qubits:
+            padded_layer.append(pygsti.baseobjs.label.LabelTupWithArgs.init('Gu3', (q,), args=(0.0, 0.0, 0.0)))
+
+    return padded_layer
+
+
+def update_u3_parameters(layer, p, q, qubits):
+    """
+    Takes a layer containing u3 gates, and finds a new layer containing
+    u3 gates that implements p * layer * q (p followed by layer followed by
+    q), where p and q are vectors  describing layers of paulis.
+
+    """
+    used_qubits = []
+
+    new_layer = []
+    n = len(qubits)
+
+    for g in layer:
+        assert(g.name == 'Gu3')
+        (theta, phi, lamb) = (float(g.args[0]), float(g.args[1]), float(g.args[2]))
+        qubit_index = qubits.index(g.qubits[0])
+        if p[qubit_index] == 2:   # Z gate preceeding the layer
+            lamb = lamb + np.pi
+        if q[qubit_index] == 2:   # Z gate following the layer
+            phi = phi + np.pi
+        if p[n + qubit_index] == 2:  # X gate preceeding the layer
+            theta = theta - np.pi
+            phi = phi
+            lamb = -lamb - np.pi
+        if q[n + qubit_index] == 2:  # X gate following the layer
+            theta = theta - np.pi
+            phi = -phi - np.pi
+            lamb = lamb
+
+        new_args = (mod_2pi(theta), mod_2pi(phi), mod_2pi(lamb))
+        new_label = pygsti.baseobjs.label.LabelTupWithArgs.init('Gu3', g.qubits[0], args=new_args)
+        new_layer.append(new_label)
+        used_qubits.append(g.qubits[0])
+    
+#     for qubit_index, qubit in enumerate(qubits):
+#         if qubit in used_qubits:
+#             continue
+
+#         # Insert twirled idle on unpadded qubit
+#         (theta, phi, lamb) = (0.0, 0.0, 0.0)
+#         if p[qubit_index] == 2:   # Z gate preceeding the layer
+#             lamb = lamb + np.pi
+#         if q[qubit_index] == 2:   # Z gate following the layer
+#             phi = phi + np.pi
+#         if p[n + qubit_index] == 2:  # X gate preceeding the layer
+#             theta = theta - np.pi
+#             phi = phi
+#             lamb = -lamb - np.pi
+#         if q[n + qubit_index] == 2:  # X gate following the layer
+#             theta = theta - np.pi
+#             phi = -phi - np.pi
+#             lamb = lamb
+        
+#         new_args = (mod_2pi(theta), mod_2pi(phi), mod_2pi(lamb))
+#         new_label = pygsti.baseobjs.label.LabelTupWithArgs.init('Gu3', qubit, args=new_args)
+#         new_layer.append(new_label)
+#         used_qubits.append(qubit)
+#
+#    assert(set(used_qubits) == set(qubits))
+
+    return new_layer
+
+def hamming_distance_counts(data_dict, circ, idealout):
+    nQ = len(circ.line_labels)  # number of qubits
+    assert(nQ == len(idealout))
+    total = np.sum(list(data_dict.values()))
+    hamming_distance_counts = np.zeros(nQ + 1, float)
+    for outcome, counts in data_dict.items():
+        hamming_distance_counts[pygsti.tools.rbtools.hamming_distance(outcome, idealout)] += counts
+    return hamming_distance_counts
+
+def adjusted_success_probability(hamming_distance_counts):
+    if np.sum(hamming_distance_counts) == 0.: 
+        return 0.
+    else:
+        hamming_distance_pdf = np.array(hamming_distance_counts) / np.sum(hamming_distance_counts)
+        adjSP = np.sum([(-1 / 2)**n * hamming_distance_pdf[n] for n in range(len(hamming_distance_pdf))])
+        return adjSP
+
+def effective_polarization(hamming_distance_counts):
+    n = len(hamming_distance_counts) - 1 
+    asp = adjusted_success_probability(hamming_distance_counts)
+    
+    return (4**n * asp - 1)/(4**n - 1)
+
+def polarization_to_fidelity(p, n): 
+    return 1 - (4**n - 1)*(1 - p)/4**n
+
+def fidelity_to_polarization(f, n):
+    return 1 - (4**n)*(1 - f)/(4**n - 1)
+
+def predicted_process_fidelity(mirror_circuit_effective_pols, reference_effective_pols, n):
+    a = np.mean(mirror_circuit_effective_pols)
+    c = np.mean(reference_effective_pols)
+    if c <= 0.:
+        return np.nan  # raise ValueError("Reference effective polarization zero or smaller! Cannot estimate the process fidelity")
+    elif a <= 0:
+        return 0.
+    else:
+        return polarization_to_fidelity(np.sqrt(a / c), n)
+    
+def predicted_process_fidelity_with_spam_error(mirror_circuit_effective_pols, reference_effective_pols, n):
+    a = np.mean(mirror_circuit_effective_pols)
+    c = np.mean(reference_effective_pols)
+    if c <= 0.:
+        return np.nan  # raise ValueError("Reference effective polarization zero or smaller! Cannot estimate the process fidelity")
+    elif a <= 0:
+        return 0.
+    else:
+        return polarization_to_fidelity(np.sqrt(c * a), n)
+
+# Below code created by Sonny Rappaport based off of the above code provided by Tim Proctor. 
+
+def pygsti_string(qc):
+    """
+    Takes a qiskit circuit (that has been transpiled to the u3 and cx basis), and returns a pygsti-readable string. 
+    """
+
+    pygstr = ""
+
+    for gate in qc.data:
+
+        gate_str = "[" 
+
+        gate_name = gate[0].name
+        gate_qubits = [] 
+
+        for qubit in gate[1]:
+
+            # qubits look something like Qubit(QuantumRegister(4,'q1'),1)
+            gate_qubits.append(qubit._index)
+        
+        gate_parameters = gate[0].params
+
+        if gate_name == "cx":
+            gate_str = gate_str + "Gcnot"
+        elif gate_name == "u3":
+            gate_str = gate_str + "Gu3"
+        else:
+            # print(f"Passed {gate_name}")
+
+            continue
+
+        for parameter in gate_parameters:
+
+            gate_str = gate_str + f";{parameter}" 
+
+        for qubit in gate_qubits:
+
+            gate_str = gate_str + f":Q{qubit}"
+
+        gate_str = gate_str + "]"
+
+        pygstr = pygstr + gate_str
+
+    pygstr = pygstr + "@(" + f"{''.join(f'Q{i},' for i in range(qc.num_qubits))[:-1]}" + ")"
+
+    return pygstr.strip()
+
+def convert_to_mirror_circuit(qc):
+    """
+    Takes an arbitrary quantum circuit, transpiles it to CX and u3, and eventually returns a mirror circuit version of the quantum circuit. 
+
+    The mirror circuit will have a random initial state in addition to random paulis. 
+    """
+
+    qc = transpile(qc, basis_gates=["cx","u3"])
+
+    pygstr = pygsti_string(qc)
+
+    pygsti_circuit = pygsti.circuits.Circuit(pygstr)
+
+    mcs, bss = sample_mirror_circuits(pygsti_circuit, num_mcs=1)
+    print(bss)
+
+    qiskit_circuit = QuantumCircuit.from_qasm_str(mcs[0].convert_to_openqasm())
+    
+    # this circuit is made out of u3 and cx gates by pygsti default
+    return qiskit_circuit , bss[0]