fix bug in exchange calculation

8c18865f · Kostas Vilkelis · 4ac2691a · 8c18865f
Commit 8c18865f authored 1 year ago by Kostas Vilkelis
--- a/examples/testInterface.ipynb
+++ b/examples/testInterface.ipynb
@@ -171,7 +171,7 @@
    "    dict\n",
    "        Mean-field tb model.\n",
    "    \"\"\"\n",
-    "    exchange = {vec: -1*(int_model.get(vec, 0) + rhoTbfunc(vec)) for vec in frozenset(int_model)}\n",
+    "    exchange = {vec: -1*(int_model.get(vec, 0)*rhoTbfunc(vec)) for vec in frozenset(int_model)}\n",
    "    # direct interaction is missing\n",
    "    return exchange"
   ]

 %% 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)
 ```