diff --git a/codes/mf.py b/codes/mf.py
index 8d3efb7f3238a6050e8b88c7608f22635f1dd690..e3f9c3a5816040cb82eb3ceee2969abe6047179e 100644
--- a/codes/mf.py
+++ b/codes/mf.py
@@ -26,37 +26,6 @@ def density_matrix(kham, fermi_energy):
     return density_matrix_kgrid
 
 
-def fermi_on_grid(kham, filling):
-    """
-     Compute the Fermi energy on a grid of k-points.
-
-     Parameters
-     ----------
-     hkfunc : function
-         Function that returns the Hamiltonian at a given k-point.
-     nk : int
-         Number of k-points in the grid.
-     Returns
-     -------
-    fermi_energy : float
-         Fermi energy
-    """
-
-    vals = np.linalg.eigvalsh(kham)
-
-    norbs = vals.shape[-1]
-    vals_flat = np.sort(vals.flatten())
-    ne = len(vals_flat)
-    ifermi = int(round(ne * filling / norbs))
-    if ifermi >= ne:
-        return vals_flat[-1]
-    elif ifermi == 0:
-        return vals_flat[0]
-    else:
-        fermi = (vals_flat[ifermi - 1] + vals_flat[ifermi]) / 2
-        return fermi
-
-
 def meanfield(density_matrix_tb, h_int, n=2):
     """
     Compute the mean-field in k-space.
@@ -95,3 +64,34 @@ def meanfield(density_matrix_tb, h_int, n=2):
         for vec in frozenset(h_int)
     }
     return add_tb(direct, exchange)
+
+
+def fermi_on_grid(kham, filling):
+    """
+     Compute the Fermi energy on a grid of k-points.
+
+     Parameters
+     ----------
+     hkfunc : function
+         Function that returns the Hamiltonian at a given k-point.
+     nk : int
+         Number of k-points in the grid.
+     Returns
+     -------
+    fermi_energy : float
+         Fermi energy
+    """
+
+    vals = np.linalg.eigvalsh(kham)
+
+    norbs = vals.shape[-1]
+    vals_flat = np.sort(vals.flatten())
+    ne = len(vals_flat)
+    ifermi = int(round(ne * filling / norbs))
+    if ifermi >= ne:
+        return vals_flat[-1]
+    elif ifermi == 0:
+        return vals_flat[0]
+    else:
+        fermi = (vals_flat[ifermi - 1] + vals_flat[ifermi]) / 2
+        return fermi
diff --git a/codes/params/test_params.py b/codes/params/test_params.py
index 66a31da1e03043cc7a2782fb92a7b63112c1ca56..c97d3f18dcb889613dcede54bd4ee4940e419414 100644
--- a/codes/params/test_params.py
+++ b/codes/params/test_params.py
@@ -1,6 +1,6 @@
 # %%
 from codes.params.rparams import mf_to_rparams, rparams_to_mf
-from codes.kwant_helper.utils import generate_guess
+from codes.tb.utils import generate_guess
 from codes.tb.tb import compare_dicts
 import pytest
 
diff --git a/codes/tb/test_tb.py b/codes/tb/test_tb.py
index 6fee47cee4af4e21522cf59b5f02b210c287d16a..bb42c6251c19de16a51c93c7c32ae9e25bfd4be1 100644
--- a/codes/tb/test_tb.py
+++ b/codes/tb/test_tb.py
@@ -1,8 +1,8 @@
 # %%
 import numpy as np
 from codes.tb.tb import compare_dicts
-from codes.kwant_helper import utils
 import itertools as it
+from codes.tb.utils import generate_guess
 from codes.tb.transforms import kfunc_to_tb, tb_to_kfunc, tb_to_kham, tb_to_khamvector
 import pytest
 
@@ -39,7 +39,7 @@ def test_tbkham_transform():
         (-1, -1),
     )
     ndof = 10
-    h_0 = utils.generate_guess(vectors, ndof)
+    h_0 = generate_guess(vectors, ndof)
 
     assert np.allclose(
         tb_to_kham(h_0, nk=nk, ndim=2), tb_to_khamvector(h_0, nk=nk, ndim=2)
diff --git a/codes/test_graphene.py b/codes/test_graphene.py
index 921428ef58b304e21423c1a7fef4361a9b277d3a..9ff39c0d99e08056f21bfe90d174bdb9c8b4003e 100644
--- a/codes/test_graphene.py
+++ b/codes/test_graphene.py
@@ -4,7 +4,7 @@ from codes.model import Model
 from codes.solvers import solver
 from codes import kwant_examples
 from codes.kwant_helper import utils
-from codes.tb.utils import compute_gap
+from codes.tb.utils import compute_gap, generate_guess
 from codes.tb.tb import add_tb
 import pytest
 
@@ -46,7 +46,7 @@ def gap_prediction(U, V):
     nk = 20
 
     h_int = utils.builder_to_tb(int_builder, params)
-    guess = utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
+    guess = generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
     model = Model(h_0, h_int, filling)
 
     mf_sol = solver(model, guess, nk=nk, optimizer_kwargs={"verbose": True, "M": 0})