diff --git a/tangelo/algorithms/variational/vqe_solver.py b/tangelo/algorithms/variational/vqe_solver.py
index 9d703b3e3..81f19c5cd 100644
--- a/tangelo/algorithms/variational/vqe_solver.py
+++ b/tangelo/algorithms/variational/vqe_solver.py
@@ -21,12 +21,11 @@
 
 from enum import Enum
 import numpy as np
-from openfermion.ops.operators.qubit_operator import QubitOperator
 
 from tangelo.helpers.utils import HiddenPrints
 from tangelo.linq import get_backend, Circuit
 from tangelo.linq.helpers.circuits.measurement_basis import measurement_basis_gates
-from tangelo.toolboxes.operators import count_qubits, FermionOperator
+from tangelo.toolboxes.operators import count_qubits, FermionOperator, 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
 from tangelo.toolboxes.post_processing.bootstrapping import get_resampled_frequencies
diff --git a/tangelo/linq/__init__.py b/tangelo/linq/__init__.py
index 65304e50b..c826cad75 100644
--- a/tangelo/linq/__init__.py
+++ b/tangelo/linq/__init__.py
@@ -18,3 +18,4 @@
 from .simulator import get_backend
 from .target.backend import get_expectation_value_from_frequencies_oneterm
 from .target import backend_info, Backend
+from .noisy_simulation import NoiseModel
diff --git a/tangelo/linq/target/backend.py b/tangelo/linq/target/backend.py
index 4281201a7..c5754ec67 100644
--- a/tangelo/linq/target/backend.py
+++ b/tangelo/linq/target/backend.py
@@ -37,10 +37,10 @@
 import numpy as np
 from scipy import stats
 from bitarray import bitarray
-from openfermion.ops import QubitOperator
 
 from tangelo.linq import Gate, Circuit
 from tangelo.linq.helpers.circuits.measurement_basis import measurement_basis_gates
+from tangelo.toolboxes.operators import QubitOperator
 
 
 def get_expectation_value_from_frequencies_oneterm(term, frequencies):
@@ -48,8 +48,7 @@ def get_expectation_value_from_frequencies_oneterm(term, frequencies):
     the result of a state-preparation.
 
     Args:
-        term (openfermion-style QubitOperator object): a qubit operator, with
-            only a single term.
+        term (QubitOperator): a single-term qubit operator.
         frequencies (dict): histogram of frequencies of measurements (assumed
             to be in lsq-first format).
 
@@ -82,8 +81,7 @@ def get_variance_from_frequencies_oneterm(term, frequencies):
     """Return the variance of the expectation value of a single-term qubit-operator, given
     the result of a state-preparation.
     Args:
-        term (openfermion-style QubitOperator object): a qubit operator, with
-            only a single term.
+        term (QubitOperator): a single-term qubit operator.
         frequencies (dict): histogram of frequencies of measurements (assumed
             to be in lsq-first format).
     Returns:
@@ -315,8 +313,7 @@ def get_expectation_value(self, qubit_operator, state_prep_circuit, initial_stat
         actual QPU.
 
         Args:
-            qubit_operator (openfermion-style QubitOperator class): a qubit
-                operator.
+            qubit_operator (QubitOperator): the qubit operator.
             state_prep_circuit (Circuit): an abstract circuit used for state preparation.
             initial_statevector (array): The initial statevector for the simulation
             desired_meas_result (str): The mid-circuit measurement results to select for.
@@ -376,8 +373,7 @@ def get_variance(self, qubit_operator, state_prep_circuit, initial_statevector=N
         actual QPU.
 
         Args:
-            qubit_operator (openfermion-style QubitOperator class): a qubit
-                operator.
+            qubit_operator (QubitOperator): the qubit operator.
             state_prep_circuit (Circuit): an abstract circuit used for state preparation.
             initial_statevector (list/array) : A valid statevector in the format
                 supported by the target backend.
@@ -431,8 +427,7 @@ def get_standard_error(self, qubit_operator, state_prep_circuit, initial_stateve
         actual QPU.
 
         Args:
-            qubit_operator (openfermion-style QubitOperator class): a qubit
-                operator.
+            qubit_operator (QubitOperator): the qubit operator.
             state_prep_circuit (Circuit): an abstract circuit used for state preparation.
             initial_statevector (list/array): A valid statevector in the format
                 supported by the target backend.
@@ -452,7 +447,7 @@ def _get_expectation_value_from_statevector(self, qubit_operator, state_prep_cir
         this function directly, please call "get_expectation_value" instead.
 
         Args:
-            qubit_operator (openfermion-style QubitOperator class): a qubit operator.
+            qubit_operator (QubitOperator): the qubit operator.
             state_prep_circuit (Circuit): an abstract circuit used for state preparation (only pure states).
             initial_statevector (array): The initial state of the system
 
@@ -503,7 +498,7 @@ def _get_expectation_value_from_frequencies(self, qubit_operator, state_prep_cir
         using the frequencies of observable states.
 
         Args:
-            qubit_operator (openfermion-style QubitOperator class): a qubitoperator.
+            qubit_operator (QubitOperator): the qubit operator.
             state_prep_circuit (Circuit): an abstract circuit used for state preparation.
             initial_statevector (array): The initial state of the system
             desired_meas_result (str): The mid-circuit measurement results to select for.
@@ -551,7 +546,7 @@ def _get_variance_from_frequencies(self, qubit_operator, state_prep_circuit, ini
         using the frequencies of observable states.
 
         Args:
-            qubit_operator (openfermion-style QubitOperator class): a qubit operator.
+            qubit_operator (QubitOperator): the qubit operator.
             state_prep_circuit (Circuit): an abstract circuit used for state preparation.
             initial_statevector (list/array) : A valid statevector in the format
                 supported by the target backend.
@@ -599,8 +594,7 @@ def get_expectation_value_from_frequencies_oneterm(term, frequencies):
         the result of a state-preparation.
 
         Args:
-            term (openfermion-style QubitOperator object): a qubit operator, with
-                only a single term.
+            term (QubitOperator): a single-term qubit operator
             frequencies (dict): histogram of frequencies of measurements (assumed
                 to be in lsq-first format).
 
@@ -617,8 +611,7 @@ def get_variance_from_frequencies_oneterm(term, frequencies):
         the result of a state-preparation.
 
         Args:
-            term (openfermion-style QubitOperator object): a qubit operator, with
-                only a single term.
+            term (QubitOperator): a single-term qubit operator.
             frequencies (dict): histogram of frequencies of measurements (assumed
                 to be in lsq-first format).
 
diff --git a/tangelo/linq/tests/test_simulator.py b/tangelo/linq/tests/test_simulator.py
index 0a85ad26d..8917c5363 100644
--- a/tangelo/linq/tests/test_simulator.py
+++ b/tangelo/linq/tests/test_simulator.py
@@ -22,15 +22,14 @@
 import time
 
 import numpy as np
-from openfermion.ops import QubitOperator
 from openfermion import load_operator, get_sparse_operator
 
+from tangelo.toolboxes.operators import QubitOperator
 from tangelo.linq import Gate, Circuit, get_backend
 from tangelo.linq.translator import translate_circuit as translate_c
 from tangelo.linq.gate import PARAMETERIZED_GATES
-from tangelo.helpers.utils import installed_simulator, installed_sv_simulator, installed_backends
 from tangelo.linq.target.backend import Backend, get_expectation_value_from_frequencies_oneterm
-from tangelo.helpers.utils import assert_freq_dict_almost_equal
+from tangelo.helpers.utils import installed_simulator, installed_sv_simulator, installed_backends, assert_freq_dict_almost_equal
 
 path_data = os.path.dirname(os.path.abspath(__file__)) + '/data'
 
@@ -425,7 +424,7 @@ def test_get_exp_value_mixed_state_desired_measurement_with_shots(self):
         """ Get expectation value of mixed state by post-selecting on desired measurement."""
         qubit_operator = QubitOperator("X0 X1") + QubitOperator("Y0 Y1") + QubitOperator("Z0 Z1") + QubitOperator("X0 Y1", 1j)
 
-        ham = get_sparse_operator(qubit_operator).toarray()
+        ham = get_sparse_operator(qubit_operator.to_openfermion()).toarray()
         exact_sv = np.array([0.+0.j, 0.+0.j, 0.87758256+0.j, -0.47942554+0.j])
         exact_exp = np.vdot(exact_sv, ham @ exact_sv)
 
diff --git a/tangelo/linq/tests/test_simulator_noisy.py b/tangelo/linq/tests/test_simulator_noisy.py
index 2d79408cc..19be463c9 100644
--- a/tangelo/linq/tests/test_simulator_noisy.py
+++ b/tangelo/linq/tests/test_simulator_noisy.py
@@ -19,11 +19,11 @@
 import unittest
 
 import numpy as np
-from openfermion.ops import QubitOperator
 
 from tangelo.linq import Gate, Circuit, get_backend, backend_info
 from tangelo.linq.noisy_simulation import NoiseModel, get_qiskit_noise_dict
 from tangelo.helpers.utils import default_simulator, installed_backends, assert_freq_dict_almost_equal
+from tangelo.toolboxes.operators import QubitOperator
 
 # Noisy simulation: circuits, noise models, references
 cn1 = Circuit([Gate('X', target=0)])
diff --git a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
index 0d1631fad..2f5e4dd8d 100644
--- a/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
+++ b/tangelo/toolboxes/ansatz_generator/ansatz_utils.py
@@ -22,12 +22,11 @@
 
 import numpy as np
 from openfermion.ops import FermionOperator as ofFermionOperator
-from openfermion.ops import InteractionOperator as ofInteractionOperator
 from openfermion.ops import QubitOperator as ofQubitOperator
 
 from tangelo.linq import Circuit, Gate
 from tangelo.toolboxes.operators import FermionOperator, QubitOperator
-from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping, get_fermion_operator
+from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
 
 
 def pauli_op_to_gate(index, op, inverse=False):
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_ansatz_util.py b/tangelo/toolboxes/ansatz_generator/tests/test_ansatz_util.py
index 03818f694..9e4080566 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_ansatz_util.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_ansatz_util.py
@@ -44,7 +44,7 @@ def test_trotterize_fermion_input(self):
                                                          n_electrons=mol_H4_sto3g.n_active_electrons,
                                                          up_then_down=True)
 
-            ham_mat = get_sparse_operator(qubit_hamiltonian).toarray()
+            ham_mat = get_sparse_operator(qubit_hamiltonian.to_openfermion()).toarray()
             evolve_exact = expm(-1j * time * ham_mat) @ refwave
 
             options = {"up_then_down": True,
@@ -75,7 +75,7 @@ def test_trotterize_qubit_input(self):
                                                          n_electrons=mol_H4_sto3g.n_active_electrons,
                                                          up_then_down=True)
 
-            ham_mat = get_sparse_operator(qubit_hamiltonian).toarray()
+            ham_mat = get_sparse_operator(qubit_hamiltonian.to_openfermion()).toarray()
             evolve_exact = expm(-1j * time * ham_mat) @ refwave
 
             tcircuit, phase = trotterize(qubit_hamiltonian, trotter_order=1, n_trotter_steps=1, time=time,
@@ -103,7 +103,7 @@ def test_trotterize_different_order_and_steps(self):
                                                      n_electrons=mol_H4_sto3g.n_active_electrons,
                                                      up_then_down=True)
 
-        ham_mat = get_sparse_operator(qubit_hamiltonian).toarray()
+        ham_mat = get_sparse_operator(qubit_hamiltonian.to_openfermion()).toarray()
         evolve_exact = expm(-1j * time * ham_mat) @ refwave
 
         for trotter_order, n_trotter_steps in [(1, 1), (2, 1), (1, 2), (4, 1), (6, 1)]:
@@ -144,7 +144,7 @@ def test_trotterize_fermionic_input_different_times(self):
             total_fermion_operator += fermion_operators[i]
             qubit_hamiltonian = fermion_to_qubit_mapping(fermion_operator=fermion_operators[i],
                                                          mapping=mapping)
-            ham_mat = get_sparse_operator(qubit_hamiltonian, n_qubits=4).toarray()
+            ham_mat = get_sparse_operator(qubit_hamiltonian.to_openfermion(), n_qubits=4).toarray()
             evolve_exact = expm(-1j * time[next(iter(fermion_operators[i].terms))] * ham_mat) @ evolve_exact
 
         # Apply trotter-suzuki steps using different times for each term
@@ -171,7 +171,7 @@ def test_trotterize_qubit_input_different_times(self):
         # Exactly evolve for each time step
         evolve_exact = refwave
         for i in range(3):
-            ham_mat = get_sparse_operator(qubit_operator_list[i], n_qubits=4).toarray()
+            ham_mat = get_sparse_operator(qubit_operator_list[i].to_openfermion(), n_qubits=4).toarray()
             evolve_exact = expm(-1j * time[next(iter(qubit_operator_list[i].terms))] * ham_mat) @ evolve_exact
 
         # Apply trotter-suzuki with different times for each qubit operator term
@@ -219,7 +219,7 @@ def test_controlled_time_evolution_by_phase_estimation(self):
         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).toarray()
+        ham_mat = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, wavefunction = eigh(ham_mat)
 
         # Append four qubits in the zero state to eigenvector 9
@@ -275,9 +275,9 @@ def test_derangement_circuit_by_estimating_pauli_string(self):
         qu_op = QubitOperator("X0 Y1", 0.125) + QubitOperator("Y1 Y2", 0.125) + QubitOperator("Z2 Z3", 0.125)
         pa = QubitOperator("X0 X1 X2 X3", 1)
 
-        ham_mat = get_sparse_operator(qu_op).toarray()
+        ham_mat = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, wavefunction = eigh(ham_mat)
-        pamat = get_sparse_operator(pa).toarray()
+        pamat = get_sparse_operator(pa.to_openfermion()).toarray()
 
         mixed_wave = np.sqrt(3)/2*wavefunction[:, -1] + 1/2*wavefunction[:, 0]
         mixed_wave_3 = np.kron(np.kron(mixed_wave, mixed_wave), mixed_wave)
diff --git a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
index a54703d54..eff018a61 100644
--- a/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/tests/test_vsqs.py
@@ -18,11 +18,11 @@
 import numpy as np
 from openfermion import load_operator
 
+from tangelo.linq import get_backend
 from tangelo.molecule_library import mol_H2_sto3g
 from tangelo.toolboxes.operators import QubitOperator
 from tangelo.toolboxes.qubit_mappings import jordan_wigner, symmetry_conserving_bravyi_kitaev
 from tangelo.toolboxes.ansatz_generator import VSQS
-from tangelo.linq import get_backend
 from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
 
 # For openfermion.load_operator function.
@@ -71,17 +71,21 @@ def test_vsqs_H2(self):
         self.assertAlmostEqual(energy, -1.1372701255155757, delta=1e-6)
 
     def test_vsqs_H4_doublecation(self):
-        """Verify closed-shell VSQS functionalities for H4 2+ by using saved qubit hamiltonian and initial hamiltonian"""
+        """ Verify closed-shell VSQS functionalities for H4 2+ by using saved qubit and initial hamiltonians """
 
         var_params = [-2.53957674, 0.72683888, 1.08799500, 0.49836183,
                       -0.23020698, 0.93278630, 0.50591026, 0.50486903]
 
-        # Build qubit hamiltonian for energy evaluation
-        qubit_hamiltonian = load_operator("mol_H4_doublecation_minao_qubitham_jw_b.data", data_directory=pwd_this_test+"/data", plain_text=True)
-        initial_hamiltonian = load_operator("mol_H4_doublecation_minao_init_qubitham_jw_b.data", data_directory=pwd_this_test+"/data", plain_text=True)
-        reference_state = get_reference_circuit(8, 2, "jw", up_then_down=True, spin=0)
+        # Build qubit hamiltonian for energy evaluation, loading them with Openfermion.
+        qubit_ofhamiltonian = load_operator("mol_H4_doublecation_minao_qubitham_jw_b.data", data_directory=pwd_this_test+"/data", plain_text=True)
+        initial_ofhamiltonian = load_operator("mol_H4_doublecation_minao_init_qubitham_jw_b.data", data_directory=pwd_this_test+"/data", plain_text=True)
+
+        # Load into Tangelo operators
+        qubit_hamiltonian = QubitOperator.from_openfermion(qubit_ofhamiltonian)
+        initial_hamiltonian = QubitOperator.from_openfermion(initial_ofhamiltonian)
 
         # Build circuit
+        reference_state = get_reference_circuit(8, 2, "jw", up_then_down=True, spin=0)
         vsqs_ansatz = VSQS(qubit_hamiltonian=qubit_hamiltonian, h_init=initial_hamiltonian, reference_state=reference_state,
                            intervals=5, time=5, trotter_order=2)
         vsqs_ansatz.build_circuit()
diff --git a/tangelo/toolboxes/ansatz_generator/uccgd.py b/tangelo/toolboxes/ansatz_generator/uccgd.py
index c3ef43a05..4996ea35f 100644
--- a/tangelo/toolboxes/ansatz_generator/uccgd.py
+++ b/tangelo/toolboxes/ansatz_generator/uccgd.py
@@ -24,10 +24,10 @@
 
 import numpy as np
 from openfermion import hermitian_conjugated
-from openfermion import FermionOperator as ofFermionOperator
 from itertools import combinations_with_replacement, product
 
 from tangelo.linq import Circuit
+from tangelo.toolboxes.operators import FermionOperator
 from tangelo.toolboxes.ansatz_generator.ansatz import Ansatz
 from tangelo.toolboxes.ansatz_generator.ansatz_utils import get_exponentiated_qubit_operator_circuit
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
@@ -180,7 +180,7 @@ def _get_qubit_operator(self):
         Returns:
             QubitOperator: qubit-encoded elements of the UCCGD
         """
