remove old functions and remove kvector ending on new functions

codes/mf.py

@@ -2,36 +2,6 @@ import numpy as np
 from codes.tb.tb import addTb
 
 
-def densityMatrixGenerator(hkfunc, E_F):
-    """
-    Generate a function that returns the density matrix at a given k-point.
-
-    Parameters
-    ----------
-    hkfunc : function
-        Function that return Hamiltonian at a given k-point.
-    E_F : float
-        Fermi level
-
-    Returns
-    -------
-    densityMatrixFunc : function
-        Returns a density matrix at a given k-point (kx, kx, ...)
-    """
-
-    def densityMatrixFunc(k):
-        hk = hkfunc(k)
-        vals, vecs = np.linalg.eigh(hk)
-        unocc_vals = vals > E_F
-        occ_vecs = vecs
-        occ_vecs[:, unocc_vals] = 0
-
-        # Outter products between eigenvectors
-        return occ_vecs @ occ_vecs.T.conj()
-
-    return densityMatrixFunc
-
-
 def densityMatrix(kham, E_F):
     """
     Parameters
@@ -56,24 +26,7 @@ def densityMatrix(kham, E_F):
     return densityMatrixKgrid
 
 
-def fermiOnGridkvector(kham, filling):
-
-    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 fermiOnGrid(hkfunc, filling, nK=100, ndim=1):  # need to extend to 2D
+def fermiOnGrid(kham, filling):
     """
     Compute the Fermi energy on a grid of k-points.
 
@@ -88,17 +41,8 @@ def fermiOnGrid(hkfunc, filling, nK=100, ndim=1):  # need to extend to 2D
     E_F : float
         Fermi energy
     """
-    ks = np.linspace(-np.pi, np.pi, nK, endpoint=False)
-    if ndim == 1:
-        hkarray = np.array([hkfunc(k) for k in ks])
-    if ndim == 2:
-        hkarray = np.array([[hkfunc((kx, ky)) for kx in ks] for ky in ks])
-    elif ndim > 2:
-        raise NotImplementedError(
-            "Fermi energy calculation is not implemented for ndim > 2"
-        )
 
-    vals = np.linalg.eigvalsh(hkarray)
+    vals = np.linalg.eigvalsh(kham)
 
     norbs = vals.shape[-1]
     vals_flat = np.sort(vals.flatten())
@@ -113,32 +57,7 @@ def fermiOnGrid(hkfunc, filling, nK=100, ndim=1):  # need to extend to 2D
         return fermi
 
 
-def meanFieldFFTkvector(densityMatrixTb, h_int, n=2):
-
-    localKey = tuple(np.zeros((n,), dtype=int))
-
-    direct = {
-        localKey: np.sum(
-            np.array(
-                [
-                    np.diag(
-                        np.einsum("pp,pn->n", densityMatrixTb[localKey], h_int[vec])
-                    )
-                    for vec in frozenset(h_int)
-                ]
-            ),
-            axis=0,
-        )
-    }
-
-    exchange = {
-        vec: -1 * h_int.get(vec, 0) * densityMatrixTb[vec]  # / (2 * np.pi)#**2
-        for vec in frozenset(h_int)
-    }
-    return addTb(direct, exchange)
-
-
-def meanField(densityMatrix, int_model, n=2, nK=100):
+def meanField(densityMatrixTb, h_int, n=2):
     """
     Compute the mean-field in k-space.
 
