diff --git a/codes/kwant_helper/utils.py b/codes/kwant_helper/utils.py
index 6e74a4e9f237cb564c93edd97d54a8dd3dbb3f66..9f7e2149d42a2c612a8688a9372b32e054f264e0 100644
--- a/codes/kwant_helper/utils.py
+++ b/codes/kwant_helper/utils.py
@@ -202,22 +202,22 @@ def generate_vectors(cutoff, dim):
     return [*product(*([[*range(-cutoff, cutoff + 1)]] * dim))]
 
 
-def calc_gap(vals, E_F):
+def calc_gap(vals, fermi_energy):
     """
-    Compute gap.
-
-    Parameters:
-    -----------
-    vals : nd-array
-        Eigenvalues on a k-point grid.
-    E_F : float
-        Fermi energy.
-
-    Returns:
-    --------
-    gap : float
-        Indirect gap.
+     Compute gap.
+
+     Parameters:
+     -----------
+     vals : nd-array
+         Eigenvalues on a k-point grid.
+    fermi_energy : float
+         Fermi energy.
+
+     Returns:
+     --------
+     gap : float
+         Indirect gap.
     """
-    emax = np.max(vals[vals <= E_F])
-    emin = np.min(vals[vals > E_F])
+    emax = np.max(vals[vals <= fermi_energy])
+    emin = np.min(vals[vals > fermi_energy])
     return np.abs(emin - emax)
diff --git a/codes/mf.py b/codes/mf.py
index 6ae70ed164f0ed1450944659caba75b21bf5cd3a..8d3efb7f3238a6050e8b88c7608f22635f1dd690 100644
--- a/codes/mf.py
+++ b/codes/mf.py
@@ -1,45 +1,45 @@
 import numpy as np
-from codes.tb.tb import addTb
+from codes.tb.tb import add_tb
 
 
-def densityMatrix(kham, E_F):
+def density_matrix(kham, fermi_energy):
     """
-    Parameters
-    ----------
-    kham : npndarray
-        Hamiltonian in k-space of shape (len(dim), norbs, norbs)
+     Parameters
+     ----------
+     kham : npndarray
+         Hamiltonian in k-space of shape (len(dim), norbs, norbs)
 
-    E_F : float
-        Fermi level
+    fermi_energy : float
+         Fermi level
 
-    Returns
-    -------
-    densityMatrixKgrid : np.ndarray
-        Density matrix in k-space.
+     Returns
+     -------
+     density_matrix_kgrid : np.ndarray
+         Density matrix in k-space.
 
     """
     vals, vecs = np.linalg.eigh(kham)
-    unocc_vals = vals > E_F
+    unocc_vals = vals > fermi_energy
     occ_vecs = vecs
     np.moveaxis(occ_vecs, -1, 2)[unocc_vals, :] = 0
-    densityMatrixKgrid = occ_vecs @ np.moveaxis(occ_vecs, -1, -2).conj()
-    return densityMatrixKgrid
+    density_matrix_kgrid = occ_vecs @ np.moveaxis(occ_vecs, -1, -2).conj()
+    return density_matrix_kgrid
 
 
-def fermiOnGrid(kham, filling):
+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
-    -------
-    E_F : float
-        Fermi energy
+     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)
@@ -57,13 +57,13 @@ def fermiOnGrid(kham, filling):
         return fermi
 
 
-def meanField(densityMatrixTb, h_int, n=2):
+def meanfield(density_matrix_tb, h_int, n=2):
     """
     Compute the mean-field in k-space.
 
     Parameters
     ----------
-    densityMatrix : dict
+    density_matrix : dict
         Density matrix in real-space tight-binding format.
     int_model : dict
         Interaction tb model.
@@ -74,14 +74,14 @@ def meanField(densityMatrixTb, h_int, n=2):
         Mean-field tb model.
     """
 
-    localKey = tuple(np.zeros((n,), dtype=int))
+    local_key = tuple(np.zeros((n,), dtype=int))
 
     direct = {
-        localKey: np.sum(
+        local_key: np.sum(
             np.array(
                 [
                     np.diag(
-                        np.einsum("pp,pn->n", densityMatrixTb[localKey], h_int[vec])
+                        np.einsum("pp,pn->n", density_matrix_tb[local_key], h_int[vec])
                     )
                     for vec in frozenset(h_int)
                 ]
@@ -91,7 +91,7 @@ def meanField(densityMatrixTb, h_int, n=2):
     }
 
     exchange = {
-        vec: -1 * h_int.get(vec, 0) * densityMatrixTb[vec]  # / (2 * np.pi)#**2
+        vec: -1 * h_int.get(vec, 0) * density_matrix_tb[vec]  # / (2 * np.pi)#**2
         for vec in frozenset(h_int)
     }
-    return addTb(direct, exchange)
+    return add_tb(direct, exchange)
diff --git a/codes/model.py b/codes/model.py
index 60c7e697837a1051c250bc487799a175edf57b5c..0f7086fab2e7bd7c06180700f529e7f37e356003 100644
--- a/codes/model.py
+++ b/codes/model.py
@@ -1,10 +1,10 @@
 # %%
-from codes.tb.tb import addTb
-from codes.tb.transforms import tb2khamvector, ifftn2tb
+from codes.tb.tb import add_tb
+from codes.tb.transforms import tb_to_khamvector, ifftn_to_tb
 from codes.mf import (
-    densityMatrix,
-    fermiOnGrid,
-    meanField,
+    density_matrix,
+    fermi_on_grid,
+    meanfield,
 )
 import numpy as np
 from scipy.fftpack import ifftn
@@ -16,10 +16,10 @@ class Model:
         self.h_int = h_int
         self.filling = filling
 
-        _firstKey = list(h_0)[0]
-        self._ndim = len(_firstKey)
-        self._size = h_0[_firstKey].shape[0]
-        self._localKey = tuple(np.zeros((self._ndim,), dtype=int))
+        _first_key = list(h_0)[0]
+        self._ndim = len(_first_key)
+        self._size = h_0[_first_key].shape[0]
+        self._local_key = tuple(np.zeros((self._ndim,), dtype=int))
 
         def _check_hermiticity(h):
             # assert hermiticity of the Hamiltonian
@@ -32,19 +32,19 @@ class Model:
         _check_hermiticity(h_0)
         _check_hermiticity(h_int)
 
-    def calculateEF(self):
-        self.EF = fermiOnGrid(self.kham, self.filling)
+    def calculate_EF(self):
+        self.EF = fermi_on_grid(self.kham, self.filling)
 
-    def makeDensityMatrixTb(self, mf_model, nK=200):
-        self.kham = tb2khamvector(addTb(self.h_0, mf_model), nK=nK, ndim=self._ndim)
-        self.calculateEF()
-        return ifftn2tb(
-            ifftn(densityMatrix(self.kham, self.EF), axes=np.arange(self._ndim))
+    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)
+        self.calculate_EF()
+        return ifftn_to_tb(
+            ifftn(density_matrix(self.kham, self.EF), axes=np.arange(self._ndim))
         )
 
-    def mfield(self, mf_model, nK=200):
-        densityMatrixTb = self.makeDensityMatrixTb(mf_model, nK=nK)
-        return addTb(
-            meanField(densityMatrixTb, self.h_int, n=self._ndim),
-            {self._localKey: -self.EF * np.eye(self._size)},
+    def mfield(self, mf_model, nk=200):
+        density_matrix_tb = self.make_density_matrix_tb(mf_model, nk=nk)
+        return add_tb(
+            meanfield(density_matrix_tb, self.h_int, n=self._ndim),
+            {self._local_key: -self.EF * np.eye(self._size)},
         )
diff --git a/codes/observables.py b/codes/observables.py
index ace8f03f25bae0bc940dae1254545208b23a9c08..52621e0e4e177995a205f38a1d68f9da9e749c40 100644
--- a/codes/observables.py
+++ b/codes/observables.py
@@ -1,12 +1,13 @@
 import numpy as np
 
-def expectationValue(densityMatrix, observable):
+
+def expectation_value(density_matrix, observable):
     """
     Compute the expectation value of an observable with respect to a density matrix.
 
     Parameters
     ----------
-    densityMatrix : dict
+    density_matrix : dict
         Density matrix in tight-binding format.
     observable : dict
         Observable in tight-binding format.
@@ -16,4 +17,9 @@ def expectationValue(densityMatrix, observable):
     float
         Expectation value.
     """
-    return np.sum([np.trace(densityMatrix[k] @ observable[k]) for k in frozenset(densityMatrix) & frozenset(observable)])
\ No newline at end of file
+    return np.sum(
+        [
+            np.trace(density_matrix[k] @ observable[k])
+            for k in frozenset(density_matrix) & frozenset(observable)
+        ]
+    )
diff --git a/codes/params/rparams.py b/codes/params/rparams.py
index 1fe03f97b9c2877b04aaa10e99301d80745afb5b..74c50887209d1f5514cd966681c68dc586bfa857 100644
--- a/codes/params/rparams.py
+++ b/codes/params/rparams.py
@@ -1,7 +1,12 @@
-from codes.params.matrixShaping import (complex_to_real, hop_dict_to_flat,
-                             real_to_complex, flat_to_hop_dict)
+from codes.params.matrixShaping import (
+    complex_to_real,
+    hop_dict_to_flat,
+    real_to_complex,
+    flat_to_hop_dict,
+)
 
-def mf2rParams(mf_model):
+
+def mf_to_rparams(mf_model):
     """
     Convert a mean-field tight-binding model to a set of real parameters.
 
@@ -15,9 +20,10 @@ def mf2rParams(mf_model):
     dict
         Real parameters.
     """
-    return complex_to_real(hop_dict_to_flat(mf_model)) # placeholder for now
+    return complex_to_real(hop_dict_to_flat(mf_model))  # placeholder for now
+
 
-def rParams2mf(rParams, keyList, size):
+def rparams_to_mf(rParams, key_list, size):
     """
     Extract mean-field tight-binding model from a set of real parameters.
 
@@ -34,5 +40,5 @@ def rParams2mf(rParams, keyList, size):
         Mean-field tight-binding model.
     """
 
-    flatMatrix = real_to_complex(rParams)
-    return flat_to_hop_dict(flatMatrix, (len(keyList), size, size), keyList)
\ No newline at end of file
+    flat_matrix = real_to_complex(rParams)
+    return flat_to_hop_dict(flat_matrix, (len(key_list), size, size), key_list)
diff --git a/codes/params/test_params.py b/codes/params/test_params.py
index d765bcf13c09aa8e830a545eb49a121138cfc614..66a31da1e03043cc7a2782fb92a7b63112c1ca56 100644
--- a/codes/params/test_params.py
+++ b/codes/params/test_params.py
@@ -1,17 +1,19 @@
 # %%
-from codes.params.rparams import mf2rParams, rParams2mf
+from codes.params.rparams import mf_to_rparams, rparams_to_mf
 from codes.kwant_helper.utils import generate_guess
-from codes.tb.tb import compareDicts
+from codes.tb.tb import compare_dicts
 import pytest
-repeatNumber = 10
+
+repeat_number = 10
 
 # %%
 ndof = 10
 vectors = ((0, 0), (1, 0), (-1, 0), (0, 1), (0, -1), (1, -1), (-1, 1), (1, 1), (-1, -1))
 
-@pytest.mark.repeat(repeatNumber)
+
+@pytest.mark.repeat(repeat_number)
 def test_parametrisation():
     mf_guess = generate_guess(vectors, ndof)
-    mf_params = mf2rParams(mf_guess)
-    mf_new = rParams2mf(mf_params, vectors, ndof)
-    compareDicts(mf_guess, mf_new)
\ No newline at end of file
+    mf_params = mf_to_rparams(mf_guess)
+    mf_new = rparams_to_mf(mf_params, vectors, ndof)
+    compare_dicts(mf_guess, mf_new)
diff --git a/codes/solvers.py b/codes/solvers.py
index c0acf567eaff6e0b775ef06a2ba5fd0d075da65a..ce83449348c37550ef3c93d57b8b1f2d96f4bc0d 100644
--- a/codes/solvers.py
+++ b/codes/solvers.py
@@ -1,11 +1,11 @@
-from codes.params.rparams import mf2rParams, rParams2mf
+from codes.params.rparams import mf_to_rparams, rparams_to_mf
 import scipy
 from functools import partial
-from codes.tb.tb import addTb
+from codes.tb.tb import add_tb
 import numpy as np
 
 
-def cost(mf_param, Model, nK=100):
+def cost(mf_param, Model, nk=100):
     """
     Define the cost function for fixed point iteration.
     The cost function is the difference between the input mean-field real space parametrisation
@@ -17,18 +17,18 @@ def cost(mf_param, Model, nK=100):
         The mean-field real space parametrisation.
     Model : Model
         The model object.
-    nK : int, optional
+    nk : int, optional
         The number of k-points to use in the grid. The default is 100.
     """
     shape = Model._size
-    mf_tb = rParams2mf(mf_param, list(Model.h_int), shape)
-    mf_tb_new = Model.mfield(mf_tb, nK=nK)
-    mf_params_new = mf2rParams(mf_tb_new)
+    mf_tb = rparams_to_mf(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)
     return mf_params_new - mf_param
 
 
 def solver(
-    Model, mf_guess, nK=100, optimizer=scipy.optimize.anderson, optimizer_kwargs={}
+    Model, mf_guess, nk=100, optimizer=scipy.optimize.anderson, optimizer_kwargs={}
 ):
     """
     Solve the mean-field self-consistent equation.
@@ -39,7 +39,7 @@ def solver(
         The model object.
     mf_guess : numpy.array
         The initial guess for the mean-field tight-binding model.
-    nK : int, optional
+    nk : int, optional
         The number of k-points to use in the grid. The default is 100.
     optimizer : scipy.optimize, optional
         The optimizer to use to solve for fixed-points. The default is scipy.optimize.anderson.
@@ -53,11 +53,11 @@ def solver(
     """
 
     shape = Model._size
-    mf_params = mf2rParams(mf_guess)
-    f = partial(cost, Model=Model, nK=nK)
-    result = rParams2mf(
+    mf_params = mf_to_rparams(mf_guess)
+    f = partial(cost, Model=Model, nk=nk)
+    result = rparams_to_mf(
         optimizer(f, mf_params, **optimizer_kwargs), list(Model.h_int), shape
     )
-    Model.calculateEF()
-    localKey = tuple(np.zeros((Model._ndim,), dtype=int))
-    return addTb(result, {localKey: -Model.EF * np.eye(shape)})
+    Model.calculate_EF()
+    local_key = tuple(np.zeros((Model._ndim,), dtype=int))
+    return add_tb(result, {local_key: -Model.EF * np.eye(shape)})
diff --git a/codes/tb/tb.py b/codes/tb/tb.py
index b744ebcc674b23d7c940fcaf5e11027cbec55d99..c943d405a95408353f8f392d5c18c5e758eed4a8 100644
--- a/codes/tb/tb.py
+++ b/codes/tb/tb.py
@@ -1,16 +1,17 @@
 import numpy as np
 
-def addTb(tb1, tb2):
+
+def add_tb(tb1, tb2):
     """
     Add up two tight-binding models together.
-    
+
     Parameters:
     -----------
     tb1 : dict
         Tight-binding model.
     tb2 : dict
         Tight-binding model.
-    
+
     Returns:
     --------
     dict
@@ -18,7 +19,7 @@ def addTb(tb1, tb2):
     """
     return {k: tb1.get(k, 0) + tb2.get(k, 0) for k in frozenset(tb1) | frozenset(tb2)}
 
-def compareDicts(dict1, dict2):
+
+def compare_dicts(dict1, dict2):
     for key in dict1.keys():
         assert np.allclose(dict1[key], dict2[key])
-
diff --git a/codes/tb/test_tb.py b/codes/tb/test_tb.py
index a4d442bae6152d2d2b661e399f6e7b670b54106c..6fee47cee4af4e21522cf59b5f02b210c287d16a 100644
--- a/codes/tb/test_tb.py
+++ b/codes/tb/test_tb.py
@@ -1,31 +1,31 @@
 # %%
 import numpy as np
-from codes.tb.tb import compareDicts
+from codes.tb.tb import compare_dicts
 from codes.kwant_helper import utils
 import itertools as it
-from codes.tb.transforms import kfunc2tb, tb2kfunc, tb2kham, tb2khamvector
+from codes.tb.transforms import kfunc_to_tb, tb_to_kfunc, tb_to_kham, tb_to_khamvector
 import pytest
 
-repeatNumber = 10
+repeat_number = 10
 
 # %%
 ndim = 2
-maxOrder = 5
-matrixSize = 5
-nK = 10
+max_order = 5
+matrix_size = 5
+nk = 10
 
 
-@pytest.mark.repeat(repeatNumber)
+@pytest.mark.repeat(repeat_number)
 def test_fourier():
-    keys = [np.arange(-maxOrder + 1, maxOrder) for i in range(ndim)]
+    keys = [np.arange(-max_order + 1, max_order) for i in range(ndim)]
     keys = it.product(*keys)
-    h_0 = {key: (np.random.rand(matrixSize, matrixSize) - 0.5) * 2 for key in keys}
-    kfunc = tb2kfunc(h_0)
-    tb_new = kfunc2tb(kfunc, nK, ndim=ndim)
-    compareDicts(h_0, tb_new)
+    h_0 = {key: (np.random.rand(matrix_size, matrix_size) - 0.5) * 2 for key in keys}
+    kfunc = tb_to_kfunc(h_0)
+    tb_new = kfunc_to_tb(kfunc, nk, ndim=ndim)
+    compare_dicts(h_0, tb_new)
 
 
-@pytest.mark.repeat(repeatNumber)
+@pytest.mark.repeat(repeat_number)
 def test_tbkham_transform():
     vectors = (
         (0, 0),
@@ -41,4 +41,6 @@ def test_tbkham_transform():
     ndof = 10
     h_0 = utils.generate_guess(vectors, ndof)
 
-    assert np.allclose(tb2kham(h_0, nK=nK, ndim=2), tb2khamvector(h_0, nK=nK, ndim=2))
+    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/tb/transforms.py b/codes/tb/transforms.py
index dbd64db699af5d78adf4804e810c01bbfc9a5594..6df1adad1dc85ec2220fb9a4748e8d1ea378290d 100644
--- a/codes/tb/transforms.py
+++ b/codes/tb/transforms.py
@@ -3,7 +3,7 @@ from scipy.fftpack import ifftn
 import itertools as it
 
 
-def tb2khamvector(h_0, nK, ndim):
+def tb_to_khamvector(h_0, nk, ndim):
     """
     Real-space tight-binding model to hamiltonian on k-space grid.
 
@@ -11,7 +11,7 @@ def tb2khamvector(h_0, nK, ndim):
     ----------
     h_0 : dict
         A dictionary with real-space vectors as keys and complex np.arrays as values.
-    nK : int
+    nk : int
         Number of k-points along each direction.
     ndim : int
         Number of dimensions.
@@ -23,8 +23,8 @@ def tb2khamvector(h_0, nK, ndim):
 
     """
 
-    ks = np.linspace(-np.pi, np.pi, nK, endpoint=False)
-    ks = np.concatenate((ks[nK // 2 :], ks[: nK // 2]), axis=0)  # shift for ifft
+    ks = np.linspace(-np.pi, np.pi, nk, endpoint=False)
+    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()))
@@ -40,7 +40,7 @@ def tb2khamvector(h_0, nK, ndim):
     return np.sum(tb_array * k_dependency, axis=0)
 
 
-def tb2kfunc(h_0):
+def tb_to_kfunc(h_0):
     """
     Fourier transforms a real-space tight-binding model to a k-space function.
 
@@ -64,13 +64,13 @@ def tb2kfunc(h_0):
     return bloch_ham
 
 
-def ifftn2tb(ifftArray):
+def ifftn_to_tb(ifft_array):
     """
     Converts an array from ifftn to a tight-binding model format.
 
     Parameters
     ----------
-    ifftArray : ndarray
+    ifft_array : ndarray
         An array obtained from ifftn.
 
     Returns
@@ -79,14 +79,14 @@ def ifftn2tb(ifftArray):
         A dictionary with real-space vectors as keys and complex np.arrays as values.
     """
 
-    size = ifftArray.shape[:-2]
+    size = ifft_array.shape[:-2]
 
     keys = [np.arange(-size[0] // 2 + 1, size[0] // 2) for i in range(len(size))]
     keys = it.product(*keys)
-    return {tuple(k): ifftArray[tuple(k)] for k in keys}
+    return {tuple(k): ifft_array[tuple(k)] for k in keys}
 
 
-def kfunc2kham(nk, hk, dim, return_ks=False, hermitian=True):
+def kfunc_to_kham(nk, hk, dim, return_ks=False, hermitian=True):
     """
     Evaluates Hamiltonian on a k-point grid.
 
@@ -132,13 +132,13 @@ def kfunc2kham(nk, hk, dim, return_ks=False, hermitian=True):
         return ham.reshape(*shape)
 
 
-def tb2kham(h_0, nK=200, ndim=1):
-    kfunc = tb2kfunc(h_0)
-    kham = kfunc2kham(nK, kfunc, ndim)
+def tb_to_kham(h_0, nk=200, ndim=1):
+    kfunc = tb_to_kfunc(h_0)
+    kham = kfunc_to_kham(nk, kfunc, ndim)
     return kham
 
 
-def kfunc2tb(kfunc, nSamples, ndim=1):
+def kfunc_to_tb(kfunc, n_samples, ndim=1):
     """
     Applies FFT on a k-space function to obtain a real-space components.
 
@@ -146,7 +146,7 @@ def kfunc2tb(kfunc, nSamples, ndim=1):
     ----------
     kfunc : function
         A function that takes a k-space vector and returns a np.array.
-    nSamples : int
+    n_samples : int
         Number of samples to take in k-space.
 
     Returns
@@ -157,23 +157,23 @@ def kfunc2tb(kfunc, nSamples, ndim=1):
     """
 
     ks = np.linspace(
-        -np.pi, np.pi, nSamples, endpoint=False
-    )  # for now nSamples need to be even such that 0 is in the middle
+        -np.pi, np.pi, n_samples, endpoint=False
+    )  # for now n_samples need to be even such that 0 is in the middle
     ks = np.concatenate(
-        (ks[nSamples // 2 :], ks[: nSamples // 2]), axis=0
+        (ks[n_samples // 2 :], ks[: n_samples // 2]), axis=0
     )  # shift for ifft
 
     if ndim == 1:
-        kfuncOnGrid = np.array([kfunc(k) for k in ks])
+        kfunc_on_grid = np.array([kfunc(k) for k in ks])
     if ndim == 2:
-        kfuncOnGrid = np.array([[kfunc((kx, ky)) for ky in ks] for kx in ks])
+        kfunc_on_grid = np.array([[kfunc((kx, ky)) for ky in ks] for kx in ks])
     if ndim > 2:
         raise NotImplementedError("n > 2 not implemented")
-    ifftnArray = ifftn(kfuncOnGrid, axes=np.arange(ndim))
-    return ifftn2tb(ifftnArray)
+    ifftn_array = ifftn(kfunc_on_grid, axes=np.arange(ndim))
+    return ifftn_to_tb(ifftn_array)
 
 
-def hk2h_0(hk, hopping_vecs, ks=None):
+def hk_to_h0(hk, hopping_vecs, ks=None):
     """
     Extract hopping matrices from Bloch Hamiltonian.
 
diff --git a/codes/tb/utils.py b/codes/tb/utils.py
index 065690645df6e3b9f26c3fbffc7c59b4667bdc64..1d00831c97c7cfcaba577e6088144037f97a3b9d 100644
--- a/codes/tb/utils.py
+++ b/codes/tb/utils.py
@@ -1,32 +1,33 @@
 import numpy as np
-from codes.tb.transforms import tb2kfunc
+from codes.tb.transforms import tb_to_kfunc
 import itertools as it
 
-def compute_gap(h, E_F=0, n=100):
+
+def compute_gap(h, fermi_energy=0, n=100):
     """
-    Compute gap.
+     Compute gap.
 
-    Parameters:
-    -----------
-    h : dict
-    Tight-binding model for which to compute the gap.
-    E_F : float
-    Fermi energy.
-    n : int
-    Number of k-points to sample along each dimension.
+     Parameters:
+     -----------
+     h : dict
+     Tight-binding model for which to compute the gap.
+    fermi_energy : float
+     Fermi energy.
+     n : int
+     Number of k-points to sample along each dimension.
 
-    Returns:
-    --------
-    gap : float
-    Indirect gap.
+     Returns:
+     --------
+     gap : float
+     Indirect gap.
     """
     ndim = len(list(h)[0])
-    hkfunc = tb2kfunc(h)
-    kArray = np.linspace(0, 2*np.pi, n)
-    kGrid = list(it.product(*[kArray for i in range(ndim)]))
-    kham = np.array([hkfunc(k) for k in kGrid]) 
+    hkfunc = tb_to_kfunc(h)
+    k_array = np.linspace(0, 2 * np.pi, n)
+    kgrid = list(it.product(*[k_array for i in range(ndim)]))
+    kham = np.array([hkfunc(k) for k in kgrid])
     vals = np.linalg.eigvalsh(kham)
 
-    emax = np.max(vals[vals <= E_F])
-    emin = np.min(vals[vals > E_F])
-    return np.abs(emin - emax)
\ No newline at end of file
+    emax = np.max(vals[vals <= fermi_energy])
+    emin = np.min(vals[vals > fermi_energy])
+    return np.abs(emin - emax)
diff --git a/codes/test_graphene.py b/codes/test_graphene.py
index 59769cc10a4159c28fe9ec9663eb44be0aabcb6d..1ebcd41ae1e313e7b120713756ef69e91380cb35 100644
--- a/codes/test_graphene.py
+++ b/codes/test_graphene.py
@@ -5,10 +5,10 @@ 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.tb import addTb
+from codes.tb.tb import add_tb
 import pytest
 
-repeatNumber = 10
+repeat_number = 10
 # %%
 graphene_builder, int_builder = kwant_examples.graphene_extended_hubbard()
 h_0 = utils.builder2tb(graphene_builder)
@@ -31,39 +31,39 @@ def gap_prediction(U, V):
     params = {"U": U, "V": V}
 
     # Compare to phase diagram in https://arxiv.org/pdf/1204.4531.pdf
-    upperPhaseLine = 0.181 * U + 0.416
-    lowerPhaseLine = 1.707 * U - 3.823
-    triplePoint = (2.78, 0.92)
+    upper_phase_line = 0.181 * U + 0.416
+    lower_phase_line = 1.707 * U - 3.823
+    triple_point = (2.78, 0.92)
 
     gapped = False
-    if triplePoint < (U, V):
+    if triple_point < (U, V):
         gapped = True
-    elif (upperPhaseLine < V) | (lowerPhaseLine > V):
+    elif (upper_phase_line < V) | (lower_phase_line > V):
         gapped = True
 
     # the mean-field calculation
     filling = 2
-    nK = 20
+    nk = 20
 
     h_int = utils.builder2tb(int_builder, params)
     guess = utils.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})
-    gap = compute_gap(addTb(h_0, mf_sol), n=100)
+    mf_sol = solver(model, guess, nk=nk, optimizer_kwargs={"verbose": True, "M": 0})
+    gap = compute_gap(add_tb(h_0, mf_sol), n=100)
 
     # Check if the gap is predicted correctly
     if gap > 0.1:
-        gappedPredicted = True
+        gapped_predicted = True
     else:
-        gappedPredicted = False
+        gapped_predicted = False
     assert (
-        gapped == gappedPredicted
+        gapped == gapped_predicted
     ), f"Mean-field theory failed to predict the gap for U = {U}, V = {V}"
 
 
 # %%
-@pytest.mark.repeat(repeatNumber)
+@pytest.mark.repeat(repeat_number)
 def test_gap():
     U = np.random.uniform(0, 2)
     V = np.random.uniform(0, 1)
diff --git a/examples/tEnergyTest.ipynb b/examples/tEnergyTest.ipynb
index 6cf315bfa09455995477d2483e33ffbf3f20da15..80a5a8b637792a64cf266cee32cabec2e2ba8e79 100644
--- a/examples/tEnergyTest.ipynb
+++ b/examples/tEnergyTest.ipynb
@@ -64,11 +64,11 @@
     "    )\n",
     "    vals, vecs = np.linalg.eigh(tb_mf_k)\n",
     "    EF = utils.get_fermi_energy(vals, filling)\n",
-    "    densityMatrix = hf.density_matrix(vals, vecs, EF)\n",
+    "    density_matrix = hf.density_matrix(vals, vecs, EF)\n",
     "\n",
-    "    return tb_mf_k, densityMatrix, _model.H_int\n",
+    "    return tb_mf_k, density_matrix, _model.H_int\n",
     "\n",
-    "tb_mf0, densityMatrix0, H_int0 = groundstate(U0)\n",
+    "tb_mf0, density_matrix0, H_int0 = groundstate(U0)\n",
     "mf0 = tb_mf0 - hamiltonians_0\n",
     "\n",
     "\n",
@@ -77,10 +77,10 @@
     "    hamiltonian = hamiltonians_0 + mf0 * alpha\n",
     "    vals, vecs = np.linalg.eigh(hamiltonian)\n",
     "    EF = utils.get_fermi_energy(vals, filling)\n",
-    "    densityMatrix = hf.density_matrix(vals, vecs, EF)\n",
+    "    density_matrix = hf.density_matrix(vals, vecs, EF)\n",
     "    return hf.total_energy(\n",
     "        hamiltonians_0 + np.sign(alpha) * mf0,\n",
-    "        densityMatrix,\n",
+    "        density_matrix,\n",
     "    )"
    ]
   },
diff --git a/examples/testInterface.ipynb b/examples/testInterface.ipynb
index 0bec84ae85628698bc0d6718df44e337c1dfafb7..c3bf96057e4027d357a8979d7ba0868f87b73b61 100644
--- a/examples/testInterface.ipynb
+++ b/examples/testInterface.ipynb
@@ -49,7 +49,7 @@
       "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
       "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
       "Cell \u001b[0;32mIn[16], line 10\u001b[0m\n\u001b[1;32m      8\u001b[0m mf_model \u001b[38;5;241m=\u001b[39m {(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m0\u001b[39m) : guess}\n\u001b[1;32m      9\u001b[0m tb \u001b[38;5;241m=\u001b[39m Model(h_0, h_int, filling)\n\u001b[0;32m---> 10\u001b[0m mf_sol \u001b[38;5;241m=\u001b[39m \u001b[43msolver\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmf_model\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnK\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m20\u001b[39;49m\u001b[43m)\u001b[49m \u001b[38;5;66;03m# does not converge as fast, something is wrong in 2D\u001b[39;00m\n",
-      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/solvers.py:59\u001b[0m, in \u001b[0;36msolver\u001b[0;34m(Model, mf_guess, nK, optimizer, optimizer_kwargs)\u001b[0m\n\u001b[1;32m     56\u001b[0m mf_params \u001b[38;5;241m=\u001b[39m mf2rParams(mf_guess)\n\u001b[1;32m     57\u001b[0m f \u001b[38;5;241m=\u001b[39m partial(cost, Model\u001b[38;5;241m=\u001b[39mModel, nK\u001b[38;5;241m=\u001b[39mnK)\n\u001b[1;32m     58\u001b[0m result \u001b[38;5;241m=\u001b[39m rParams2mf(\n\u001b[0;32m---> 59\u001b[0m     \u001b[43moptimizer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmf_params\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43moptimizer_kwargs\u001b[49m\u001b[43m)\u001b[49m, \u001b[38;5;28mlist\u001b[39m(Model\u001b[38;5;241m.\u001b[39mh_int), shape\n\u001b[1;32m     60\u001b[0m )\n\u001b[1;32m     61\u001b[0m Model\u001b[38;5;241m.\u001b[39mcalculateEF(nK\u001b[38;5;241m=\u001b[39mnK)\n\u001b[1;32m     62\u001b[0m localKey \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mtuple\u001b[39m(np\u001b[38;5;241m.\u001b[39mzeros((Model\u001b[38;5;241m.\u001b[39m_ndim,), dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mint\u001b[39m))\n",
+      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/solvers.py:59\u001b[0m, in \u001b[0;36msolver\u001b[0;34m(Model, mf_guess, nK, optimizer, optimizer_kwargs)\u001b[0m\n\u001b[1;32m     56\u001b[0m mf_params \u001b[38;5;241m=\u001b[39m mf2rParams(mf_guess)\n\u001b[1;32m     57\u001b[0m f \u001b[38;5;241m=\u001b[39m partial(cost, Model\u001b[38;5;241m=\u001b[39mModel, nK\u001b[38;5;241m=\u001b[39mnK)\n\u001b[1;32m     58\u001b[0m result \u001b[38;5;241m=\u001b[39m rParams2mf(\n\u001b[0;32m---> 59\u001b[0m     \u001b[43moptimizer\u001b[49m\u001b[43m(\u001b[49m\u001b[43mf\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmf_params\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43moptimizer_kwargs\u001b[49m\u001b[43m)\u001b[49m, \u001b[38;5;28mlist\u001b[39m(Model\u001b[38;5;241m.\u001b[39mh_int), shape\n\u001b[1;32m     60\u001b[0m )\n\u001b[1;32m     61\u001b[0m Model\u001b[38;5;241m.\u001b[39mcalculateEF(nK\u001b[38;5;241m=\u001b[39mnK)\n\u001b[1;32m     62\u001b[0m local_key \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mtuple\u001b[39m(np\u001b[38;5;241m.\u001b[39mzeros((Model\u001b[38;5;241m.\u001b[39m_ndim,), dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mint\u001b[39m))\n",
       "File \u001b[0;32m<string>:6\u001b[0m, in \u001b[0;36manderson\u001b[0;34m(F, xin, iter, alpha, w0, M, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, **kw)\u001b[0m\n",
       "File \u001b[0;32m~/mambaforge/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:214\u001b[0m, in \u001b[0;36mnonlin_solve\u001b[0;34m(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception)\u001b[0m\n\u001b[1;32m    212\u001b[0m \u001b[38;5;66;03m# Line search, or Newton step\u001b[39;00m\n\u001b[1;32m    213\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m line_search:\n\u001b[0;32m--> 214\u001b[0m     s, x, Fx, Fx_norm_new \u001b[38;5;241m=\u001b[39m \u001b[43m_nonlin_line_search\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mFx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdx\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    215\u001b[0m \u001b[43m                                                \u001b[49m\u001b[43mline_search\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    216\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    217\u001b[0m     s \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1.0\u001b[39m\n",
       "File \u001b[0;32m~/mambaforge/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:293\u001b[0m, in \u001b[0;36m_nonlin_line_search\u001b[0;34m(func, x, Fx, dx, search_type, rdiff, smin)\u001b[0m\n\u001b[1;32m    290\u001b[0m     s, phi1, phi0 \u001b[38;5;241m=\u001b[39m scalar_search_wolfe1(phi, derphi, tmp_phi[\u001b[38;5;241m0\u001b[39m],\n\u001b[1;32m    291\u001b[0m                                          xtol\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1e-2\u001b[39m, amin\u001b[38;5;241m=\u001b[39msmin)\n\u001b[1;32m    292\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m search_type \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124marmijo\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[0;32m--> 293\u001b[0m     s, phi1 \u001b[38;5;241m=\u001b[39m \u001b[43mscalar_search_armijo\u001b[49m\u001b[43m(\u001b[49m\u001b[43mphi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtmp_phi\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43mtmp_phi\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    294\u001b[0m \u001b[43m                                   \u001b[49m\u001b[43mamin\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msmin\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    296\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m s \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    297\u001b[0m     \u001b[38;5;66;03m# XXX: No suitable step length found. Take the full Newton step,\u001b[39;00m\n\u001b[1;32m    298\u001b[0m     \u001b[38;5;66;03m#      and hope for the best.\u001b[39;00m\n\u001b[1;32m    299\u001b[0m     s \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1.0\u001b[39m\n",
@@ -57,10 +57,10 @@
       "File \u001b[0;32m~/mambaforge/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:277\u001b[0m, in \u001b[0;36m_nonlin_line_search.<locals>.phi\u001b[0;34m(s, store)\u001b[0m\n\u001b[1;32m    275\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m tmp_phi[\u001b[38;5;241m0\u001b[39m]\n\u001b[1;32m    276\u001b[0m xt \u001b[38;5;241m=\u001b[39m x \u001b[38;5;241m+\u001b[39m s\u001b[38;5;241m*\u001b[39mdx\n\u001b[0;32m--> 277\u001b[0m v \u001b[38;5;241m=\u001b[39m \u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43mxt\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    278\u001b[0m p \u001b[38;5;241m=\u001b[39m _safe_norm(v)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m\n\u001b[1;32m    279\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m store:\n",
       "File \u001b[0;32m~/mambaforge/lib/python3.11/site-packages/scipy/optimize/_nonlin.py:168\u001b[0m, in \u001b[0;36mnonlin_solve.<locals>.func\u001b[0;34m(z)\u001b[0m\n\u001b[1;32m    167\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfunc\u001b[39m(z):\n\u001b[0;32m--> 168\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m _as_inexact(\u001b[43mF\u001b[49m\u001b[43m(\u001b[49m\u001b[43m_array_like\u001b[49m\u001b[43m(\u001b[49m\u001b[43mz\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mx0\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m)\u001b[38;5;241m.\u001b[39mflatten()\n",
       "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/solvers.py:25\u001b[0m, in \u001b[0;36mcost\u001b[0;34m(mf_param, Model, nK)\u001b[0m\n\u001b[1;32m     23\u001b[0m shape \u001b[38;5;241m=\u001b[39m Model\u001b[38;5;241m.\u001b[39m_size\n\u001b[1;32m     24\u001b[0m mf_tb \u001b[38;5;241m=\u001b[39m rParams2mf(mf_param, \u001b[38;5;28mlist\u001b[39m(Model\u001b[38;5;241m.\u001b[39mh_int), shape)\n\u001b[0;32m---> 25\u001b[0m mf_tb_new \u001b[38;5;241m=\u001b[39m \u001b[43mModel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmfieldFFT\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmf_tb\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnK\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnK\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     26\u001b[0m mf_params_new \u001b[38;5;241m=\u001b[39m mf2rParams(mf_tb_new)\n\u001b[1;32m     27\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m mf_params_new \u001b[38;5;241m-\u001b[39m mf_param\n",
-      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/model.py:34\u001b[0m, in \u001b[0;36mModel.mfieldFFT\u001b[0;34m(self, mf_model, nK)\u001b[0m\n\u001b[1;32m     33\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmfieldFFT\u001b[39m(\u001b[38;5;28mself\u001b[39m, mf_model, nK\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m200\u001b[39m):\n\u001b[0;32m---> 34\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdensityMatrix \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmakeDensityMatrix\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmf_model\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     35\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m addTb(\n\u001b[1;32m     36\u001b[0m         meanFieldFFT(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdensityMatrix, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mh_int, n\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ndim, nK\u001b[38;5;241m=\u001b[39mnK),\n\u001b[1;32m     37\u001b[0m         {\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_localKey: \u001b[38;5;241m-\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mEF \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39meye(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_size)},\n\u001b[1;32m     38\u001b[0m     )\n",
-      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/model.py:20\u001b[0m, in \u001b[0;36mModel.makeDensityMatrix\u001b[0;34m(self, mf_model, nK)\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmakeDensityMatrix\u001b[39m(\u001b[38;5;28mself\u001b[39m, mf_model, nK\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m200\u001b[39m):\n\u001b[1;32m     19\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhkfunc \u001b[38;5;241m=\u001b[39m tb2kfunc(addTb(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mh_0, mf_model))\n\u001b[0;32m---> 20\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcalculateEF\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnK\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnK\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     21\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m densityMatrixGenerator(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhkfunc, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mEF)\n",
-      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/model.py:24\u001b[0m, in \u001b[0;36mModel.calculateEF\u001b[0;34m(self, nK)\u001b[0m\n\u001b[1;32m     23\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mcalculateEF\u001b[39m(\u001b[38;5;28mself\u001b[39m, nK\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m200\u001b[39m):\n\u001b[0;32m---> 24\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mEF \u001b[38;5;241m=\u001b[39m \u001b[43mfermiOnGrid\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mhkfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfilling\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnK\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnK\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mndim\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_ndim\u001b[49m\u001b[43m)\u001b[49m\n",
-      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/mf.py:54\u001b[0m, in \u001b[0;36mfermiOnGrid\u001b[0;34m(hkfunc, filling, nK, ndim)\u001b[0m\n\u001b[1;32m     52\u001b[0m     hkarray \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([hkfunc(k) \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m ks])\n\u001b[1;32m     53\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ndim \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[0;32m---> 54\u001b[0m     hkarray \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray(\u001b[43m[\u001b[49m\u001b[43m[\u001b[49m\u001b[43mhkfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mky\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mkx\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mks\u001b[49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mky\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mks\u001b[49m\u001b[43m]\u001b[49m)\n\u001b[1;32m     55\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m ndim \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[1;32m     56\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mNotImplementedError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFermi energy calculation is not implemented for ndim > 2\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
+      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/model.py:34\u001b[0m, in \u001b[0;36mModel.mfieldFFT\u001b[0;34m(self, mf_model, nK)\u001b[0m\n\u001b[1;32m     33\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmfieldFFT\u001b[39m(\u001b[38;5;28mself\u001b[39m, mf_model, nK\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m200\u001b[39m):\n\u001b[0;32m---> 34\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdensity_matrix \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmakedensity_matrix\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmf_model\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     35\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m addTb(\n\u001b[1;32m     36\u001b[0m         meanFieldFFT(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdensity_matrix, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mh_int, n\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_ndim, nK\u001b[38;5;241m=\u001b[39mnK),\n\u001b[1;32m     37\u001b[0m         {\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_local_key: \u001b[38;5;241m-\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mEF \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39meye(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_size)},\n\u001b[1;32m     38\u001b[0m     )\n",
+      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/model.py:20\u001b[0m, in \u001b[0;36mModel.makedensity_matrix\u001b[0;34m(self, mf_model, nK)\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmakedensity_matrix\u001b[39m(\u001b[38;5;28mself\u001b[39m, mf_model, nK\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m200\u001b[39m):\n\u001b[1;32m     19\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhkfunc \u001b[38;5;241m=\u001b[39m tb2kfunc(addTb(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mh_0, mf_model))\n\u001b[0;32m---> 20\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcalculateEF\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnK\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnK\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     21\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m density_matrixGenerator(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhkfunc, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mEF)\n",
+      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/model.py:24\u001b[0m, in \u001b[0;36mModel.calculateEF\u001b[0;34m(self, nK)\u001b[0m\n\u001b[1;32m     23\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mcalculateEF\u001b[39m(\u001b[38;5;28mself\u001b[39m, nK\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m200\u001b[39m):\n\u001b[0;32m---> 24\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mEF \u001b[38;5;241m=\u001b[39m \u001b[43mfermi_on_grid\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mhkfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfilling\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnK\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnK\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mndim\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_ndim\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/mf.py:54\u001b[0m, in \u001b[0;36mfermi_on_grid\u001b[0;34m(hkfunc, filling, nK, ndim)\u001b[0m\n\u001b[1;32m     52\u001b[0m     hkarray \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([hkfunc(k) \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m ks])\n\u001b[1;32m     53\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ndim \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[0;32m---> 54\u001b[0m     hkarray \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray(\u001b[43m[\u001b[49m\u001b[43m[\u001b[49m\u001b[43mhkfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mky\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mkx\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mks\u001b[49m\u001b[43m]\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mky\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mks\u001b[49m\u001b[43m]\u001b[49m)\n\u001b[1;32m     55\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m ndim \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[1;32m     56\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mNotImplementedError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFermi energy calculation is not implemented for ndim > 2\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
       "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/mf.py:54\u001b[0m, in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m     52\u001b[0m     hkarray \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([hkfunc(k) \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m ks])\n\u001b[1;32m     53\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ndim \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[0;32m---> 54\u001b[0m     hkarray \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([\u001b[43m[\u001b[49m\u001b[43mhkfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mky\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mfor\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mkx\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01min\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mks\u001b[49m\u001b[43m]\u001b[49m \u001b[38;5;28;01mfor\u001b[39;00m ky \u001b[38;5;129;01min\u001b[39;00m ks])\n\u001b[1;32m     55\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m ndim \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[1;32m     56\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mNotImplementedError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFermi energy calculation is not implemented for ndim > 2\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
       "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/mf.py:54\u001b[0m, in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m     52\u001b[0m     hkarray \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([hkfunc(k) \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m ks])\n\u001b[1;32m     53\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m ndim \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[0;32m---> 54\u001b[0m     hkarray \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([[\u001b[43mhkfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkx\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mky\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mfor\u001b[39;00m kx \u001b[38;5;129;01min\u001b[39;00m ks] \u001b[38;5;28;01mfor\u001b[39;00m ky \u001b[38;5;129;01min\u001b[39;00m ks])\n\u001b[1;32m     55\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m ndim \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[1;32m     56\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mNotImplementedError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mFermi energy calculation is not implemented for ndim > 2\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
       "File \u001b[0;32m~/Projects/kwant-scf/examples/codes/tb/transforms.py:24\u001b[0m, in \u001b[0;36mtb2kfunc.<locals>.bloch_ham\u001b[0;34m(k)\u001b[0m\n\u001b[1;32m     22\u001b[0m ham \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n\u001b[1;32m     23\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m vector \u001b[38;5;129;01min\u001b[39;00m h_0\u001b[38;5;241m.\u001b[39mkeys():\n\u001b[0;32m---> 24\u001b[0m     ham \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m h_0[vector] \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39mexp(\n\u001b[1;32m     25\u001b[0m         \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39mj \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39mdot(k, np\u001b[38;5;241m.\u001b[39marray(vector, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mfloat\u001b[39m))\n\u001b[1;32m     26\u001b[0m     )\n\u001b[1;32m     27\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ham\n",
diff --git a/profiling/graphene.py b/profiling/graphene.py
index a82e382721a2321ff6ba02bc51903c4adf1ccc09..9731d90fb0e984ab34250c976b3971819cbc77f8 100644
--- a/profiling/graphene.py
+++ b/profiling/graphene.py
@@ -12,7 +12,7 @@ graphene_builder, int_builder = kwant_examples.graphene_extended_hubbard()
 
 params = {"U": 0.5, "V": 1.1}
 filling = 2
-nK = 600
+nk = 600
 
 h_int = utils.builder2tb(int_builder, params)
 h_0 = utils.builder2tb(graphene_builder)
@@ -22,7 +22,7 @@ model = Model(h_0, h_int, filling)
 
 
 def scf_loop():
-    model.mfield(guess, nK=nK)
+    model.mfield(guess, nk=nk)
 
 
 # %% Memory profile
@@ -40,12 +40,12 @@ profiler.write_html(path="timeProfile.html")
 # %%
 number = 1
 
-timeSCF = timeit.timeit(scf_loop, number=number) / number
+time_scf = timeit.timeit(scf_loop, number=number) / number
 
-H = np.random.rand(nK, nK)
+H = np.random.rand(nk, nk)
 H += H.T.conj()
-timeDiag = timeit.timeit(lambda: np.linalg.eigh(H), number=number) / number
+time_diag = timeit.timeit(lambda: np.linalg.eigh(H), number=number) / number
 
 print(
-    f"Single SCF loop takes {timeSCF} whereas a single diagonalization of a corresponding system takes {timeDiag}"
+    f"Single SCF loop takes {time_scf} whereas a single diagonalization of a corresponding system takes {time_diag}"
 )