diff --git a/codes/interface.py b/codes/interface.py
index 77637b34ad9057343ca6b1354565b7f844349d9d..21c6edc6bb3c7c2b09c6a819c10c4cab784795b7 100644
--- a/codes/interface.py
+++ b/codes/interface.py
@@ -46,6 +46,9 @@ def find_groundstate_ham(
     solver(model, optimizer, cost_function, optimizer_kwargs)
     model.vectors=[*model.vectors, *model.tb_model.keys()]
     assert np.allclose(model.mf_k - np.moveaxis(model.mf_k, -1, -2).conj(), 0, atol=1e-15)
+    vals, _ = np.linalg.eigh(model.hamiltonians_0 + model.mf_k)
+    EF = utils.get_fermi_energy(vals, filling)
+    model.mf_k -=  EF * np.eye(model.hamiltonians_0.shape[-1])
     if return_kspace:
         return model.hamiltonians_0 + model.mf_k
     else: