diff --git a/tangelo/linq/circuit.py b/tangelo/linq/circuit.py
index e0368960c..d635d22ee 100644
--- a/tangelo/linq/circuit.py
+++ b/tangelo/linq/circuit.py
@@ -224,7 +224,7 @@ def trim_qubits(self):
         """Trim unnecessary qubits and update indices with the lowest values possible.
         """
         qubits_in_use = set().union(*self.get_entangled_indices())
-        mapping = {ind: i for i, ind in enumerate(qubits_in_use)}
+        mapping = {ind: i for i, ind in enumerate(sorted(list(qubits_in_use)))}
         for g in self._gates:
             g.target = [mapping[ind] for ind in g.target]
             if g.control:
diff --git a/tangelo/toolboxes/operators/tests/test_trim_trivial_qubits.py b/tangelo/toolboxes/operators/tests/test_trim_trivial_qubits.py
new file mode 100644
index 000000000..ed819f9a9
--- /dev/null
+++ b/tangelo/toolboxes/operators/tests/test_trim_trivial_qubits.py
@@ -0,0 +1,92 @@
+# 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.
+
+import os
+import unittest
+
+import numpy as np
+from openfermion.linalg import qubit_operator_sparse
+from openfermion.utils import load_operator
+
+from tangelo.linq import Gate, Circuit, get_backend
+from tangelo.toolboxes.operators import QubitOperator
+from tangelo.toolboxes.ansatz_generator.ansatz_utils import exp_pauliword_to_gates
+from tangelo.toolboxes.operators.trim_trivial_qubits import trim_trivial_qubits, trim_trivial_operator, trim_trivial_circuit, is_bitflip_gate
+
+
+pwd_this_test = os.path.dirname(os.path.abspath(__file__))
+
+qb_ham = load_operator("H4_JW_spinupfirst.data", data_directory=pwd_this_test+"/data", plain_text=True)
+
+# Generate reference and test circuits using single qcc generator
+ref_qcc_op = 0.2299941483397896 * 0.5 * QubitOperator("Y0 X1 X2 X3")
+qcc_op = 0.2299941483397896 * 0.5 * QubitOperator("Y1 X3 X5 X7")
+
+ref_mf_gates = [Gate("RX", 0, parameter=np.pi), Gate("X", 2)]
+
+mf_gates = [
+            Gate("RZ", 0, parameter=np.pi/2), Gate("RX", 0, parameter=3.14159),
+            Gate("RX", 1, parameter=np.pi), Gate("RZ", 2, parameter=np.pi),
+            Gate("X", 4), Gate("X", 5), Gate("Z", 6),
+            Gate("RZ", 6, parameter=np.pi), Gate("RX", 8, parameter=-3*np.pi),
+            Gate("X", 8), Gate("RZ", 9, parameter=np.pi), Gate("Z", 9)
+           ]
+
+ref_pauli_words_gates = sum((exp_pauliword_to_gates(pword, coef) for pword, coef in ref_qcc_op.terms.items()), start=[])
+pauli_words_gates = sum((exp_pauliword_to_gates(pword, coef) for pword, coef in qcc_op.terms.items()), start=[])
+
+ref_circ = Circuit(ref_mf_gates + ref_pauli_words_gates)
+circ = Circuit(mf_gates + pauli_words_gates)
+
+# Reference energy for H4 molecule with single QCC generator
+ref_value = -1.8039875664891176
+
+# Reference indices and states to be removed from system
+ref_trim_states = {0: 1, 2: 0, 4: 1, 6: 0, 8: 0, 9: 0}
+
+# Reference for which mf_gates are bitflip gates
+ref_bool = [False, True, True, False, True, True, False, False, True, True, False, False]
+
+sim = get_backend()
+
+
+class TrimTrivialQubits(unittest.TestCase):
+    def test_trim_trivial_operator(self):
+        """ Test if trimming operator returns the correct eigenvalue """
+
+        trimmed_operator = trim_trivial_operator(qb_ham, trim_states={key: ref_trim_states[key] for key in [0, 2, 4, 6]}, reindex=False)
+        self.assertAlmostEqual(np.min(np.linalg.eigvalsh(qubit_operator_sparse(trimmed_operator).todense())), ref_value, places=5)
+
+    def test_is_bitflip_gate(self):
+        """ Test if bitflip gate function correctly identifies bitflip gates """
+        self.assertEqual(ref_bool, [is_bitflip_gate(g) for g in mf_gates])
+
+    def test_trim_trivial_circuit(self):
+        """ Test if circuit trimming returns the correct circuit, states, and indices  """
+
+        trimmed_circuit, trim_states = trim_trivial_circuit(circ)
+        self.assertEqual(ref_circ._gates, trimmed_circuit._gates)
+        self.assertEqual(ref_trim_states, trim_states)
+
+    def test_trim_trivial_qubits(self):
+        """ Test if trim trivial qubit function produces correct and compatible circuits and operators """
+
+        trimmed_operator, trimmed_circuit = trim_trivial_qubits(qb_ham, circ)
+        self.assertAlmostEqual(np.min(np.linalg.eigvalsh(qubit_operator_sparse(trimmed_operator).todense())), ref_value, places=5)
+        self.assertEqual(ref_circ._gates, trimmed_circuit._gates)
+        self.assertAlmostEqual(sim.get_expectation_value(trimmed_operator, trimmed_circuit), ref_value, places=5)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tangelo/toolboxes/operators/trim_trivial_qubits.py b/tangelo/toolboxes/operators/trim_trivial_qubits.py
new file mode 100644
index 000000000..695da9ec0
--- /dev/null
+++ b/tangelo/toolboxes/operators/trim_trivial_qubits.py
@@ -0,0 +1,171 @@
+# 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.
+
+import numpy as np
+
+from tangelo.toolboxes.operators import QubitOperator, count_qubits
+from tangelo.linq import Circuit
+from tangelo.linq.helpers.circuits import pauli_string_to_of, pauli_of_to_string
+
+
+def trim_trivial_operator(qu_op, trim_states, n_qubits=None, reindex=True):
+    """
+    Calculate expectation values of a QubitOperator acting on qubits in a
+    trivial |0> or |1> state. Return a trimmed QubitOperator with updated coefficients
+
+    Args:
+        qu_op (QubitOperator): Operator to trim
+        trim_states (dict): Dictionary mapping qubit indices to states to trim, e.g. {1: 0, 3: 1}
+        n_qubits (int): Optional, number of qubits in full system
+        reindex (bool): Optional, if True, remaining qubits will be reindexed
+    Returns:
+        QubitOperator : trimmed QubitOperator with updated coefficients
+    """
+
+    qu_op_trim = QubitOperator()
+    n_qubits = count_qubits(qu_op) if n_qubits is None else n_qubits
+
+    # Calculate expectation values of trivial qubits, update coefficients
+    for op, coeff in qu_op.terms.items():
+        term = pauli_of_to_string(op, n_qubits)
+        c = np.ones(len(trim_states))
+        new_term = term
+        for i, qubit in enumerate(trim_states.keys()):
+            if term[qubit] in {'X', 'Y'}:
+                c[i] = 0
+                break
+            elif (term[qubit], trim_states[qubit]) == ('Z', 1):
+                c[i] = -1
+
+            new_term = new_term[:qubit - i] + new_term[qubit - i + 1:] if reindex else new_term[:qubit] + 'I' + new_term[qubit + 1:]
+
+        if 0 in c:
+            continue
+
+        qu_op_trim += np.prod(c) * coeff * QubitOperator(pauli_string_to_of(new_term))
+    return qu_op_trim
+
+
+def is_bitflip_gate(gate, atol=1e-5):
+    """
+    Check if a gate is a bitflip gate.
+
+    A gate is a bitflip gate if it satisfies one of the following conditions:
+    1. The gate name is either X, Y.
+    2. The gate name is RX or RY, and has a parameter that is an odd multiple of pi.
+
+    Args:
+        gate (Gate): The gate to check.
+        atol (float): Optional, the absolute tolerance for gate parameter
+
+    Returns:
+        bool: True if the gate is a single qubit bitflip gate, False otherwise.
+    """
+    if gate is None:
+        return False
+
+    if gate.name in {"X", "Y"}:
+        return True
+    elif gate.name in {"RX", "RY"}:
+        try:
+            parameter_float = float(gate.parameter)
+        except (TypeError, ValueError):
+            return False
+
+        # Check if parameter is close to an odd multiple of pi
+        return abs(parameter_float % (np.pi * 2) - np.pi) <= atol
+    else:
+        return False
+
+
+def trim_trivial_circuit(circuit):
+    """
+        Splits Circuit into entangled and unentangled components.
+        Returns entangled Circuit, and the indices and states of unentangled qubits
+
+        Args:
+            circuit (Circuit): circuit to be trimmed
+        Returns:
+            Circuit : Trimmed, entangled circuit
+            dict : dictionary mapping trimmed qubit indices to their states (0 or 1)
+
+    """
+    # Split circuit and get relevant indices
+    circs = circuit.split()
+    e_indices = circuit.get_entangled_indices()
+
+    # Find qubits with no gates applied to them, store qubit index and state |0>
+    trim_states = {}
+    for qubit_idx in set(range(circuit.width)) - set(circuit._qubit_indices):
+        trim_states[qubit_idx] = 0
+
+    circuit_new = Circuit()
+    # Go through circuit components, trim if trivial, otherwise append to new circuit
+    for i, circ in enumerate(circs):
+        if circ.width != 1 or circ.size not in (1, 2):
+            circuit_new += circ
+            continue
+
+        # Calculate state of single qubit clifford circuits, ideally this would be done with a clifford simulator
+        # for now only look at first two gate combinations typical of the QMF state in QCC methods
+        gate0, gate1 = circ._gates[:2] + [None] * (2 - circ.size)
+        gate_0_is_bitflip = is_bitflip_gate(gate0)
+        gate_1_is_bitflip = is_bitflip_gate(gate1)
+
+        if circ.size == 1:
+            if gate0.name in {"RZ", "Z"}:
+                qubit_idx = e_indices[i].pop()
+                trim_states[qubit_idx] = 0
+            elif gate0.name in {"X", "RX"} and gate_0_is_bitflip:
+                qubit_idx = e_indices[i].pop()
+                trim_states[qubit_idx] = 1
+            else:
+                circuit_new += circ
+        elif circ.size == 2:
+            if gate1.name in {"Z", "RZ"}:
+                if gate0.name in {"RZ", "Z"}:
+                    qubit_idx = e_indices[i].pop()
+                    trim_states[qubit_idx] = 0
+                else:
+                    circuit_new += circ
+            elif gate1.name in {"X", "RX"} and gate_1_is_bitflip:
+                if gate0.name in {"RX", "X"} and gate_0_is_bitflip:
+                    qubit_idx = e_indices[i].pop()
+                    trim_states[qubit_idx] = 0
+                elif gate0.name in {"Z", "RZ"}:
+                    qubit_idx = e_indices[i].pop()
+                    trim_states[qubit_idx] = 1
+                else:
+                    circuit_new += circ
+            else:
+                circuit_new += circ
+    return circuit_new, trim_states
+
+
+def trim_trivial_qubits(operator, circuit):
+    """
+        Trim circuit and operator based on expectation values calculated from
+        trivial components of the circuit.
+
+        Args:
+            operator (QubitOperator): Operator to trim
+            circuit (Circuit): circuit to be trimmed
+        Returns:
+            QubitOperator : Trimmed qubit operator
+            Circuit : Trimmed circuit
+    """
+    trimmed_circuit, trim_states = trim_trivial_circuit(circuit)
+    trimmed_operator = trim_trivial_operator(operator, trim_states, circuit.width, reindex=True)
+
+    return trimmed_operator, trimmed_circuit