diff --git a/src/gaspype/_numerics.py b/src/gaspype/_numerics.py
index 2a136aa..dc4717a 100644
--- a/src/gaspype/_numerics.py
+++ b/src/gaspype/_numerics.py
@@ -8,7 +8,7 @@ def null_space(A: FloatArray) -> FloatArray:
 
     Args:
         A: Input matrix of shape (m, n)
-    
+
     Return:
         Null space vectors as columns, shape (n, n - rank)
     """
@@ -17,5 +17,5 @@ def null_space(A: FloatArray) -> FloatArray:
     rcond = np.finfo(s.dtype).eps * max(M, N)
     tol = np.amax(s, initial=0.) * rcond
     num = np.sum(s > tol, dtype=int)
-    Q = vh[num:,:].T.conj()
-    return Q
\ No newline at end of file
+    Q = vh[num:, :].T.conj()
+    return Q
diff --git a/tests/benchmark_equalibrium.py b/tests/benchmark_equalibrium.py
index 748e98d..4c0fe3a 100644
--- a/tests/benchmark_equalibrium.py
+++ b/tests/benchmark_equalibrium.py
@@ -44,6 +44,10 @@ eq_gaspype = gp.equilibrium(fluid, t=temperatures, p=pressure)
 elapsed_gaspype = time.perf_counter() - t0
 print(f"Gaspype: {elapsed_gaspype:.4f} s")
 
+# Check if elemental balance of result is correct
+el_err = np.sum((gp.elements(eq_gaspype) - gp.elements(fluid)).get_n()**2)
+assert np.all(el_err < 1e-20)
+
 # -----------------------
 # Compare first 5 results
 # -----------------------
diff --git a/tests/test_equalibrium.py b/tests/test_equalibrium.py
new file mode 100644
index 0000000..26d9eb3
--- /dev/null
+++ b/tests/test_equalibrium.py
@@ -0,0 +1,25 @@
+import gaspype as gp
+
+
+def test_single_equilibrium():
+    # Compare equilibrium calculations to Cantera results
+
+    # gp.set_solver('system of equations')
+    # gp.set_solver('gibs minimization')
+
+    # fs = gp.fluid_system(['CH4', 'C2H6', 'C3H8', 'H2O', 'H2', 'CO2', 'CO', 'O2'])
+    fs = gp.fluid_system(['CH4', 'H2O', 'H2', 'CO2', 'CO', 'O2'])
+    # fs = gp.fluid_system([s for s in flow1.species_names if s in gps])
+
+    composition = gp.elements({'H': 2, 'O': 0, 'C': 0}, fs)
+
+    t = 1495 + 273.15  # K
+    p = 1e5  # Pa
+
+    fl = gp.equilibrium(composition, t, p)
+
+    print(fl)
+
+
+if __name__ == "__main__":
+    test_single_equilibrium()
diff --git a/tests/test_results.py b/tests/test_results.py
index 165f21a..50a4f9c 100644
--- a/tests/test_results.py
+++ b/tests/test_results.py
@@ -55,6 +55,7 @@ def test_nh3_data():
 def test_equilibrium():
     # Compare equilibrium calculations to Cycle-Tempo results
     df = pd.read_csv('tests/test_data/cycle_temp_matlab_ref.csv', sep=';', decimal=',').fillna(0)
+    #gp.set_solver('gibs minimization')
     fs = gp.fluid_system(['CH4', 'C2H6', 'C3H8', 'C4H10,n-butane', 'H2O', 'H2', 'CO2', 'CO'])
 
     for index, row in df.iterrows():
@@ -84,7 +85,7 @@ def test_equilibrium():
             result_values = gp.equilibrium(fl, t, p).array_fractions
 
             print(index, gp.get_solver(), '----')
-            print(molar_comp)
+            print('Species:        ' + ''.join(f"{s:14}" for s in fs.species))
             outp(result_values, 'Under test: ')
             outp(reference_values, 'Reference:  ')
 
@@ -123,3 +124,7 @@ def test_carbon():
             assert result_values > 0.9
         else:
             assert result_values < 1.1
+
+
+if __name__ == '__main__':
+    test_equilibrium()
diff --git a/tests/test_results_cantera.py b/tests/test_results_cantera.py
index 39ffee6..018c6e6 100644
--- a/tests/test_results_cantera.py
+++ b/tests/test_results_cantera.py
@@ -50,7 +50,7 @@ def test_equilibrium_cantera():
         print(f"Inp. Composition: {comp_dict}")
         print(f"Res. Composition: {gp_comp_result}")
         print(f"Ref. Composition: {ct_comp_result}")
-        print(f"---")
+        print("---")
         assert np.all(dev < 0.04), f"Deviateion: {dev}"
 
     assert np.mean(deviations) < 2e-4
@@ -74,4 +74,4 @@ def test_cantera():
 
 
 if __name__ == "__main__":
-    test_cantera()
\ No newline at end of file
+    test_equilibrium_cantera()