start packaging

codes/hf.py

+from scipy.signal import convolve2d
+import numpy as np
+import utils
+
+def mean_field_F(vals, vecs, E_F):
+    N_ks = vecs.shape[0]
+    unocc_vals = vals > E_F
+
+    def mf_generator(i, j):
+        occ_vecs = vecs[i, j]
+        occ_vecs[:, unocc_vals[i, j]] = 0
+        F_ij = occ_vecs @ occ_vecs.conj().T
+        return F_ij
+
+    F = np.array([[mf_generator(i, j) for i in range(N_ks)] for j in range(N_ks)])
+    return F
+
+def convolution(M1, M2):
+    cell_size = M2.shape[-1]
+    V_output = np.array(
+        [
+            [
+                convolve2d(M1[:, :, i, j], M2[:, :, i, j], boundary="wrap", mode="same")
+                for i in range(cell_size)
+            ]
+            for j in range(cell_size)
+        ]
+    )
+    V_output = np.transpose(V_output, axes=(2, 3, 0, 1))
+    return V_output
+
+def compute_mf(vals, vecs, filling, nk, H_int):
+    H0_int = H_int[0,0]
+    E_F = utils.get_fermi_energy(vals, filling)
+    F = mean_field_F(vals, vecs, E_F=E_F)
+    rho = np.diag(np.average(F, axis=(0, 1)))
+    exchange_mf = convolution(F, H_int) * nk ** (-2)
+    direct_mf = np.diag(np.einsum("i,ij->j", rho, H0_int))
+
+    return direct_mf - exchange_mf
\ No newline at end of file

codes/utils.py

+import numpy as np
+import kwant
+
+def get_fermi_energy(vals, filling):
+    norbs = vals.shape[-1]
+    vals_flat = np.sort(vals.flatten())
+    ne = len(vals_flat)
+    ifermi = int(round(ne * filling / norbs))
+    if ifermi >= ne:
+        ifermi = ne - 1
+    sorte = np.sort(vals_flat)  # sorted eigenvalues
+    if ifermi == 0:
+        return sorte[0]
+    fermi = (sorte[ifermi - 1] + sorte[ifermi]) / 2.0  # fermi energy
+    return fermi
+
+def syst2hamiltonian(kxs, kys, syst, params={}):
+    def h_k(kx, ky):
+        return syst.hamiltonian_submatrix(params={**params, **dict(k_x=kx, k_y=ky)})
+    return np.array(
+        [[h_k(kx, ky) for kx in kxs] for ky in kys]
+    )
+
+def potential2hamiltonian(
+    syst, lattice, func_onsite, func_hop, ks, params={}, max_neighbor=1
+):
+    V = kwant.Builder(kwant.TranslationalSymmetry(*lattice.prim_vecs))
+    V[syst.sites()] = func_onsite
+    for neighbors in range(max_neighbor):
+        V[lattice.neighbors(neighbors + 1)] = func_hop
+    wrapped_V = kwant.wraparound.wraparound(V).finalized()
+    return syst2hamiltonian(kxs=ks, kys=ks, syst=wrapped_V, params=params)
+
+def generate_guess(max_neighbor, norbs, lattice):
+    n_sub = len(lattice.sublattices)
+    guess = np.zeros((n_sub + max_neighbor, 2, norbs, norbs))
+    for i in range(n_sub):
+        # Real part
+        guess[i, 0] = np.random.rand(norbs, norbs)
+        guess[i, 0] += guess[i, 0].T
+        # Imag part
+        guess[i, 1] = np.random.rand(norbs, norbs)
+        guess[i, 1] -= guess[i, 1].T
+    for neighbor in range(max_neighbor):
+        # Real part
+        guess[n_sub + neighbor, 0] = np.random.rand(norbs, norbs)
+        # Imag part
+        guess[n_sub + neighbor, 1] = np.random.rand(norbs, norbs)
+    return guess
+
+def extract_hopping_vectors(builder):
+    keep=None
+    deltas=[]
+    for hop, val in builder.hopping_value_pairs():
+        a, b=hop
+        b_dom = builder.symmetry.which(b)
+        # Throw away part that is in the remaining translation direction, so we get
+        # an element of 'sym' which is being wrapped
+        b_dom = np.array([t for i, t in enumerate(b_dom) if i != keep])
+        deltas.append(b_dom)
+    return np.asarray(deltas)
+
+def generate_scf_syst(max_neighbor, syst, lattice):
+    subs = np.array(lattice.sublattices)
+    def scf_onsite(site, mat):
+        idx = np.where(subs == site.family)[0][0]
+        return mat[idx, 0] + 1j * mat[idx, 1]
+
+    scf = kwant.Builder(kwant.TranslationalSymmetry(*lattice.prim_vecs))
+    scf[syst.sites()] = scf_onsite
+    for neighbor in range(max_neighbor):
+
+        def scf_hopping(site1, site2, mat):
+            return (
+                mat[len(lattice.sublattices) + neighbor, 0]
+                + 1j * mat[len(lattice.sublattices) + neighbor, 1]
+            )
+
+        scf[lattice.neighbors(neighbor + 1)] = scf_hopping
+    deltas = extract_hopping_vectors(scf)
+    wrapped_scf = kwant.wraparound.wraparound(scf).finalized()
+    return wrapped_scf, deltas
\ No newline at end of file