diff --git a/tests/test_molecule.py b/tests/test_molecule.py
index d35d532..4cf4ecc 100644
--- a/tests/test_molecule.py
+++ b/tests/test_molecule.py
@@ -9,8 +9,8 @@
 from qibo.hamiltonians import SymbolicHamiltonian
 from qibo.symbols import X, Y, Z
 
-from qibochem.driver.molecule import Molecule
 from qibochem.driver import hamiltonian
+from qibochem.driver.molecule import Molecule
 from qibochem.measurement.expectation import expectation
 
 
@@ -42,47 +42,82 @@ def test_run_pyscf_molecule_xyz():
 
 
 def test_fermionic_hamiltonian():
-
     h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
     h2.run_pyscf()
     h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
-    
+
     # a^\dagger_0 a_0
-    ref_one_body_tensor = np.array([[-1.27785301,  0.,          0.,          0.        ],
-                                    [ 0.,         -1.27785301,  0.,          0.        ],
-                                    [ 0.,          0.,         -0.4482997,   0.        ],
-                                    [ 0.,          0.,          0.,         -0.4482997 ]])
+    ref_one_body_tensor = np.array(
+        [
+            [-1.27785301, 0.0, 0.0, 0.0],
+            [0.0, -1.27785301, 0.0, 0.0],
+            [0.0, 0.0, -0.4482997, 0.0],
+            [0.0, 0.0, 0.0, -0.4482997],
+        ]
+    )
 
     assert np.isclose(h2_ferm_ham[()], 0.7559674441714287)
     assert np.allclose(h2_ferm_ham.one_body_tensor, ref_one_body_tensor)
 
-def test_parse_pauli_string_1():
 
-    pauli_string = ((0, 'X'), (1, 'Y'))
+def test_parse_pauli_string_1():
+    pauli_string = ((0, "X"), (1, "Y"))
     qibo_pauli_string = hamiltonian.parse_pauli_string(pauli_string, 0.5)
-    ref_pauli_string = '0.5*X0*Y1'
+    ref_pauli_string = "0.5*X0*Y1"
     assert str(qibo_pauli_string) == ref_pauli_string
 
-def test_parse_pauli_string_2():
 
+def test_parse_pauli_string_2():
     qibo_pauli_string = hamiltonian.parse_pauli_string(None, 0.1)
-    assert str(qibo_pauli_string) == '0.1'
+    assert str(qibo_pauli_string) == "0.1"
 
-def test_qubit_hamiltonian():
 
+def test_qubit_hamiltonian():
     h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
     h2.run_pyscf()
     h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
     h2_qubit_ham_jw = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "jw")
     h2_qubit_ham_bk = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "bk")
-    ref_h2_qubit_ham_jw = {(): -0.04207897647782238, ((0, 'Z'),): 0.17771287465139918, ((1, 'Z'),): 0.1777128746513992, ((2, 'Z'),): -0.24274280513140478, ((3, 'Z'),): -0.24274280513140478, ((0, 'Z'), (1, 'Z')): 0.17059738328801052, ((0, 'Z'), (2, 'Z')): 0.12293305056183809, ((0, 'Z'), (3, 'Z')): 0.16768319457718972, ((1, 'Z'), (2, 'Z')): 0.16768319457718972, ((1, 'Z'), (3, 'Z')): 0.12293305056183809, ((2, 'Z'), (3, 'Z')): 0.17627640804319608, ((0, 'X'), (1, 'X'), (2, 'Y'), (3, 'Y')): -0.04475014401535165, ((0, 'X'), (1, 'Y'), (2, 'Y'), (3, 'X')): 0.04475014401535165, ((0, 'Y'), (1, 'X'), (2, 'X'), (3, 'Y')): 0.04475014401535165, ((0, 'Y'), (1, 'Y'), (2, 'X'), (3, 'X')): -0.04475014401535165}
-    ref_h2_qubit_ham_bk = {((0, 'Z'),): 0.17771287465139923, (): -0.04207897647782244, ((0, 'Z'), (1, 'Z')): 0.17771287465139918, ((2, 'Z'),): -0.24274280513140484, ((1, 'Z'), (2, 'Z'), (3, 'Z')): -0.24274280513140484, ((0, 'Y'), (1, 'Z'), (2, 'Y')): 0.04475014401535165, ((0, 'X'), (1, 'Z'), (2, 'X')): 0.04475014401535165, ((0, 'X'), (1, 'Z'), (2, 'X'), (3, 'Z')): 0.04475014401535165, ((0, 'Y'), (1, 'Z'), (2, 'Y'), (3, 'Z')): 0.04475014401535165, ((1, 'Z'),): 0.17059738328801052, ((0, 'Z'), (2, 'Z')): 0.12293305056183809, ((0, 'Z'), (1, 'Z'), (2, 'Z')): 0.16768319457718972, ((0, 'Z'), (1, 'Z'), (2, 'Z'), (3, 'Z')): 0.16768319457718972, ((0, 'Z'), (2, 'Z'), (3, 'Z')): 0.12293305056183809, ((1, 'Z'), (3, 'Z')): 0.17627640804319608}
+    ref_h2_qubit_ham_jw = {
+        (): -0.04207897647782238,
+        ((0, "Z"),): 0.17771287465139918,
+        ((1, "Z"),): 0.1777128746513992,
+        ((2, "Z"),): -0.24274280513140478,
+        ((3, "Z"),): -0.24274280513140478,
+        ((0, "Z"), (1, "Z")): 0.17059738328801052,
+        ((0, "Z"), (2, "Z")): 0.12293305056183809,
+        ((0, "Z"), (3, "Z")): 0.16768319457718972,
+        ((1, "Z"), (2, "Z")): 0.16768319457718972,
+        ((1, "Z"), (3, "Z")): 0.12293305056183809,
+        ((2, "Z"), (3, "Z")): 0.17627640804319608,
+        ((0, "X"), (1, "X"), (2, "Y"), (3, "Y")): -0.04475014401535165,
+        ((0, "X"), (1, "Y"), (2, "Y"), (3, "X")): 0.04475014401535165,
+        ((0, "Y"), (1, "X"), (2, "X"), (3, "Y")): 0.04475014401535165,
+        ((0, "Y"), (1, "Y"), (2, "X"), (3, "X")): -0.04475014401535165,
+    }
+    ref_h2_qubit_ham_bk = {
+        ((0, "Z"),): 0.17771287465139923,
+        (): -0.04207897647782244,
+        ((0, "Z"), (1, "Z")): 0.17771287465139918,
+        ((2, "Z"),): -0.24274280513140484,
+        ((1, "Z"), (2, "Z"), (3, "Z")): -0.24274280513140484,
+        ((0, "Y"), (1, "Z"), (2, "Y")): 0.04475014401535165,
+        ((0, "X"), (1, "Z"), (2, "X")): 0.04475014401535165,
+        ((0, "X"), (1, "Z"), (2, "X"), (3, "Z")): 0.04475014401535165,
+        ((0, "Y"), (1, "Z"), (2, "Y"), (3, "Z")): 0.04475014401535165,
+        ((1, "Z"),): 0.17059738328801052,
+        ((0, "Z"), (2, "Z")): 0.12293305056183809,
+        ((0, "Z"), (1, "Z"), (2, "Z")): 0.16768319457718972,
+        ((0, "Z"), (1, "Z"), (2, "Z"), (3, "Z")): 0.16768319457718972,
+        ((0, "Z"), (2, "Z"), (3, "Z")): 0.12293305056183809,
+        ((1, "Z"), (3, "Z")): 0.17627640804319608,
+    }
 
     keylist_jw = h2_qubit_ham_jw.terms.keys()
     keylist_bk = h2_qubit_ham_bk.terms.keys()