-        fermion_op = ofFermionOperator()
+        fermion_op = FermionOperator()
         n_mos = self.n_spinorbitals // 2
         p = -1
         for indices in combinations_with_replacement(range(n_mos), 4):
@@ -188,8 +188,8 @@ def _get_qubit_operator(self):
                 u, w, v, t = indices
                 p = p + 1
                 for sig, tau in product(range(2), repeat=2):
-                    c_op = ofFermionOperator(((2*t+sig, 1), (2*v+tau, 1), (2*w+tau, 0), (2*u+sig, 0)), self.var_params[p])
-                    c_op += ofFermionOperator(((2*v+sig, 1), (2*t+tau, 1), (2*u+tau, 0), (2*w+sig, 0)), self.var_params[p])
+                    c_op = FermionOperator(((2*t+sig, 1), (2*v+tau, 1), (2*w+tau, 0), (2*u+sig, 0)), self.var_params[p])
+                    c_op += FermionOperator(((2*v+sig, 1), (2*t+tau, 1), (2*u+tau, 0), (2*w+sig, 0)), self.var_params[p])
                     fermion_op += c_op - hermitian_conjugated(c_op)
 
         qubit_op = fermion_to_qubit_mapping(fermion_operator=fermion_op,
@@ -198,6 +198,9 @@ def _get_qubit_operator(self):
                                             n_electrons=self.n_electrons,
                                             up_then_down=self.up_then_down)
 
