diff --git a/codes/kwant_helper/utils.py b/codes/kwant_helper/utils.py
index 9b537edcec4001da6dcabf4ad6dc54ee6a977232..88e6261136af87e1e5e2e0e076821906cf947acf 100644
--- a/codes/kwant_helper/utils.py
+++ b/codes/kwant_helper/utils.py
@@ -146,93 +146,6 @@ def builder2h_0(builder, params={}, return_data=False):
         return h_0
 
 
-def model2hk(h_0):
-    """
-    Build Bloch Hamiltonian.
-
-    Paramters:
-    ----------
-    nk : int
-        Number of k-points along each direction.
-    h_0 : dictionary
-        Must have the following structure:
-            - Keys are tuples for each hopping vector (in units of lattice vectors).
-            - Values are hopping matrices.
-    return_ks : bool
-        Return k-points.
-
-    Returns:
-    --------
-    ham : nd.array
-        Hamiltonian evaluated on a k-point grid from k-points
-        along each direction evaluated from zero to 2*pi.
-        The indices are ordered as [k_1, ... , k_n, i, j], where
-        `k_m` corresponding to the k-point element along each
-        direction and `i` and `j` are the internal degrees of freedom.
-    ks : 1D-array
-        List of k-points over all directions. Only returned if `return_ks=True`.
-
-    Returns:
-    --------
-    bloch_ham : function
-        Evaluates the Hamiltonian at a given k-point.
-    """
-    assert (
-        len(next(iter(h_0))) > 0
-    ), "Zero-dimensional system. The Hamiltonian is simply h_0[()]."
-
-    def bloch_ham(k):
-        ham = 0
-        for vector in h_0.keys():
-            ham += h_0[vector] * np.exp(-1j * np.dot(k, np.array(vector, dtype=float)))
-        return ham
-
-    return bloch_ham
-
-
-def kgrid_hamiltonian(nk, hk, dim, return_ks=False, hermitian=True):
-    """
-    Evaluates Hamiltonian on a k-point grid.
-
-    Paramters:
-    ----------
-    nk : int
-        Number of k-points along each direction.
-    hk : function
-        Calculates the Hamiltonian at a given k-point.
-    return_ks : bool
-        If `True`, returns k-points.
-
-    Returns:
-    --------
-    ham : nd.array
-        Hamiltonian evaluated on a k-point grid from k-points
-        along each direction evaluated from zero to 2*pi.
-        The indices are ordered as [k_1, ... , k_n, i, j], where
-        `k_m` corresponding to the k-point element along each
-        direction and `i` and `j` are the internal degrees of freedom.
-    ks : 1D-array
-        List of k-points over all directions. Only returned if `return_ks=True`.
-    """
-    ks = 2 * np.pi * np.linspace(0, 1, nk, endpoint=False)
-
-    k_pts = np.tile(ks, dim).reshape(dim, nk)
-
-    ham = []
-    for k in product(*k_pts):
-        ham.append(hk(k))
-    ham = np.array(ham)
-    if hermitian:
-        assert np.allclose(
-            ham, np.transpose(ham, (0, 2, 1)).conj()
-        ), "Tight-binding provided is non-Hermitian. Not supported yet"
-    shape = (*[nk] * dim, ham.shape[-1], ham.shape[-1])
-    if return_ks:
-        return ham.reshape(*shape), ks
-    else:
-        return ham.reshape(*shape)
-
-
 def build_interacting_syst(builder, lattice, func_onsite, func_hop, max_neighbor=1):
     """
     Construct an auxiliary `kwant` system to build Hamiltonian matrix.
@@ -311,53 +224,6 @@ def generate_vectors(cutoff, dim):
     return [*product(*([[*range(-cutoff, cutoff + 1)]] * dim))]
 
 
-def hk2h_0(hk, hopping_vecs, ks=None):
-    """
-    Extract hopping matrices from Bloch Hamiltonian.
-
-    Parameters:
-    -----------
-    hk : nd-array
-        Bloch Hamiltonian matrix hk[k_x, ..., k_n, i, j]
-    h_0 : dict
-        Tight-binding model of non-interacting systems.
-    h_int : dict
-        Tight-binding model for interacting Hamiltonian.
-    ks : 1D-array
-        Set of k-points. Repeated for all directions. If the system is finite, `ks=None`.
-
-    Returns:
-    --------
-    scf_model : dict
-        TIght-binding model of Hartree-Fock solution.
-    """
-    if ks is not None:
-        ndim = len(hk.shape) - 2
-        dk = np.diff(ks)[0]
-        nk = len(ks)
-        k_pts = np.tile(ks, ndim).reshape(ndim, nk)
-        k_grid = np.array(np.meshgrid(*k_pts))
-        k_grid = k_grid.reshape(k_grid.shape[0], np.prod(k_grid.shape[1:]))
-        hk = hk.reshape(np.prod(hk.shape[:ndim]), *hk.shape[-2:])
-
-        hopps = (
-            np.einsum(
-                "ij,jkl->ikl",
-                np.exp(1j * np.einsum("ij,jk->ik", hopping_vecs, k_grid)),
-                hk,
-            )
-            * (dk / (2 * np.pi)) ** ndim
-        )
-
-        h_0 = {}
-        for i, vector in enumerate(hopping_vecs):
-            h_0[tuple(vector)] = hopps[i]
-
-        return h_0
-    else:
-        return {(): hk}
-
-
 def calc_gap(vals, E_F):
     """
     Compute gap.
