Interface refactoring

Find a design which is more intuitive and minimal.

codes/kwant_helper/utils.py

 import numpy as np
 import kwant
-from itertools import product
 from scipy.sparse import coo_array
+from itertools import product
 import inspect
 from copy import copy
 
 
-def get_fermi_energy(vals, filling):
-    """
-    Compute Fermi energy for a given filling factor.
-
-    vals : nd-array
-        Collection of eigenvalues on a grid.
-    filling : int
-        Number of electrons per cell.
-    """
-    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 builder2tb_model(builder, params={}, return_data=False):
+def builder_to_tb(builder, params={}, return_data=False):
     """
     Constructs a tight-binding model dictionary from a `kwant.Builder`.
 
@@ -43,7 +21,7 @@ def builder2tb_model(builder, params={}, return_data=False):
 
     Returns:
     --------
-    tb_model : dict
+    h_0 : dict
         Tight-binding model of non-interacting systems.
     data : dict
         Data with sites and number of orbitals. Only if `return_data=True`.
@@ -52,7 +30,7 @@ def builder2tb_model(builder, params={}, return_data=False):
     # Extract information from builder
     dims = len(builder.symmetry.periods)
     onsite_idx = tuple([0] * dims)
-    tb_model = {}
+    h_0 = {}
     sites_list = [*builder.sites()]
     norbs_list = [site[0].norbs for site in builder.sites()]
     positions_list = [site[0].pos for site in builder.sites()]
@@ -73,14 +51,14 @@ def builder2tb_model(builder, params={}, return_data=False):
                     _params[arg] = params[arg]
             val = val(site, **_params)
             data = val.flatten()
-        except:
+        except Exception:
             data = val.flatten()
-        if onsite_idx in tb_model:
-            tb_model[onsite_idx] += coo_array(
+        if onsite_idx in h_0:
+            h_0[onsite_idx] += coo_array(
                 (data, (row, col)), shape=(norbs_tot, norbs_tot)
             ).toarray()
         else:
-            tb_model[onsite_idx] = coo_array(
+            h_0[onsite_idx] = coo_array(
                 (data, (row, col)), shape=(norbs_tot, norbs_tot)
             ).toarray()
     # Hopping matrices
@@ -101,24 +79,37 @@ def builder2tb_model(builder, params={}, return_data=False):
                     _params[arg] = params[arg]
             val = val(a, b, **_params)
             data = val.flatten()
-        except:
+        except Exception:
             data = val.flatten()
-        if tuple(b_dom) in tb_model:
-            tb_model[tuple(b_dom)] += coo_array(
+        if tuple(b_dom) in h_0:
+            h_0[tuple(b_dom)] += coo_array(
                 (data, (row, col)), shape=(norbs_tot, norbs_tot)
             ).toarray()
             if np.linalg.norm(b_dom) == 0:
-                tb_model[tuple(b_dom)] += (
+                h_0[tuple(b_dom)] += (
+                    coo_array((data, (row, col)), shape=(norbs_tot, norbs_tot))
+                    .toarray()
+                    .T.conj()
+                )
+            else:
+                # Hopping vector in the opposite direction
+                h_0[tuple(-b_dom)] += (
                     coo_array((data, (row, col)), shape=(norbs_tot, norbs_tot))
                     .toarray()
                     .T.conj()
                 )
         else:
-            tb_model[tuple(b_dom)] = coo_array(
+            h_0[tuple(b_dom)] = coo_array(
                 (data, (row, col)), shape=(norbs_tot, norbs_tot)
             ).toarray()
             if np.linalg.norm(b_dom) == 0:
-                tb_model[tuple(b_dom)] += (
+                h_0[tuple(b_dom)] += (
+                    coo_array((data, (row, col)), shape=(norbs_tot, norbs_tot))
+                    .toarray()
+                    .T.conj()
+                )
+            else:
+                h_0[tuple(-b_dom)] = (
                     coo_array((data, (row, col)), shape=(norbs_tot, norbs_tot))
                     .toarray()
                     .T.conj()
@@ -128,84 +119,9 @@ def builder2tb_model(builder, params={}, return_data=False):
         data = {}
         data["norbs"] = norbs_list
         data["positions"] = positions_list
-        return tb_model, data
-    else:
-        return tb_model
-
-
-def dict2hk(tb_dict):
-    """
-    Build Bloch Hamiltonian.
-
-    Returns:
-    --------
-    bloch_ham : function
-        Evaluates the Hamiltonian at a given k-point.
-    """
-
-    def bloch_ham(k):
-        ham = 0
-        for vector in tb_dict.keys():
-            ham += tb_dict[vector] * np.exp(
-                1j * np.dot(k, np.array(vector, dtype=float))
-            )
-            if np.linalg.norm(np.array(vector)) > 0:
-                ham += tb_dict[vector].T.conj() * np.exp(
-                    -1j * np.dot(k, np.array(vector))
-                )
-        return ham
-
-    return bloch_ham
-
-
-def kgrid_hamiltonian(nk, tb_model, return_ks=False):
-    """
-    Evaluates Hamiltonian on a k-point grid.
-
-    Paramters:
-    ----------
-    nk : int
-        Number of k-points along each direction.
-    tb_model : dictionary
-        Must have the following structure:
-            - Keys are tuples for each hopping vector (in units of lattice vectors).
-            - Values are hopping matrices.
-    return_ks : bool
-        Return k-points.
-
-    Returns:
-    --------
-    ham : nd.array
-        Hamiltonian evaluated on a k-point grid from k-points
-        along each direction evaluated from zero to 2*pi.
-        The indices are ordered as [k_1, ... , k_n, i, j], where
-        `k_m` corresponding to the k-point element along each
-        direction and `i` and `j` are the internal degrees of freedom.
-    ks : 1D-array
-        List of k-points over all directions. Only returned if `return_ks=True`.
-    """
-    dim = len(next(iter(tb_model)))
-    if dim == 0:
-        if return_ks:
-            return syst[next(iter(tb_model))], None
-        else:
-            return syst[next(iter(tb_model))]
-    else:
-        hk = dict2hk(tb_model)
-
-    ks = 2 * np.pi * np.linspace(0, 1, nk, endpoint=False)
-
-    k_pts = np.tile(ks, dim).reshape(dim, nk)
-
-    ham = []
-    for k in product(*k_pts):
-        ham.append(hk(k))
-    ham = np.array(ham)
-    shape = (*[nk] * dim, ham.shape[-1], ham.shape[-1])
-    if return_ks:
-        return ham.reshape(*shape), ks
+        return h_0, data
     else:
-        return ham.reshape(*shape)
+        return h_0
 
 
 def build_interacting_syst(builder, lattice, func_onsite, func_hop, max_neighbor=1):
@@ -235,105 +151,3 @@ def build_interacting_syst(builder, lattice, func_onsite, func_hop, max_neighbor
     for neighbors in range(max_neighbor):
         int_builder[lattice.neighbors(neighbors + 1)] = func_hop
     return int_builder
-
-
-def generate_guess(tb_model, int_model, scale=0.1):
-    """
-    tb_model : dict
-        Tight-binding model of non-interacting systems.
-    int_model : dict
-        Tight-binding model for interacting Hamiltonian.
-    scale : float
-        The scale of the guess. Maximum absolute value of each element of the guess.
-
-    Returns:
-    --------
-    guess : tb dictionary
-        TB guess.
-    """
-    ndof = tb_model[next(iter(tb_model))].shape[-1]
-    guess = {}
-    vectors = frozenset(tb_model) | frozenset(int_model)
-    for vector in vectors:
-        amplitude = np.random.rand(ndof, ndof)
-        phase = 2 * np.pi * np.random.rand(ndof, ndof)
-        rand_hermitian = amplitude * np.exp(1j * phase)
-        guess[vector] = rand_hermitian * scale
-        op_vector = tuple(-np.array(vector))
-        if op_vector==vector:
-            guess[vector] += guess[vector].T.conj()
-            guess[vector] /= 2
-        else:
-            guess[op_vector] = guess[vector].T.conj()
-    return guess
-
-
-def hk2tb_model(hk, tb_model, int_model, ks=None):
-    """
-    Extract hopping matrices from Bloch Hamiltonian.
-
-    Parameters:
-    -----------
-    hk : nd-array
-        Bloch Hamiltonian matrix hk[k_x, ..., k_n, i, j]
-    tb_model : dict
-        Tight-binding model of non-interacting systems.
-    int_model : dict
-        Tight-binding model for interacting Hamiltonian.
-    ks : 1D-array
-        Set of k-points. Repeated for all directions. If the system is finite, `ks=None`.
-
-    Returns:
-    --------
-    scf_model : dict
-        TIght-binding model of Hartree-Fock solution.
-    """
-    if ks is not None:
-        hopping_vecs = np.unique(
-            np.array([*tb_model.keys(), *int_model.keys()]), axis=0
-        )
-        ndim = len(hk.shape) - 2
-        dk = np.diff(ks)[0]
-        nk = len(ks)
-        k_pts = np.tile(ks, ndim).reshape(ndim, nk)
-        k_grid = np.array(np.meshgrid(*k_pts))
-        k_grid = k_grid.reshape(k_grid.shape[0], np.prod(k_grid.shape[1:]))
-        hk = hk.reshape(np.prod(hk.shape[:ndim]), *hk.shape[-2:])
-
-        hopps = (
-            np.einsum(
-                "ij,jkl->ikl",
-                np.exp(1j * np.einsum("ij,jk->ik", hopping_vecs, k_grid)),
-                hk,
-            )
-            * (dk / (2 * np.pi)) ** ndim
-        )
-
-        tb_model = {}
-        for i, vector in enumerate(hopping_vecs):
-            tb_model[tuple(vector)] = hopps[i]
-
-        return tb_model
-    else:
-        return {(): hk}
-
-
-def calc_gap(vals, E_F):
-    """
-    Compute gap.
-
-    Parameters:
-    -----------
-    vals : nd-array
-        Eigenvalues on a k-point grid.
-    E_F : float
-        Fermi energy.
-
-    Returns:
-    --------
-    gap : float
-        Indirect gap.
-    """
-    emax = np.max(vals[vals <= E_F])
-    emin = np.min(vals[vals > E_F])
-    return np.abs(emin - emax)