+        print(type(fermion_op))
+        print(type(qubit_op))
+
         # Cast all coefs to floats (rotations angles are real)
         for key in qubit_op.terms:
             qubit_op.terms[key] = float(qubit_op.terms[key].imag)
diff --git a/tangelo/toolboxes/ansatz_generator/vsqs.py b/tangelo/toolboxes/ansatz_generator/vsqs.py
index 80c1f8597..d1bf6259c 100644
--- a/tangelo/toolboxes/ansatz_generator/vsqs.py
+++ b/tangelo/toolboxes/ansatz_generator/vsqs.py
@@ -16,12 +16,11 @@
 Adiabatic State Preparation (ASP) inspired ansatz as described in https://arxiv.org/abs/2003.09913."""
 
 import numpy as np
-from openfermion import QubitOperator as ofQubitOperator
 
 from .ansatz import Ansatz
 from .ansatz_utils import get_exponentiated_qubit_operator_circuit
-from tangelo.toolboxes.operators import FermionOperator
 from tangelo.linq import Circuit
+from tangelo.toolboxes.operators import FermionOperator, QubitOperator
 from tangelo.toolboxes.qubit_mappings.mapping_transform import get_qubit_number, fermion_to_qubit_mapping
 from tangelo.toolboxes.qubit_mappings.statevector_mapping import get_reference_circuit
 