@@ -377,35 +243,3 @@ def calc_gap(vals, E_F):
     emax = np.max(vals[vals <= E_F])
     emin = np.min(vals[vals > E_F])
     return np.abs(emin - emax)
-
-
-def matrix_to_flat(matrix):
-    """
-    Flatten the upper triangle of a collection of matrices.
-
-    Parameters:
-    -----------
-    matrix : nd-array
-        Array with shape (..., n, n)
-    """
-    return matrix[..., *np.triu_indices(matrix.shape[-1])].flatten()
-
-
-def flat_to_matrix(flat, shape):
-    """
-    Undo `matrix_to_flat`.
-
-    Parameters:
-    -----------
-    flat : 1d-array
-        Output from `matrix_to_flat`.
-    shape : 1d-array
-        Shape of the input from `matrix_to_flat`.
-    """
-    matrix = np.zeros(shape, dtype=complex)
-    matrix[..., *np.triu_indices(shape[-1])] = flat.reshape(*shape[:-2], -1)
-    indices = np.arange(shape[-1])
-    diagonal = matrix[..., indices, indices]
-    matrix += np.moveaxis(matrix, -1, -2).conj()
-    matrix[..., indices, indices] -= diagonal
-    return matrix
\ No newline at end of file
diff --git a/codes/tb/transforms.py b/codes/tb/transforms.py
index 4bfe6b398fa6b2f6faff2d50ce11a69a8120f6bc..dbd64db699af5d78adf4804e810c01bbfc9a5594 100644
--- a/codes/tb/transforms.py
+++ b/codes/tb/transforms.py
@@ -171,3 +171,50 @@ def kfunc2tb(kfunc, nSamples, ndim=1):
         raise NotImplementedError("n > 2 not implemented")
     ifftnArray = ifftn(kfuncOnGrid, axes=np.arange(ndim))
     return ifftn2tb(ifftnArray)
+
+
+def hk2h_0(hk, hopping_vecs, ks=None):
+    """
+    Extract hopping matrices from Bloch Hamiltonian.
+
+    Parameters:
+    -----------
+    hk : nd-array
+        Bloch Hamiltonian matrix hk[k_x, ..., k_n, i, j]
+    h_0 : dict
+        Tight-binding model of non-interacting systems.
+    h_int : dict
+        Tight-binding model for interacting Hamiltonian.
+    ks : 1D-array
+        Set of k-points. Repeated for all directions. If the system is finite, `ks=None`.
+
+    Returns:
+    --------
+    scf_model : dict
+        TIght-binding model of Hartree-Fock solution.
+    """
+    if ks is not None:
+        ndim = len(hk.shape) - 2
+        dk = np.diff(ks)[0]
+        nk = len(ks)
+        k_pts = np.tile(ks, ndim).reshape(ndim, nk)
+        k_grid = np.array(np.meshgrid(*k_pts))
+        k_grid = k_grid.reshape(k_grid.shape[0], np.prod(k_grid.shape[1:]))
+        hk = hk.reshape(np.prod(hk.shape[:ndim]), *hk.shape[-2:])
+
+        hopps = (
+            np.einsum(
+                "ij,jkl->ikl",
+                np.exp(1j * np.einsum("ij,jk->ik", hopping_vecs, k_grid)),
+                hk,
+            )
+            * (dk / (2 * np.pi)) ** ndim
+        )
+
+        h_0 = {}
+        for i, vector in enumerate(hopping_vecs):
+            h_0[tuple(vector)] = hopps[i]
+
+        return h_0
+    else:
+        return {(): hk}