diff --git a/codes/model.py b/codes/model.py
index 0f7086fab2e7bd7c06180700f529e7f37e356003..70c06bc3d47504d740fd321f591be7d4bcf1af1f 100644
--- a/codes/model.py
+++ b/codes/model.py
@@ -33,18 +33,20 @@ class Model:
         _check_hermiticity(h_int)
 
     def calculate_EF(self):
-        self.EF = fermi_on_grid(self.kham, self.filling)
+        self.fermi_energy = fermi_on_grid(self.kham, self.filling)
 
-    def make_density_matrix_tb(self, mf_model, nk=200):
-        self.kham = tb_to_khamvector(add_tb(self.h_0, mf_model), nk=nk, ndim=self._ndim)
+    def make_density_matrix_tb(self, mf_tb, nk=200):
+        self.kham = tb_to_khamvector(add_tb(self.h_0, mf_tb), nk=nk, ndim=self._ndim)
         self.calculate_EF()
         return ifftn_to_tb(
-            ifftn(density_matrix(self.kham, self.EF), axes=np.arange(self._ndim))
+            ifftn(
+                density_matrix(self.kham, self.fermi_energy), axes=np.arange(self._ndim)
+            )
         )
 
-    def mfield(self, mf_model, nk=200):
-        density_matrix_tb = self.make_density_matrix_tb(mf_model, nk=nk)
+    def mfield(self, mf_tb, nk=200):
+        density_matrix_tb = self.make_density_matrix_tb(mf_tb, nk=nk)
         return add_tb(
             meanfield(density_matrix_tb, self.h_int, n=self._ndim),
-            {self._local_key: -self.EF * np.eye(self._size)},
+            {self._local_key: -self.fermi_energy * np.eye(self._size)},
         )
diff --git a/codes/params/param_transforms.py b/codes/params/param_transforms.py
index b681de8617a2c044664cba7437f040ea4efa8225..d236b73be1d4a507f9f98c9d357ee1471cf804c8 100644
--- a/codes/params/param_transforms.py
+++ b/codes/params/param_transforms.py
@@ -1,13 +1,13 @@
 import numpy as np
 
 
-def tb_to_flat(hop_dict):
+def tb_to_flat(tb):
     """
     Convert a hermitian tight-binding dictionary to flat complex matrix.
 
     Parameters:
     -----------
-    hop_dict : dict with nd-array elements
+    tb : dict with nd-array elements
         Hermitian tigh-binding dictionary
 
     Returns:
@@ -15,14 +15,12 @@ def tb_to_flat(hop_dict):
     flat : complex 1d numpy array
         Flattened tight-binding dictionary
     """
-    N = len(hop_dict.keys()) // 2 + 1
-    sorted_vals = np.array(list(hop_dict.values()))[
-        np.lexsort(np.array(list(hop_dict.keys())).T)
-    ]
+    N = len(tb.keys()) // 2 + 1
+    sorted_vals = np.array(list(tb.values()))[np.lexsort(np.array(list(tb.keys())).T)]
     return sorted_vals[:N].flatten()
 
 
-def flat_to_tb(flat, shape, hop_dict_keys):
+def flat_to_tb(flat, shape, tb_keys):
     """
     Reverse operation to `tb_to_flat` that takes a flat complex 1d array
     and return the tight-binding dictionary.
@@ -32,24 +30,24 @@ def flat_to_tb(flat, shape, hop_dict_keys):
     flat : dict with nd-array elements
         Hermitian tigh-binding dictionary
     shape : tuple
-        shape of the hop_dict elements
-    hop_dict_keys : iterable
-        original hop_dict key elements
+        shape of the tb elements
+    tb_keys : iterable
+        original tb key elements
 
     Returns:
     -----------
-    hop_dict : dict
+    tb : dict
         tight-binding dictionary
     """
     matrix = np.zeros(shape, dtype=complex)
-    N = len(hop_dict_keys) // 2 + 1
+    N = len(tb_keys) // 2 + 1
     matrix[:N] = flat.reshape(N, *shape[1:])
     matrix[N:] = np.moveaxis(matrix[-(N + 1) :: -1], -1, -2).conj()
 
-    hop_dict_keys = np.array(list(hop_dict_keys))
-    sorted_keys = hop_dict_keys[np.lexsort(hop_dict_keys.T)]
-    hop_dict = dict(zip(map(tuple, sorted_keys), matrix))
-    return hop_dict
+    tb_keys = np.array(list(tb_keys))
+    sorted_keys = tb_keys[np.lexsort(tb_keys.T)]
+    tb = dict(zip(map(tuple, sorted_keys), matrix))
+    return tb
 
 
 def complex_to_real(z):