@@ -68,7 +67,7 @@ def __init__(self, molecule=None, mapping="jw", up_then_down=False, intervals=2,
 
         if molecule is None:
             self.qubit_hamiltonian = qubit_hamiltonian
-            if not isinstance(h_init, ofQubitOperator):
+            if not isinstance(h_init, QubitOperator):
                 raise ValueError("When providing a qubit hamiltonian, an initial qubit Hamiltonian must also be provided")
             self.h_init = h_init
             if not isinstance(reference_state, Circuit):
@@ -102,7 +101,7 @@ def qu_op_to_list(qu_op):
             self.stride = 2
             self.n_h_nav = 0
         else:
-            if isinstance(self.h_nav, ofQubitOperator):
+            if isinstance(self.h_nav, QubitOperator):
                 self.stride = 3
                 self.h_nav_list = qu_op_to_list(self.h_nav)
                 self.n_h_nav = len(self.h_nav_list)
diff --git a/tangelo/toolboxes/circuits/tests/test_diagonal_coulomb.py b/tangelo/toolboxes/circuits/tests/test_diagonal_coulomb.py
index 6e37a89c8..bcf29ec78 100644
--- a/tangelo/toolboxes/circuits/tests/test_diagonal_coulomb.py
+++ b/tangelo/toolboxes/circuits/tests/test_diagonal_coulomb.py
@@ -14,7 +14,7 @@
 
 import unittest
 
-from openfermion import get_sparse_operator, linalg
+from openfermion import linalg, get_sparse_operator
 
 from tangelo.linq import get_backend, Circuit
 from tangelo.molecule_library import mol_H4_sto3g
@@ -24,7 +24,7 @@
 sim = get_backend(target="cirq")
 
 
-class diagonal_coulomb_Test(unittest.TestCase):
+class DiagonalCoulombTest(unittest.TestCase):
 
     def test_orbital_rotations(self):
         """Test calculating energy expectation value of H4 hamiltonian by decomposing into diagional terms"""
diff --git a/tangelo/toolboxes/circuits/tests/test_discrete_clock.py b/tangelo/toolboxes/circuits/tests/test_discrete_clock.py
index df8059058..5ae0ef431 100644
--- a/tangelo/toolboxes/circuits/tests/test_discrete_clock.py
+++ b/tangelo/toolboxes/circuits/tests/test_discrete_clock.py
@@ -15,17 +15,16 @@
 import unittest
 import math
 
-from openfermion import get_sparse_operator
 import numpy as np
 from scipy.linalg import expm
+from openfermion import get_sparse_operator
 
-from tangelo.linq import get_backend, Circuit
+from tangelo.molecule_library import mol_H2_sto3g
 from tangelo.helpers.utils import installed_backends
+from tangelo.linq import get_backend, Circuit
 from tangelo.linq.helpers.circuits.statevector import StateVector
-from tangelo.toolboxes.operators.operators import QubitOperator
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
-from tangelo.toolboxes.ansatz_generator.ansatz_utils import get_qft_circuit, trotterize
-from tangelo.molecule_library import mol_H2_sto3g
+from tangelo.toolboxes.ansatz_generator.ansatz_utils import trotterize
 from tangelo.toolboxes.circuits.discrete_clock import get_discrete_clock_circuit
 from tangelo.toolboxes.circuits.grid_circuits import get_psquared_circuit, get_xsquared_circuit
 
@@ -45,7 +44,7 @@ def test_time_independant_hamiltonian(self):
         qu_op = fermion_to_qubit_mapping(mol_H2_sto3g.fermionic_hamiltonian, "scbk", mol_H2_sto3g.n_active_sos, mol_H2_sto3g.n_active_electrons,
                                          True, 0)
 
-        ham = get_sparse_operator(qu_op).toarray()
+        ham = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, vecs = np.linalg.eigh(ham)
         vec = (vecs[:, 0] + vecs[:, 1])/np.sqrt(2)
 
@@ -62,8 +61,8 @@ def trotter_func(t0, time, n_trotter_steps, control):
             sv_circuit = sv.initializing_circuit()
 
             for k in [2, 3]:
-                taylor_circuit = get_discrete_clock_circuit(trotter_func=trotter_func, trotter_kwargs={}, time=time, mp_order=k, n_state_qus=2,
-                                                            n_time_steps=4)
+                taylor_circuit = get_discrete_clock_circuit(trotter_func=trotter_func, trotter_kwargs={},
+                                                            time=time, mp_order=k, n_state_qus=2, n_time_steps=4)
                 _, v = sim.simulate(sv_circuit + taylor_circuit, return_statevector=True)
                 n_ancilla = 2 + math.ceil(np.log2(k+2))
                 len_ancilla = 2**n_ancilla
@@ -108,9 +107,9 @@ def trotter_func(t0, time, n_trotter_steps, control, dx, qubit_list):
 
             for k in [2, 3]:
                 qubit_list = list(reversed(range(n_qubits))) if statevector_order == "lsq_first" else list((range(n_qubits)))
-                taylor_circuit = get_discrete_clock_circuit(trotter_func=trotter_func, trotter_kwargs={"dx": dx, "qubit_list": qubit_list}, time=time,
-                                                            mp_order=k, n_state_qus=6,
-                                                            n_time_steps=2)
+                taylor_circuit = get_discrete_clock_circuit(trotter_func=trotter_func,
+                                                            trotter_kwargs={"dx": dx, "qubit_list": qubit_list}, time=time,
+                                                            mp_order=k, n_state_qus=6, n_time_steps=2)
                 _, v = sim.simulate(sv_circuit + taylor_circuit, return_statevector=True)
                 n_ancilla = 1 + math.ceil(np.log2(k+2))
                 len_ancilla = 2**n_ancilla
diff --git a/tangelo/toolboxes/circuits/tests/test_grid.py b/tangelo/toolboxes/circuits/tests/test_grid.py
index ae2df84df..4f96297e0 100644
--- a/tangelo/toolboxes/circuits/tests/test_grid.py
+++ b/tangelo/toolboxes/circuits/tests/test_grid.py
@@ -13,11 +13,8 @@
 # limitations under the License.
 
 import unittest
-import math
 
-from openfermion import get_sparse_operator
 import numpy as np
-from scipy.linalg import expm
 
 from tangelo.linq import get_backend
 from tangelo.linq.helpers.circuits.statevector import StateVector
diff --git a/tangelo/toolboxes/circuits/tests/test_lcu.py b/tangelo/toolboxes/circuits/tests/test_lcu.py
index 0af757d74..692251134 100644
--- a/tangelo/toolboxes/circuits/tests/test_lcu.py
+++ b/tangelo/toolboxes/circuits/tests/test_lcu.py
@@ -15,9 +15,9 @@
 import unittest
 import math
 
-from openfermion import get_sparse_operator
 import numpy as np
 from scipy.linalg import expm
+from openfermion import get_sparse_operator
 
 from tangelo.linq import get_backend
 from tangelo.helpers.utils import installed_backends
@@ -36,7 +36,7 @@
 sim_cirq = get_backend("cirq")
 
 
-class lcu_Test(unittest.TestCase):
+class LCUTest(unittest.TestCase):
 
     def test_get_truncated_taylor_series(self):
         """Test time-evolution of truncated Taylor series for different orders and times"""
@@ -44,7 +44,7 @@ def test_get_truncated_taylor_series(self):
         qu_op = fermion_to_qubit_mapping(mol_H2_sto3g.fermionic_hamiltonian, "scbk", mol_H2_sto3g.n_active_sos, mol_H2_sto3g.n_active_electrons,
                                          True, 0)
         n_qubits_qu_op = math.ceil(math.log2(len(qu_op.terms)))
-        ham = get_sparse_operator(qu_op).toarray()
+        ham = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, vecs = np.linalg.eigh(ham)
         vec = (vecs[:, 0] + vecs[:, 1])/np.sqrt(2)
 
@@ -65,21 +65,23 @@ def test_get_truncated_taylor_series(self):
                 v = v.reshape([4, len_ancilla])[:, 0] if statevector_order == "lsq_first" else v.reshape([len_ancilla, 4])[0, :]
                 self.assertAlmostEqual(1, np.abs(v.conj().dot(exact)), delta=3.e-1**k)
 
-        # Raise ValueError if Taylor series order is less than 1 or greater than 4 or imaginary coefficients in qubit operator
+        # Raise ValueError if Taylor series order is less than 1 or greater than 4
+        # or imaginary coefficients in qubit operator
         self.assertRaises(ValueError, get_truncated_taylor_series, qu_op, 0, time)
         self.assertRaises(ValueError, get_truncated_taylor_series, qu_op * 1j, 2, time)
 
     def test_get_oaa_lcu_circuit(self):
-        """Test time-evolution of truncated Taylor series for order k = 3 passing explicitly calculated qubit operator exponential"""
+        """Test time-evolution of truncated Taylor series for order k = 3 passing explicitly calculated
+        qubit operator exponential"""
 
-        qu_op = fermion_to_qubit_mapping(mol_H2_sto3g.fermionic_hamiltonian, "scbk", mol_H2_sto3g.n_active_sos, mol_H2_sto3g.n_active_electrons,
-                                         True, 0)
+        qu_op = fermion_to_qubit_mapping(mol_H2_sto3g.fermionic_hamiltonian, "scbk",
+                                         mol_H2_sto3g.n_active_sos, mol_H2_sto3g.n_active_electrons, True, 0)
         time = 0.5
         # Generate explicit qubit operator exponential
         exp_qu_op = 1 + -1j*qu_op*time + (-1j*qu_op*time)**2/2 + (-1j*qu_op*time)**3/6
         exp_qu_op.compress()
         n_qubits_qu_op = math.ceil(math.log2(len(exp_qu_op.terms)))
-        ham = get_sparse_operator(qu_op).toarray()
+        ham = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, vecs = np.linalg.eigh(ham)
         vec = (vecs[:, 0] + vecs[:, 1])/np.sqrt(2)
 
@@ -109,7 +111,7 @@ def test_controlled_time_evolution_by_phase_estimation_for_get_truncated_taylor_
         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).toarray()
