Skip to content
Snippets Groups Projects

Interface refactoring

Merged Interface refactoring
interface-refactoring into main
Merged Kostas Vilkelis requested to merge interface-refactoring into main
Compare
and
Show latest version
26 files
+ 743
− 1623
Compare changes
  • Side-by-side
  • Inline
Files
26
codes/kwant_helper/utils.py
+ 3
− 279
Antonio Manesco
Last comment by Johanna Zijderveld
@@ -6,29 +6,7 @@ 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 builder2h_0(builder, params={}, return_data=False):
def builder_to_tb(builder, params={}, return_data=False):
"""
Constructs a tight-binding model dictionary from a `kwant.Builder`.
@@ -73,7 +51,7 @@ def builder2h_0(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 h_0:
h_0[onsite_idx] += coo_array(
@@ -101,7 +79,7 @@ def builder2h_0(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 h_0:
h_0[tuple(b_dom)] += coo_array(
@@ -146,93 +124,6 @@ def builder2h_0(builder, params={}, return_data=False):
return h_0
def model2hk(h_0):
"""
Build Bloch Hamiltonian.
Paramters:
----------
nk : int
Number of k-points along each direction.
h_0 : 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`.
Returns:
--------
bloch_ham : function
Evaluates the Hamiltonian at a given k-point.
"""
assert (
len(next(iter(h_0))) > 0
), "Zero-dimensional system. The Hamiltonian is simply h_0[()]."
def bloch_ham(k):
ham = 0
for vector in h_0.keys():
ham += h_0[vector] * np.exp(-1j * np.dot(k, np.array(vector, dtype=float)))
return ham
return bloch_ham
def kgrid_hamiltonian(nk, hk, dim, return_ks=False, hermitian=True):
"""
Evaluates Hamiltonian on a k-point grid.
Paramters:
----------
nk : int
Number of k-points along each direction.
hk : function
Calculates the Hamiltonian at a given k-point.
return_ks : bool
If `True`, returns 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`.
"""
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)
if hermitian:
assert np.allclose(
ham, np.transpose(ham, (0, 2, 1)).conj()
), "Tight-binding provided is non-Hermitian. Not supported yet"
shape = (*[nk] * dim, ham.shape[-1], ham.shape[-1])
if return_ks:
return ham.reshape(*shape), ks
else:
return ham.reshape(*shape)
def build_interacting_syst(builder, lattice, func_onsite, func_hop, max_neighbor=1):
"""
Construct an auxiliary `kwant` system to build Hamiltonian matrix.
@@ -260,170 +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(vectors, ndof, scale=1):
"""
vectors : list
List of hopping vectors.
ndof : int
Number internal degrees of freedom (orbitals),
scale : float
The scale of the guess. Maximum absolute value of each element of the guess.
Returns:
--------
guess : tb dictionary
Guess in the form of a tight-binding model.
"""
guess = {}
for vector in vectors:
if vector not in guess.keys():
amplitude = scale * np.random.rand(ndof, ndof)
phase = 2 * np.pi * np.random.rand(ndof, ndof)
rand_hermitian = amplitude * np.exp(1j * phase)
if np.linalg.norm(np.array(vector)) == 0:
rand_hermitian += rand_hermitian.T.conj()
rand_hermitian /= 2
guess[vector] = rand_hermitian
else:
guess[vector] = rand_hermitian
guess[tuple(-np.array(vector))] = rand_hermitian.T.conj()
return guess
def generate_vectors(cutoff, dim):
"""
Generates hopping vectors up to a cutoff.
Parameters:
-----------
cutoff : int
Maximum distance along each direction.
dim : int
Dimension of the vectors.
Returns:
--------
List of hopping vectors.
"""
return [*product(*([[*range(-cutoff, cutoff + 1)]] * dim))]
def hk2h_0(hk, hopping_vecs, ks=None):
"""
Extract hopping matrices from Bloch Hamiltonian.
Parameters:
-----------
hk : nd-array
Bloch Hamiltonian matrix hk[k_x, ..., k_n, i, j]
h_0 : dict
Tight-binding model of non-interacting systems.
h_int : 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:
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
)
h_0 = {}
for i, vector in enumerate(hopping_vecs):
h_0[tuple(vector)] = hopps[i]
return h_0
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)
def matrix_to_flat(matrix):
"""
Flatten the upper triangle of a collection of matrices.
Parameters:
-----------
matrix : nd-array
Array with shape (..., n, n)
"""
return matrix[..., *np.triu_indices(matrix.shape[-1])].flatten()
def flat_to_matrix(flat, shape):
"""
Undo `matrix_to_flat`.
Parameters:
-----------
flat : 1d-array
Output from `matrix_to_flat`.
shape : 1d-array
Shape of the input from `matrix_to_flat`.
"""
matrix = np.zeros(shape, dtype=complex)
matrix[..., *np.triu_indices(shape[-1])] = flat.reshape(*shape[:-2], -1)
indices = np.arange(shape[-1])
diagonal = matrix[..., indices, indices]
matrix += np.moveaxis(matrix, -1, -2).conj()
matrix[..., indices, indices] -= diagonal
return matrix
def complex_to_real(z):
"""
Split real and imaginary parts of a complex array.
Parameters:
-----------
z : array
"""
return np.concatenate((np.real(z), np.imag(z)))
def real_to_complex(z):
"""
Undo `complex_to_real`.
"""
return z[: len(z) // 2] + 1j * z[len(z) // 2 :]