Interface refactoring

Find a design which is more intuitive and minimal.

examples/testInterface.ipynb

 %% Cell type:code id: tags:
 
 ``` python
 import numpy as np
 from scipy.integrate import quad, quad_vec
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def tb2kfunc(tb_model):
     """
     Fourier transforms a real-space tight-binding model to a k-space function.
 
     Parameters
     ----------
     tb_model : dict
         A dictionary with real-space vectors as keys and complex np.arrays as values.
 
     Returns
     -------
     function
         A function that takes a k-space vector and returns a complex np.array.
     """
 
     def bloch_ham(k):
         ham = 0
         for vector in tb_model.keys():
             ham += tb_model[vector] * np.exp(
                 -1j * np.dot(k, np.array(vector, dtype=float))
             )
         return ham
 
     return bloch_ham
 
 
 def kfunc2tbfunc(kfunc):
     """
     Inverse Fourier transforms a k-space function to a real-space function.
 
     Parameters
     ----------
     kfunc : function
         A function that takes a k-space vector and returns a np.array.
 
     Returns
     -------
     function
         A function that takes a real-space integer vector and returns a np.array.
     """
     def tbfunc(vector):
 
         def integrand(k):
             return kfunc(k) * np.exp(1j * np.dot(k, np.array(vector, dtype=float)))
 
         return quad_vec(integrand, -np.pi, np.pi)[0]
 
     return tbfunc
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def addTb(tb1, tb2):
     """
     Add up two tight-binding models together.
 
     Parameters:
     -----------
     tb1 : dict
         Tight-binding model.
     tb2 : dict
         Tight-binding model.
 
     Returns:
     --------
     dict
         Sum of the two tight-binding models.
     """
     return {k: tb1.get(k, 0) + tb2.get(k, 0) for k in frozenset(tb1) | frozenset(tb2)}
 
 def generateRho(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
     -------
     rho : function
         Density matrix at the same k-point as hk
     """
     def rhofunc(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 rhofunc
 
 def fermiOnGrid(hkfunc, filling, Nk=100):
     """
     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
     """
 
     ks = np.linspace(-np.pi, np.pi, Nk, endpoint=False)
     hkarray = np.array([hkfunc(k) for k in ks])
     vals = np.linalg.eigvalsh(hkarray)
 
     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 totalEnergy(rho, hk):
     def integrand(k):
         return np.real(np.trace(rho(k) @ hk(k)))
     return quad(integrand, -np.pi, np.pi)[0]
 
 def meanField(rhoTbfunc, int_model):
     """
     Compute the mean-field in real space.
 
     Parameters
     ----------
     rhoTbfunc : function
         Function that returns the density matrix at a given real-space vector.
     int_model : dict
         Interaction tb model.
 
     Returns
     -------
     dict
         Mean-field tb model.
     """
-    exchange = {vec: -1*(int_model.get(vec, 0) + rhoTbfunc(vec)) for vec in frozenset(int_model)}
+    exchange = {vec: -1*(int_model.get(vec, 0)*rhoTbfunc(vec)) for vec in frozenset(int_model)}
     # direct interaction is missing
     return exchange
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 U0 = 1
 nk = 100
 filling = 2
 
 hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
 tb_model = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
 mf_model = {} # placeholder for guess
 int_model = {
     (0,): U0 * np.kron(np.eye(2), np.ones((2, 2))),
 }
 
 # attempting to do a single iteration of mean-field
 hkfunc = tb2kfunc(addTb(tb_model, mf_model))
 EF = fermiOnGrid(hkfunc, filling, nk)
 rhofunc = generateRho(hkfunc, EF)
 rhoTbfunc = kfunc2tbfunc(rhofunc)
 mf_model = meanField(rhoTbfunc, int_model)
 ```