+        ham_mat = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, wavefunction = np.linalg.eigh(ham_mat)
 
         # Kronecker product 13 qubits in the zero state to eigenvector 9 to account for ancilla qubits
@@ -143,7 +145,7 @@ def test_controlled_time_evolution_by_phase_estimation_for_get_oaa_lcu_circuit(s
         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).toarray()
+        ham_mat = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, wavefunction = np.linalg.eigh(ham_mat)
 
         # break time into 6 parts so 1-norm is less than 2. i.e. can use Oblivious Amplitude Amplification
diff --git a/tangelo/toolboxes/circuits/tests/test_mp.py b/tangelo/toolboxes/circuits/tests/test_mp.py
index 7d8c43ffd..ec2c09739 100644
--- a/tangelo/toolboxes/circuits/tests/test_mp.py
+++ b/tangelo/toolboxes/circuits/tests/test_mp.py
@@ -13,11 +13,10 @@
 # limitations under the License.
 
 import unittest
-import math
 
-from openfermion import get_sparse_operator
 import numpy as np
 from scipy.linalg import expm
+from openfermion import get_sparse_operator
 
 from tangelo.linq import get_backend
 from tangelo.helpers.utils import installed_backends
@@ -41,10 +40,11 @@ class MultiProductTest(unittest.TestCase):
     def test_time_evolution(self):
         """Test time-evolution of multi-product circuit for different orders"""
 