diff --git a/codes/params/rparams.py b/codes/params/rparams.py
index ffb770e23e8fbd8cd7c1481fb6f588ee65001ea2..7acbbe3259415e61e3e8c059b244e265df29a3bc 100644
--- a/codes/params/rparams.py
+++ b/codes/params/rparams.py
@@ -1,4 +1,4 @@
-from codes.params.matrixShaping import (
+from codes.params.param_transforms import (
     complex_to_real,
     tb_to_flat,
     real_to_complex,
@@ -6,13 +6,13 @@ from codes.params.matrixShaping import (
 )
 
 
-def mf_to_rparams(mf_model):
+def tb_to_rparams(tb):
     """
     Convert a mean-field tight-binding model to a set of real parameters.
 
     Parameters
     ----------
-    mf_model : dict
+    tb : dict
         Mean-field tight-binding model.
 
     Returns
@@ -20,10 +20,10 @@ def mf_to_rparams(mf_model):
     dict
         Real parameters.
     """
-    return complex_to_real(tb_to_flat(mf_model))  # placeholder for now
+    return complex_to_real(tb_to_flat(tb))  # placeholder for now
 
 
-def rparams_to_mf(rParams, key_list, size):
+def rparams_to_tb(rParams, key_list, size):
     """
     Extract mean-field tight-binding model from a set of real parameters.
 
diff --git a/codes/params/test_params.py b/codes/params/test_params.py
index c97d3f18dcb889613dcede54bd4ee4940e419414..9fcf1bff1fd3c754ae03ad775c4d18f4cf6b64d6 100644
--- a/codes/params/test_params.py
+++ b/codes/params/test_params.py
@@ -1,5 +1,5 @@
 # %%
-from codes.params.rparams import mf_to_rparams, rparams_to_mf
+from codes.params.rparams import tb_to_rparams, rparams_to_tb
 from codes.tb.utils import generate_guess
 from codes.tb.tb import compare_dicts
 import pytest
@@ -14,6 +14,6 @@ vectors = ((0, 0), (1, 0), (-1, 0), (0, 1), (0, -1), (1, -1), (-1, 1), (1, 1), (
 @pytest.mark.repeat(repeat_number)
 def test_parametrisation():
     mf_guess = generate_guess(vectors, ndof)
-    mf_params = mf_to_rparams(mf_guess)
-    mf_new = rparams_to_mf(mf_params, vectors, ndof)
+    mf_params = tb_to_rparams(mf_guess)
+    mf_new = rparams_to_tb(mf_params, vectors, ndof)
     compare_dicts(mf_guess, mf_new)
diff --git a/codes/solvers.py b/codes/solvers.py
index ce83449348c37550ef3c93d57b8b1f2d96f4bc0d..9bbd4b1eaf5e41d6d54c6d5827dfcf9450c060d2 100644
--- a/codes/solvers.py
+++ b/codes/solvers.py
@@ -1,4 +1,4 @@
-from codes.params.rparams import mf_to_rparams, rparams_to_mf
+from codes.params.rparams import tb_to_rparams, rparams_to_tb
 import scipy
 from functools import partial
 from codes.tb.tb import add_tb
@@ -21,9 +21,9 @@ def cost(mf_param, Model, nk=100):
         The number of k-points to use in the grid. The default is 100.
     """
     shape = Model._size
-    mf_tb = rparams_to_mf(mf_param, list(Model.h_int), shape)
+    mf_tb = rparams_to_tb(mf_param, list(Model.h_int), shape)
     mf_tb_new = Model.mfield(mf_tb, nk=nk)
-    mf_params_new = mf_to_rparams(mf_tb_new)
+    mf_params_new = tb_to_rparams(mf_tb_new)
     return mf_params_new - mf_param
 
 
@@ -53,11 +53,11 @@ def solver(
     """
 
     shape = Model._size
-    mf_params = mf_to_rparams(mf_guess)
+    mf_params = tb_to_rparams(mf_guess)
     f = partial(cost, Model=Model, nk=nk)
-    result = rparams_to_mf(
+    result = rparams_to_tb(
         optimizer(f, mf_params, **optimizer_kwargs), list(Model.h_int), shape
     )
     Model.calculate_EF()
     local_key = tuple(np.zeros((Model._ndim,), dtype=int))
-    return add_tb(result, {local_key: -Model.EF * np.eye(shape)})
+    return add_tb(result, {local_key: -Model.fermi_energy * np.eye(shape)})
diff --git a/codes/tb/transforms.py b/codes/tb/transforms.py
index 6df1adad1dc85ec2220fb9a4748e8d1ea378290d..288a412bc03bec17c9b7c7242ff6d8551e62a8c8 100644
--- a/codes/tb/transforms.py
+++ b/codes/tb/transforms.py
@@ -3,13 +3,13 @@ from scipy.fftpack import ifftn
 import itertools as it
 
 
-def tb_to_khamvector(h_0, nk, ndim):
+def tb_to_khamvector(tb, nk, ndim):
     """
     Real-space tight-binding model to hamiltonian on k-space grid.
 
     Parameters
     ----------
-    h_0 : dict
+    tb : dict
         A dictionary with real-space vectors as keys and complex np.arrays as values.
     nk : int
         Number of k-points along each direction.
@@ -27,9 +27,9 @@ def tb_to_khamvector(h_0, nk, ndim):
     ks = np.concatenate((ks[nk // 2 :], ks[: nk // 2]), axis=0)  # shift for ifft
     kgrid = np.meshgrid(ks, ks, indexing="ij")
 
-    num_keys = len(list(h_0.keys()))
-    tb_array = np.array(list(h_0.values()))
-    keys = np.array(list(h_0.keys()))
+    num_keys = len(list(tb.keys()))
+    tb_array = np.array(list(tb.values()))
+    keys = np.array(list(tb.keys()))
 
     k_dependency = np.exp(-1j * np.tensordot(keys, kgrid, 1))[
         (...,) + (np.newaxis,) * 2
@@ -40,13 +40,13 @@ def tb_to_khamvector(h_0, nk, ndim):
     return np.sum(tb_array * k_dependency, axis=0)
 
 
-def tb_to_kfunc(h_0):
+def tb_to_kfunc(tb):
     """
     Fourier transforms a real-space tight-binding model to a k-space function.
 
     Parameters
     ----------
-    h_0 : dict
+    tb : dict
         A dictionary with real-space vectors as keys and complex np.arrays as values.
 
     Returns
@@ -55,13 +55,13 @@ def tb_to_kfunc(h_0):
         A function that takes a k-space vector and returns a complex np.array.
     """
 
-    def bloch_ham(k):
+    def kfunc(k):
         ham = 0
-        for vector in h_0.keys():
-            ham += h_0[vector] * np.exp(-1j * np.dot(k, np.array(vector, dtype=float)))
+        for vector in tb.keys():
+            ham += tb[vector] * np.exp(-1j * np.dot(k, np.array(vector, dtype=float)))
         return ham
 
-    return bloch_ham
+    return kfunc
 
 
 def ifftn_to_tb(ifft_array):
@@ -117,19 +117,19 @@ def kfunc_to_kham(nk, hk, dim, return_ks=False, hermitian=True):
 
     k_pts = np.tile(ks, dim).reshape(dim, nk)
 
-    ham = []
+    kham = []
     for k in it.product(*k_pts):
-        ham.append(hk(k))
-    ham = np.array(ham)
+        kham.append(hk(k))
+    kham = np.array(kham)
     if hermitian:
         assert np.allclose(
-            ham, np.transpose(ham, (0, 2, 1)).conj()
+            kham, np.transpose(kham, (0, 2, 1)).conj()
         ), "Tight-binding provided is non-Hermitian. Not supported yet"
-    shape = (*[nk] * dim, ham.shape[-1], ham.shape[-1])
+    shape = (*[nk] * dim, kham.shape[-1], kham.shape[-1])
     if return_ks:
-        return ham.reshape(*shape), ks
+        return kham.reshape(*shape), ks
     else:
-        return ham.reshape(*shape)
+        return kham.reshape(*shape)
 
 
 def tb_to_kham(h_0, nk=200, ndim=1):
@@ -173,14 +173,14 @@ def kfunc_to_tb(kfunc, n_samples, ndim=1):
     return ifftn_to_tb(ifftn_array)
 
 
-def hk_to_h0(hk, hopping_vecs, ks=None):
+def kham_to_tb(kham, hopping_vecs, ks=None):
     """
     Extract hopping matrices from Bloch Hamiltonian.
 
     Parameters:
     -----------
-    hk : nd-array
-        Bloch Hamiltonian matrix hk[k_x, ..., k_n, i, j]
+    kham : nd-array
+        Bloch Hamiltonian matrix kham[k_x, ..., k_n, i, j]
     h_0 : dict
         Tight-binding model of non-interacting systems.
     h_int : dict
@@ -191,22 +191,22 @@ def hk_to_h0(hk, hopping_vecs, ks=None):
     Returns:
     --------
     scf_model : dict
-        TIght-binding model of Hartree-Fock solution.
+        Tight-binding model of Hartree-Fock solution.
     """
     if ks is not None:
-        ndim = len(hk.shape) - 2
+        ndim = len(kham.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:])
+        kham = kham.reshape(np.prod(kham.shape[:ndim]), *kham.shape[-2:])
 
         hopps = (
             np.einsum(
                 "ij,jkl->ikl",
                 np.exp(1j * np.einsum("ij,jk->ik", hopping_vecs, k_grid)),
-                hk,
+                kham,
             )
             * (dk / (2 * np.pi)) ** ndim
         )
@@ -217,4 +217,4 @@ def hk_to_h0(hk, hopping_vecs, ks=None):
 
         return h_0
     else:
-        return {(): hk}
+        return {(): kham}