-    jw_array = np.array([h2_qubit_ham_jw.terms[key] for key in keylist_jw])    
-    bk_array = np.array([h2_qubit_ham_bk.terms[key] for key in keylist_bk])    
-    
+    jw_array = np.array([h2_qubit_ham_jw.terms[key] for key in keylist_jw])
+    bk_array = np.array([h2_qubit_ham_bk.terms[key] for key in keylist_bk])
+
     keylist_ref_jw = ref_h2_qubit_ham_jw.keys()
     keylist_ref_bk = ref_h2_qubit_ham_bk.keys()
     ref_jw_array = np.array([ref_h2_qubit_ham_jw[key] for key in keylist_ref_jw])
@@ -90,7 +125,7 @@ def test_qubit_hamiltonian():
 
     assert np.allclose(jw_array, ref_jw_array)
     assert np.allclose(bk_array, ref_bk_array)
- 
+
 
 def test_symbolic_hamiltonian():
     h2 = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
@@ -98,11 +133,28 @@ def test_symbolic_hamiltonian():
     h2_ferm_ham = hamiltonian.fermionic_hamiltonian(h2.oei, h2.tei, h2.e_nuc)
     h2_qubit_ham = hamiltonian.qubit_hamiltonian(h2_ferm_ham, "jw")
     h2_sym_ham = hamiltonian.symbolic_hamiltonian(h2_qubit_ham)
-    ref_ham = -0.0420789764778224 - 0.0447501440153516*X(0)*X(1)*Y(2)*Y(3) + 0.0447501440153516*X(0)*Y(1)*Y(2)*X(3) + 0.0447501440153516*Y(0)*X(1)*X(2)*Y(3) - 0.0447501440153516*Y(0)*Y(1)*X(2)*X(3) + 0.177712874651399*Z(0) + 0.170597383288011*Z(0)*Z(1) + 0.122933050561838*Z(0)*Z(2) + 0.16768319457719*Z(0)*Z(3) + 0.177712874651399*Z(1) + 0.16768319457719*Z(1)*Z(2) + 0.122933050561838*Z(1)*Z(3) - 0.242742805131405*Z(2) + 0.176276408043196*Z(2)*Z(3) - 0.242742805131405*Z(3)
+    ref_ham = (
+        -0.0420789764778224
+        - 0.0447501440153516 * X(0) * X(1) * Y(2) * Y(3)
+        + 0.0447501440153516 * X(0) * Y(1) * Y(2) * X(3)
+        + 0.0447501440153516 * Y(0) * X(1) * X(2) * Y(3)
+        - 0.0447501440153516 * Y(0) * Y(1) * X(2) * X(3)
+        + 0.177712874651399 * Z(0)
+        + 0.170597383288011 * Z(0) * Z(1)
+        + 0.122933050561838 * Z(0) * Z(2)
+        + 0.16768319457719 * Z(0) * Z(3)
+        + 0.177712874651399 * Z(1)
+        + 0.16768319457719 * Z(1) * Z(2)
+        + 0.122933050561838 * Z(1) * Z(3)
+        - 0.242742805131405 * Z(2)
+        + 0.176276408043196 * Z(2) * Z(3)
+        - 0.242742805131405 * Z(3)
+    )
     ref_sym_ham = SymbolicHamiltonian(ref_ham)
 
     assert np.allclose(h2_sym_ham.matrix, ref_sym_ham.matrix)
 
+
 @pytest.mark.skip(reason="psi4 doesn't offer pip install, so needs to be installed through conda or manually.")
 def test_run_psi4():
     """PSI4 driver"""