-        qu_op = fermion_to_qubit_mapping(mol_H2_sto3g.fermionic_hamiltonian, "scbk", mol_H2_sto3g.n_active_sos, mol_H2_sto3g.n_active_electrons,
+        qu_op = fermion_to_qubit_mapping(mol_H2_sto3g.fermionic_hamiltonian, "scbk",
+                                         mol_H2_sto3g.n_active_sos, mol_H2_sto3g.n_active_electrons,
                                          True, 0)
 
-        ham = get_sparse_operator(qu_op).toarray()
+        ham = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, vecs = np.linalg.eigh(ham)
         vec = (vecs[:, 0] + vecs[:, 1])/np.sqrt(2)
 
@@ -81,7 +81,7 @@ def test_controlled_time_evolution_by_phase_estimation(self):
         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).toarray()
+        ham_mat = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, wavefunction = np.linalg.eigh(ham_mat)
 
         # Kronecker product 13 qubits in the zero state to eigenvector 9 to account for ancilla qubits
@@ -93,7 +93,8 @@ def test_controlled_time_evolution_by_phase_estimation(self):
 
         pe_circuit = get_qft_circuit(qubit_list)
         for i, qubit in enumerate(qubit_list):
-            pe_circuit += get_multi_product_circuit(operator=qu_op, n_state_qus=4, order=5, time=-(2*np.pi)*2**i, control=qubit)
+            pe_circuit += get_multi_product_circuit(operator=qu_op, n_state_qus=4, order=5,
+                                                    time=-(2*np.pi)*2**i, control=qubit)
         pe_circuit += get_qft_circuit(qubit_list, inverse=True)
 
         freqs, _ = sim_cirq.simulate(pe_circuit, initial_statevector=wave_9)
diff --git a/tangelo/toolboxes/circuits/tests/test_qsp.py b/tangelo/toolboxes/circuits/tests/test_qsp.py
index 73c76b650..bb69c9b06 100644
--- a/tangelo/toolboxes/circuits/tests/test_qsp.py
+++ b/tangelo/toolboxes/circuits/tests/test_qsp.py
@@ -14,14 +14,14 @@
 
 import unittest
 
-from openfermion import get_sparse_operator
 import numpy as np
 from scipy.linalg import expm
+from openfermion import get_sparse_operator
 
 from tangelo.linq import get_backend, backend_info
 from tangelo.helpers.utils import installed_backends
 from tangelo.linq.helpers.circuits.statevector import StateVector
-from tangelo.toolboxes.operators.operators import QubitOperator
+from tangelo.toolboxes.operators import QubitOperator
 from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping
 from tangelo.toolboxes.ansatz_generator.ansatz_utils import get_qft_circuit
 from tangelo.molecule_library import mol_H2_sto3g
