ignore build files

meanfi/model.py

 import numpy as np
+from typing import Tuple
+
 
 from meanfi.mf import (
     density_matrix,
@@ -76,6 +78,25 @@ class Model:
         _check_hermiticity(h_0)
         _check_hermiticity(h_int)
 
+    def density_matrix(self, rho: _tb_type, nk: int = 20) -> Tuple[_tb_type, float]:
+        """Computes the density matrix from a given initial density matrix.
+
+        Parameters
+        ----------
+        rho :
+            Initial density matrix tight-binding dictionary.
+        nk :
+            Number of k-points in a grid to sample the Brillouin zone along each dimension.
+            If the system is 0-dimensional (finite), this parameter is ignored.
+
+        Returns
+        -------
+        :
+            Density matrix tight-binding dictionary.
+        """
+        mf = meanfield(rho, self.h_int)
+        return density_matrix(add_tb(self.h_0, mf), self.filling, nk)[0]
+
     def mfield(self, mf: _tb_type, nk: int = 20) -> _tb_type:
         """Computes a new mean-field correction from a given one.
 

meanfi/solvers.py

@@ -5,11 +5,12 @@ from typing import Optional, Callable
 
 from meanfi.params.rparams import rparams_to_tb, tb_to_rparams
 from meanfi.tb.tb import add_tb, _tb_type
+from meanfi.mf import density_matrix, meanfield
 from meanfi.model import Model
 from meanfi.tb.utils import fermi_energy
 
 
-def cost(mf_param: np.ndarray, model: Model, nk: int = 20) -> np.ndarray:
+def cost_mf(mf_param: np.ndarray, model: Model, nk: int = 20) -> np.ndarray:
     """Defines the cost function for root solver.
 
     The cost function is the difference between the computed and inputted mean-field.
@@ -37,14 +38,45 @@ def cost(mf_param: np.ndarray, model: Model, nk: int = 20) -> np.ndarray:
     return mf_params_new - mf_param
 
 
-def solver(
+def cost_density(rho_params: np.ndarray, model: Model, nk: int = 20) -> np.ndarray:
+    """Defines the cost function for root solver.
+
+    The cost function is the difference between the computed and inputted density matrix
+    reduced to the hoppings present in the h_int.
+
+    Parameters
+    ----------
+    rho_params :
+        1D real array that parametrises the density matrix.
+    Model :
+        Interacting tight-binding problem definition.
+    nk :
+        Number of k-points in a grid to sample the Brillouin zone along each dimension.
+        If the system is 0-dimensional (finite), this parameter is ignored.
+
+    Returns
+    -------
+    :
+        1D real array that is the difference between the computed and inputted mean-field
+        parametrisations
+    """
+    shape = model._ndof
+    rho_red = rparams_to_tb(rho_params, list(model.h_int), shape)
+    rho_new = model.density_matrix(rho_red, nk=nk)
+    rho_red_new = {key: rho_new[key] for key in model.h_int}
+    rho_params_new = tb_to_rparams(rho_red_new)
+    return rho_params_new - rho_params
+
+
+def solver_mf(
     model: Model,
     mf_guess: np.ndarray,
     nk: int = 20,
     optimizer: Optional[Callable] = scipy.optimize.anderson,
     optimizer_kwargs: Optional[dict[str, str]] = {"M": 0},
 ) -> _tb_type:
-    """Solve for the mean-field correction through self-consistent root finding.
+    """Solve for the mean-field correction through self-consistent root finding
+    by finding the mean-field correction fixed point.
 
     Parameters
     ----------
@@ -68,9 +100,57 @@ def solver(
     """
     shape = model._ndof
     mf_params = tb_to_rparams(mf_guess)
-    f = partial(cost, model=model, nk=nk)
+    f = partial(cost_mf, model=model, nk=nk)
     result = rparams_to_tb(
         optimizer(f, mf_params, **optimizer_kwargs), list(model.h_int), shape
     )
     fermi = fermi_energy(add_tb(model.h_0, result), model.filling, nk=nk)
     return add_tb(result, {model._local_key: -fermi * np.eye(model._ndof)})
+
+
+def solver_density(
+    model: Model,
+    mf_guess: _tb_type,
+    nk: int = 20,
+    optimizer: Optional[Callable] = scipy.optimize.anderson,
+    optimizer_kwargs: Optional[dict[str, str]] = {"M": 0},
+) -> _tb_type:
+    """Solve for the mean-field correction through self-consistent root finding
+    by finding the density matrix fixed point.
+
+    Parameters
+    ----------
+    model :
+        Interacting tight-binding problem definition.
+    mf_guess :
+        The initial guess for the mean-field correction in the tight-binding dictionary format.
+    nk :
+        Number of k-points in a grid to sample the Brillouin zone along each dimension.
+        If the system is 0-dimensional (finite), this parameter is ignored.
+    optimizer :
+        The solver used to solve the fixed point iteration.
+        Default uses `scipy.optimize.anderson`.
+    optimizer_kwargs :
+        The keyword arguments to pass to the optimizer.
+
+    Returns
+    -------
+    :
+        Mean-field correction solution in the tight-binding dictionary format.
+    """
+    shape = model._ndof
+    rho_guess = density_matrix(
+        add_tb(model.h_0, mf_guess), filling=model.filling, nk=nk
+    )[0]
+    rho_guess_red = {key: rho_guess[key] for key in model.h_int}
+
+    rho_params = tb_to_rparams(rho_guess_red)
+    f = partial(cost_density, model=model, nk=nk)
+    rho_result = rparams_to_tb(
+        optimizer(f, rho_params, **optimizer_kwargs), list(model.h_int), shape
+    )
+    mf_result = meanfield(rho_result, model.h_int)
+    return mf_result
+
+
+solver = solver_density

meanfi/tests/test_graphene.py

@@ -5,6 +5,7 @@ import pytest
 from meanfi.kwant_helper import kwant_examples, utils
 from meanfi import (
     Model,
+    fermi_energy,
     solver,
     tb_to_kgrid,
     guess_tb,
@@ -78,7 +79,9 @@ def gap_prediction(U, V):
     model = Model(h_0, h_int, filling)
 
     mf_sol = solver(model, guess, nk=nk, optimizer_kwargs={"M": 0, "f_tol": 1e-8})
-    gap = compute_gap(add_tb(h_0, mf_sol), nk=200)
+    mf_full = add_tb(h_0, mf_sol)
+    EF = fermi_energy(mf_full, filling=filling, nk=nk)
+    gap = compute_gap(mf_full, fermi_energy=EF, nk=200)
 
     # Check if the gap is predicted correctly
     if gap > 0.1:

meanfi/tests/test_hubbard.py

@@ -5,6 +5,7 @@ import pytest
 from meanfi.tests.test_graphene import compute_gap
 from meanfi import (
     Model,
+    fermi_energy,
     solver,
     guess_tb,
     add_tb,
@@ -38,7 +39,9 @@ def gap_relation_hubbard(Us, nk, nk_dense, tol=1e-3):
         guess = guess_tb(frozenset(h_int), len(list(h_0.values())[0]))
         full_model = Model(h_0, h_int, filling=2)
         mf_sol = solver(full_model, guess, nk=nk)
-        _gap = compute_gap(add_tb(h_0, mf_sol), fermi_energy=0, nk=nk_dense)
+        mf_full = add_tb(h_0, mf_sol)
+        EF = fermi_energy(mf_full, filling=2, nk=nk)
+        _gap = compute_gap(mf_full, fermi_energy=EF, nk=nk_dense)
         gaps.append(_gap)
 
     fit_gap = np.polyfit(Us, np.array(gaps), 1)[0]