Interface refactoring

Find a design which is more intuitive and minimal.

codes/hf.py

@@ -104,17 +104,11 @@ def compute_mf(rho, H_int):
         local_density = np.diag(np.average(rho, axis=tuple([i for i in range(dim)])))
         exchange_mf = convolution(rho, H_int) * nk ** (-dim)
         direct_mf = np.diag(np.einsum("i,ij->j", local_density, H0_int))
-        dc_direct = local_density.T @ H0_int @ local_density
-        dc_exchange = np.einsum('kij, kji', exchange_mf, rho) * nk ** (-dim)
-        dc_energy = 0.5*(-dc_exchange + dc_direct) * np.eye(direct_mf.shape[-1])
     else:
         local_density = np.diag(rho)
         exchange_mf = rho * H_int
         direct_mf = np.diag(np.einsum("i,ij->j", local_density, H_int))
-        dc_energy_direct = np.diag(np.einsum("ij, i, j->", H_int, local_density, local_density))
-        dc_energy_cross = np.diag(np.einsum("ij, ij, ji->", H_int, rho, rho))
-        dc_energy = 2 * dc_energy_direct - dc_energy_cross
-    return direct_mf - exchange_mf# - dc_energy * np.eye(direct_mf.shape[-1])
+    return direct_mf - exchange_mf
 
 def total_energy(h, rho):
     """
@@ -132,35 +126,4 @@ def total_energy(h, rho):
     total_energy : float
         System total energy computed as tr[h@rho].
     """
-    return np.sum(np.trace(h @ rho, axis1=-1, axis2=-2)).real / np.prod(rho.shape[:-2])
-
-def updated_matrices(mf_k, model):
-    """
-    Self-consistent loop.
-
-    Parameters:
-    -----------
-    mf : nd-array
-        Mean-field correction. Same format as the initial guess.
-    H_int : nd-array
-        Interaction matrix.
-    filling: int
-        Number of electrons per cell.
-    hamiltonians_0 : nd-array
-        Non-interacting Hamiltonian. Same format as `H_int`.
-
-    Returns:
-    --------
-    mf_new : nd-array
-        Updated mean-field solution.
-    """
-    # Generate the Hamiltonian
-    hamiltonians = model.hamiltonians_0 + mf_k
-    vals, vecs = np.linalg.eigh(hamiltonians)
-    vecs = np.linalg.qr(vecs)[0]
-    E_F = utils.get_fermi_energy(vals, model.filling)
-    rho = density_matrix(vals=vals, vecs=vecs, E_F=E_F)
-    return rho, compute_mf(
-        rho=rho,
-        H_int=model.H_int) - E_F * np.eye(model.hamiltonians_0.shape[-1])
-
+    return np.sum(np.trace(h @ rho, axis1=-1, axis2=-2)).real / np.prod(rho.shape[:-2])
\ No newline at end of file