@@ -36,16 +36,16 @@
 sim_cirq = get_backend("cirq")
 
 
-class lcu_Test(unittest.TestCase):
+class QSPTest(unittest.TestCase):
 
     def test_get_qsp_circuit(self):
         """Test QSP time-evolution"""
 
-        qu_op = fermion_to_qubit_mapping(mol_H2_sto3g.fermionic_hamiltonian, "scbk", mol_H2_sto3g.n_active_sos, mol_H2_sto3g.n_active_electrons,
-                                         True, 0)
+        qu_op = fermion_to_qubit_mapping(mol_H2_sto3g.fermionic_hamiltonian, "scbk",
+                                         mol_H2_sto3g.n_active_sos, mol_H2_sto3g.n_active_electrons, True, 0)
         # need to ensure eigenvalues are between -1 and 1
         qu_op /= 1.2
-        ham = get_sparse_operator(qu_op).toarray()
+        ham = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, vecs = np.linalg.eigh(ham)
         vec = (vecs[:, 0] + vecs[:, 2])/np.sqrt(2)
 
@@ -76,7 +76,7 @@ def test_controlled_time_evolution_by_phase_estimation_for_get_qsp_circuit(self)
         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).toarray()
+        ham_mat = get_sparse_operator(qu_op.to_openfermion()).toarray()
         _, wavefunction = np.linalg.eigh(ham_mat)
 
         # Kronecker product 13 qubits in the zero state to eigenvector 9 to account for ancilla qubits
diff --git a/tangelo/toolboxes/measurements/tests/test_qubit_terms_grouping.py b/tangelo/toolboxes/measurements/tests/test_qubit_terms_grouping.py
index cabe4f1e2..132791af6 100644
--- a/tangelo/toolboxes/measurements/tests/test_qubit_terms_grouping.py
+++ b/tangelo/toolboxes/measurements/tests/test_qubit_terms_grouping.py
@@ -16,10 +16,10 @@
 import os
 
 from openfermion import load_operator
-from openfermion.ops import QubitOperator
 
 from tangelo.linq import translator, get_backend, Circuit
 from tangelo.helpers import measurement_basis_gates
+from tangelo.toolboxes.operators import QubitOperator
 from tangelo.toolboxes.measurements import group_qwc, exp_value_from_measurement_bases, \
     check_bases_commute_qwc, map_measurements_qwc
 
diff --git a/tangelo/toolboxes/operators/operators.py b/tangelo/toolboxes/operators/operators.py
index 58059af1a..6490a7a74 100644
--- a/tangelo/toolboxes/operators/operators.py
+++ b/tangelo/toolboxes/operators/operators.py
@@ -168,6 +168,20 @@ class QubitOperator(of.QubitOperator):
     replaced by our own implementation.
     """
 
+    @classmethod
+    def from_openfermion(cls, of_qop):
+        """ Enable instantiation of a QubitOperator from an openfermion QubitOperator object.
+
+        Args:
+            of_qop (openfermion QubitOperator): an existing qubit operator defined with Openfermion
+
+        Returns:
+            corresponding QubitOperator object.
+        """
+        qop = cls()
+        qop.terms = of_qop.terms.copy()
+        return qop
+
     def frobenius_norm_compression(self, epsilon, n_qubits):
         """Reduces the number of operator terms based on its Frobenius norm
         and a user-defined threshold, epsilon. The eigenspectrum of the
@@ -230,6 +244,12 @@ def qubit_indices(self):
 
         return qubit_indices
 
+    def to_openfermion(self):
+        """Converts Tangelo QubitOperator to openfermion"""
+        qu_op = of.QubitOperator()
+        qu_op.terms = self.terms.copy()
+        return qu_op
+
 
 class QubitHamiltonian(QubitOperator):
     """QubitHamiltonian objects are essentially openfermion.QubitOperator
@@ -243,8 +263,7 @@ class QubitHamiltonian(QubitOperator):
         coefficient (complex): Coefficient for this term.
         mapping (string): Mapping procedure for fermionic to qubit encoding
             (ex: "JW", "BK", etc.).
-        up_then_down (bool): Whether or not spin ordering is all up then
-            all down.
+        up_then_down (bool): Whether spin ordering is all up then all down.
 
     Properties:
         n_terms (int): Number of terms in this qubit Hamiltonian.
diff --git a/tangelo/toolboxes/qubit_mappings/statevector_mapping.py b/tangelo/toolboxes/qubit_mappings/statevector_mapping.py
index f8f6adea1..2f74d5025 100644
--- a/tangelo/toolboxes/qubit_mappings/statevector_mapping.py
+++ b/tangelo/toolboxes/qubit_mappings/statevector_mapping.py
@@ -20,7 +20,7 @@
 import warnings
 
 import numpy as np
-from openfermion import QubitOperator
+from openfermion.ops.operators import QubitOperator
 from openfermion.transforms import bravyi_kitaev_code
 from openfermion.transforms.opconversions.bravyi_kitaev_tree import _transform_ladder_operator, FenwickTree