@@ -154,6 +73,7 @@ def meanField(densityMatrix, int_model, n=2, nK=100):
     dict
         Mean-field tb model.
     """
+
     localKey = tuple(np.zeros((n,), dtype=int))
 
     direct = {
@@ -161,9 +81,9 @@ def meanField(densityMatrix, int_model, n=2, nK=100):
             np.array(
                 [
                     np.diag(
-                        np.einsum("pp,pn->n", densityMatrix[localKey], int_model[vec])
+                        np.einsum("pp,pn->n", densityMatrixTb[localKey], h_int[vec])
                     )
-                    for vec in frozenset(int_model)
+                    for vec in frozenset(h_int)
                 ]
             ),
             axis=0,
@@ -171,7 +91,7 @@ def meanField(densityMatrix, int_model, n=2, nK=100):
     }
 
     exchange = {
-        vec: -1 * int_model.get(vec, 0) * densityMatrix[vec]  # / (2 * np.pi)#**2
-        for vec in frozenset(int_model)
+        vec: -1 * h_int.get(vec, 0) * densityMatrixTb[vec]  # / (2 * np.pi)#**2
+        for vec in frozenset(h_int)
     }
     return addTb(direct, exchange)

codes/model.py

 # %%
 from codes.tb.tb import addTb
-from codes.tb.transforms import tb2kfunc, tb2kham, kdens2tbFFT, kfunc2tb, ifftn2tb
+from codes.tb.transforms import tb2kham, kdens2tbFFT, ifftn2tb
 from codes.mf import (
-    densityMatrixGenerator,
     densityMatrix,
-    fermiOnGridkvector,
-    meanFieldFFTkvector,
-    meanField,
     fermiOnGrid,
+    meanField,
 )
 import numpy as np
 
@@ -34,37 +31,19 @@ class Model:
         _check_hermiticity(h_0)
         _check_hermiticity(h_int)
 
-    def calculateEF(self, nK=200):
-        self.EF = fermiOnGrid(self.hkfunc, self.filling, nK=nK, ndim=self._ndim)
+    def calculateEF(self):
+        self.EF = fermiOnGrid(self.kham, self.filling)
 
     def makeDensityMatrix(self, mf_model, nK=200):
-        self.hkfunc = tb2kfunc(addTb(self.h_0, mf_model))
-        self.calculateEF(nK=nK)
-        return kfunc2tb(
-            densityMatrixGenerator(self.hkfunc, self.EF), nSamples=nK, ndim=self._ndim
-        )
-
-    def mfield(self, mf_model, nK=200):
-        self.densityMatrix = self.makeDensityMatrix(mf_model, nK=nK)
-        return addTb(
-            meanField(self.densityMatrix, self.h_int, n=self._ndim, nK=nK),
-            {self._localKey: -self.EF * np.eye(self._size)},
-        )
-
-    #######################
-    def calculateEFkvector(self):
-        self.EF = fermiOnGridkvector(self.kham, self.filling)
-
-    def makeDensityMatrixkvector(self, mf_model, nK=200):
         self.kham = tb2kham(addTb(self.h_0, mf_model), nK=nK, ndim=self._ndim)
-        self.calculateEFkvector()
+        self.calculateEF()
         return densityMatrix(self.kham, self.EF)
 
-    def mfieldFFTkvector(self, mf_model, nK=200):
-        densityMatrix = self.makeDensityMatrixkvector(mf_model, nK=nK)
+    def mfield(self, mf_model, nK=200):
+        densityMatrix = self.makeDensityMatrix(mf_model, nK=nK)
         densityMatrixTb = ifftn2tb(kdens2tbFFT(densityMatrix, self._ndim))
         return addTb(
-            meanFieldFFTkvector(densityMatrixTb, self.h_int, n=self._ndim),
+            meanField(densityMatrixTb, self.h_int, n=self._ndim),
             {self._localKey: -self.EF * np.eye(self._size)},
         )
 

codes/solvers.py

@@ -27,28 +27,6 @@ def cost(mf_param, Model, nK=100):
     return mf_params_new - mf_param
 
 
-def costkvector(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
-    and a new mean-field.
-
-    Parameters
-    ----------
-    mf_param : numpy.array
-        The mean-field real space parametrisation.
-    Model : Model
-        The model object.
-    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.mfieldFFTkvector(mf_tb, nK=nK)
-    mf_params_new = mf2rParams(mf_tb_new)
-    return mf_params_new - mf_param
-
-
 def solver(
     Model, mf_guess, nK=100, optimizer=scipy.optimize.anderson, optimizer_kwargs={}
 ):
@@ -80,42 +58,6 @@ def solver(
     result = rParams2mf(
         optimizer(f, mf_params, **optimizer_kwargs), list(Model.h_int), shape
     )
-    Model.calculateEF(nK=nK)
-    localKey = tuple(np.zeros((Model._ndim,), dtype=int))
-    return addTb(result, {localKey: -Model.EF * np.eye(shape)})
-
-
-def solverkvector(
-    Model, mf_guess, nK=100, optimizer=scipy.optimize.anderson, optimizer_kwargs={}
-):
-    """
-    Solve the mean-field self-consistent equation.
-
-    Parameters
-    ----------
-    Model : Model
-        The model object.
-    mf_guess : numpy.array
-        The initial guess for the mean-field tight-binding model.
-    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.
-    optimizer_kwargs : dict, optional
-        The keyword arguments to pass to the optimizer. The default is {}.
-
-    Returns
-    -------
-    result : numpy.array
-        The mean-field tight-binding model.
-    """
-
-    shape = Model._size
-    mf_params = mf2rParams(mf_guess)
-    f = partial(costkvector, Model=Model, nK=nK)
-    result = rParams2mf(
-        optimizer(f, mf_params, **optimizer_kwargs), list(Model.h_int), shape
-    )
-    Model.calculateEFkvector()
+    Model.calculateEF()
     localKey = tuple(np.zeros((Model._ndim,), dtype=int))
     return addTb(result, {localKey: -Model.EF * np.eye(shape)})

codes/test_graphene.py

 # %%
 import numpy as np
 from codes.model import Model
-from codes.solvers import solverkvector
+from codes.solvers import solver
 from codes import kwant_examples
 from codes.kwant_helper import utils
 from codes.tb.utils import compute_gap
@@ -49,9 +49,7 @@ def gap_prediction(U, V):
     guess = utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
     model = Model(h_0, h_int, filling)
 
-    mf_sol = solverkvector(
-        model, guess, nK=nK, optimizer_kwargs={"verbose": True, "M": 0}
-    )
+    mf_sol = solver(model, guess, nK=nK, optimizer_kwargs={"verbose": True, "M": 0})
     gap = compute_gap(addTb(h_0, mf_sol), n=100)
 
     # Check if the gap is predicted correctly