Skip to content
GitLab
Commit 403bff26 authored by Benoit Rossignol's avatar Benoit Rossignol
Browse files

API of fast, pycode everywhere

parent 47f2c7fb
Hide whitespace changes
Inline Side-by-side
boundstate/__init__.py
View file @ 403bff26
from . import fast, standard
\ No newline at end of file
"""
This library is used to find and compute boundstat wave function.
'standard' module use SVD decompostion.
'fast' module use eigenvalue decomposition.
It is faster that SVD but may return false positive.
"""
from . import fast, standard
boundstate/_common.py
View file @ 403bff26
......@@ -7,12 +7,14 @@ from collections import defaultdict
from scipy.spatial.distance import cdist
from scipy.sparse.csgraph import connected_components
def list_duplicates(seq):
tally = defaultdict(list)
for i,item in enumerate(seq):
for i, item in enumerate(seq):
tally[item].append(i)
return ((key,locs) for key,locs in tally.items()
if len(locs)>1)
return ((key, locs) for key, locs in tally.items()
if len(locs) > 1)
def null_space(V, eps=1e-8):
"""
......@@ -24,8 +26,8 @@ def null_space(V, eps=1e-8):
Q: 2D numpy array, same dimensions as V
"""
U, S, WH = la.svd(V)
assert(np.allclose(np.conj(np.transpose(WH)).dot(WH), np.eye(WH.shape[0]))), \
'SVD matrices not unitary'
assert(np.allclose(np.conj(np.transpose(WH)).dot(WH),
np.eye(WH.shape[0]))), 'SVD matrices not unitary'
zero_eig_val = (S < eps)
dim = len(S)
Q = np.zeros(shape=V.shape, dtype=complex)
......@@ -34,11 +36,14 @@ def null_space(V, eps=1e-8):
x_alpha = WH[alpha, :]
Q += np.dot(np.conj(x_alpha[:, None]), x_alpha[None, :])
# two checks, optionnal
#~ assert(np.allclose(np.dot(Q, Q), Q)), 'Projector is not idempotent (P**2 =/= P)'
#~ assert(np.allclose(np.dot(V, Q), np.zeros(shape=V.shape), atol=1e-4)), 'Projection on the kernel is wrong'
# ~ assert(np.allclose(np.dot(Q, Q), Q)), 'Projector is not
# idempotent (P**2 =/= P)'
# ~ assert(np.allclose(np.dot(V, Q), np.zeros(shape=V.shape), atol=1e-4)),
# 'Projection on the kernel is wrong'
return Q
def fill_zero_with_nan(array):
"""
Replace all elements equal to 0 by nan so that it is
......@@ -53,20 +58,22 @@ def fill_zero_with_nan(array):
array[i] = float('nan')
return np.reshape(array, (S[0], S[1]))
def leads_modes(H_leads, V_leads, E, tol=1e-5, need_orth=True):
"""
Compute the evanescent modes for every lead separetly
"""
X_out_leads = []
L_out_leads = []
assert len(H_leads)==len(V_leads), 'Check the number of leads'
assert len(H_leads) == len(V_leads), 'Check the number of leads'
for H, V in zip(H_leads, V_leads): #compute the modes of every lead
for H, V in zip(H_leads, V_leads): # compute the modes of every lead
L_out, X_out = extract_out_modes(H, V, E, tol=tol, need_orth=need_orth)
X_out_leads.append(X_out)
L_out_leads.append(L_out)
return L_out_leads, X_out_leads
def block_diag(*arrs):
"""
Copied from old version of scipy, v0.8.0
......@@ -126,7 +133,7 @@ def block_diag(*arrs):
bad_args = [k for k in range(len(arrs)) if arrs[k].ndim > 2]
if bad_args:
raise ValueError("arguments in the following positions have dimension "
"greater than 2: %s" % bad_args)
"greater than 2: %s" % bad_args)
shapes = np.array([a.shape for a in arrs])
out = np.zeros(np.sum(shapes, axis=0), dtype=arrs[0].dtype)
......@@ -138,6 +145,7 @@ def block_diag(*arrs):
c += cc
return out
def extract_out_modes(H, V, E, tol=1e-6, return_mat=False, need_orth=True):
"""
parameters:
......@@ -157,7 +165,8 @@ def extract_out_modes(H, V, E, tol=1e-6, return_mat=False, need_orth=True):
V_dag = V.conj().T
matrices, v, extract = setup_linsys(H, V, stabilization=(True, True))
all_lmb, evanselect, propselect, vec_gen, ord_schur = unified_eigenproblem(*(matrices))
all_lmb, evanselect,\
propselect, vec_gen, ord_schur = unified_eigenproblem(*(matrices))
lmb_inv = all_lmb[evanselect]
all_vecs = vec_gen(evanselect)
Phi = np.zeros(shape=(n, len(lmb_inv)), dtype=complex)
......@@ -169,32 +178,35 @@ def extract_out_modes(H, V, E, tol=1e-6, return_mat=False, need_orth=True):
Phi, arg_max, max_x = max_norm(Phi)
#function that test if there is a degeneracy
# function that test if there is a degeneracy
L = np.diag(1 / lmb_inv)
if need_orth:
if len(lmb_inv):
all_dupl = find_groups(lmb_inv, tol=tol)
for dupl in all_dupl: # test is there are similar eigenvalues
for dupl in all_dupl: # test is there are similar eigenvalues
Phi[:, dupl], R = np.linalg.qr(Phi[:, dupl])
temp = Phi[:, dupl].conj().T @ Phi[:, dupl]
else: L = np.zeros(shape=(len(lmb_inv), 0), dtype=complex)
else:
L = np.zeros(shape=(len(lmb_inv), 0), dtype=complex)
if return_mat:
return L, Phi, R
else:
return L, Phi
def find_groups(arr, tol=1e-6):
"""
Return:
all_dupl: list countaining lists with indices of equal eigenvalues
up to the tolerance tol
"""
#function that search if other_idx is already in all_dupl, modifies all_dupl
# function that search if other_idx is already in all_dupl, modifies
# all_dupl
def already_dupl(all_dupl, idx, idx_min):
for i, other_dupl in enumerate(all_dupl):
if any([idx in other_dupl]):
#a lready belong to another doublon, triplet...
# a lready belong to another doublon, triplet...
all_dupl[i].append(idx_min)
return True
return False
......@@ -218,6 +230,7 @@ def find_groups(arr, tol=1e-6):
all_dupl.append(dupl)
return all_dupl
def extract_kwant_matrices(sys, args=(), sparse=False):
"""
Function that takes a finalized kwant system (with at least one lead)
......@@ -251,10 +264,11 @@ def extract_kwant_matrices(sys, args=(), sparse=False):
V_leads.append(sys.leads[i].inter_cell_hopping())
coords = np.concatenate((coords, sys.lead_interfaces[i]))
norb = sys.hamiltonian_submatrix(args=args, sparse=sparse, return_norb=True)[1]
norb = sys.hamiltonian_submatrix(args=args, sparse=sparse,
return_norb=True)[1]
offsets = np.empty(norb.shape[0] + 1, int)
offsets[0] = 0
offsets[1 :] = np.cumsum(norb)
offsets[1:] = np.cumsum(norb)
ones = []
for coord in coords:
coord = int(coord)
......@@ -262,15 +276,16 @@ def extract_kwant_matrices(sys, args=(), sparse=False):
ones.append(i)
if sparse:
transf = sp.csr_matrix((np.ones(len(ones)), (range(len(ones)),ones)),
shape=(len(ones),sum(norb)))
transf = sp.csr_matrix((np.ones(len(ones)), (range(len(ones)), ones)),
shape=(len(ones), sum(norb)))
else:
# most inefficient way to create a dense array, to change at some point
transf = sp.csr_matrix((np.ones(len(ones)), (range(len(ones)),ones)),
shape=(len(ones),sum(norb))).toarray()
transf = sp.csr_matrix((np.ones(len(ones)), (range(len(ones)), ones)),
shape=(len(ones), sum(norb))).toarray()
return H_s, H_leads, V_leads, transf
def compute_wf(vals, vecs, L_leads, X_out_leads, scat_dim,
eps=1e-4, schur_modes=False):
"""
......@@ -295,7 +310,8 @@ def compute_wf(vals, vecs, L_leads, X_out_leads, scat_dim,
# When modes are degenerated, Lambda is not diagonal anymore
lmb, R = la.eig(L_out)
# R is already computed in extract_out_modes, is it faster to keep it in memory?
# R is already computed in extract_out_modes, is it faster to keep
# it in memory?
Phi = Q @ R
N = Phi.conj().T @ Phi
......@@ -305,7 +321,7 @@ def compute_wf(vals, vecs, L_leads, X_out_leads, scat_dim,
ll = np.conj(l1) * l2
N[i, j] *= ll / (1 - ll)
for i, psi in enumerate(psi_alpha_0.T): #loop if degeneracy
for i, psi in enumerate(psi_alpha_0.T): # loop if degeneracy
a = q_e[:, i]
norm = psi.conj().T @ psi
norm += a.conj().T @ N @ a
......@@ -318,15 +334,18 @@ def compute_wf(vals, vecs, L_leads, X_out_leads, scat_dim,
else:
return psi_alpha_0, R @ q_e, np.diag(lmb), Phi
def max_norm(X):
"""
Norm every column of X such that the maximum value of the column is equal to 1
Norm every column of X such that the maximum value of the column is equal
to 1
"""
X_tilde = np.zeros(X.shape, dtype=complex)
arg_max, max_x = [], []
for i, x in enumerate(X.T):
x_max = np.argmax(abs(x)) # take the maximum to avoid division by small numbers
# take the maximum to avoid division by small numbers
x_max = np.argmax(abs(x))
X_tilde[:, i] = x / x[x_max]
arg_max.append(x_max)
max_x.append(x[x_max])
......
boundstate/fast.py
View file @ 403bff26
"""
Compute bound state using eigenvalue.
Use 'scan_eig_values' to have a guest of energies of bound states.
Bound stats are located where singular value = 0.
Use 'finder' to find the exact energy of bound states in a small
range of energies.
Use 'solver' to obtain the wave-function of the bound states at a
precise energy.
'solver' may return false positive.
"""
from __future__ import division
import numpy as np
import scipy.linalg as la
......@@ -5,9 +19,10 @@ from scipy.optimize import newton
np.set_printoptions(precision=16, suppress=1)
import scipy.sparse as sp
from ._common import *
from .standard import singular_values
from .standard import _singular_values
def BS_solver(H_s, H_leads, V_leads, transf, E, eps=1e-4, sparse=False):
def solver(syst, E, args=(), eps=1e-4, sparse=False):
"""
Raise an error if there is no bound state (singular value > eps)
Parameters:
......@@ -23,12 +38,37 @@ def BS_solver(H_s, H_leads, V_leads, transf, E, eps=1e-4, sparse=False):
Notes:
This function will be useless once the BS_solver will be included in kwant
"""
H_s, H_leads, V_leads, transf = extract_kwant_matrices(syst, args=args,
sparse=sparse)
return _solver(H_s, H_leads, V_leads, transf, E, eps, sparse)
vals, vecs, L_leads, X_out_leads = eig_values(H_s, H_leads, V_leads, transf, E,
sparse=sparse, uv=True)
def _solver(H_s, H_leads, V_leads, transf, E, eps=1e-4, sparse=False):
"""
Returns the elements to compute the wavefunction of a bound state
Raise an error if there is no bound state (singular value > eps)
Parameters:
transf: numpy array or sparse csr matrix. V.dot(transf) gives the
hopping matrix between the scattering region and the leads.
Returns:
The normalized wavefunction of the bound state,
psi_alpha_0: wavefunction in the scattering region
q_e: numpy array, as defined in Eq.7
L_out: Evanescent eigenvalues
Q: part of the wavefunction in the lead.
Notes:
This function will be useless once the BS_solver will be included in kwant
"""
vals, vecs, L_leads, X_out_leads = _eig_values(H_s, H_leads, V_leads,
transf, E, sparse=sparse,
uv=True)
N = H_s.shape[0]
return compute_wf(vals, vecs, L_leads, X_out_leads, N, eps=eps)
def basic_solver(H_0, V, E):
"""
Setup left and right hand side matrices from general eigenproblem
......@@ -59,6 +99,7 @@ def basic_solver(H_0, V, E):
return lmb_eva, X_full, (A, B)
def direct_method(A, B, A_dot, lmb, X):
"""
Compute the derivative of the eigenvectors X and eigenvalues
......@@ -83,7 +124,8 @@ def direct_method(A, B, A_dot, lmb, X):
X_dot = np.zeros(shape=X.shape, dtype=complex)
for i, (l, x) in enumerate(zip(lmb, X.T)):
lhs = np.hstack(((A - l * B), -np.dot(B, x)[:, None]))
lhs = np.delete(lhs, arg_max[i], axis=1) #remove column m because of norm
# remove column m because of norm
lhs = np.delete(lhs, arg_max[i], axis=1)
rhs = -np.dot(A_dot, x)
sol = la.solve(lhs, rhs)
......@@ -91,11 +133,12 @@ def direct_method(A, B, A_dot, lmb, X):
x_dot = sol[:-1]
x_dot *= max_x[i]
x_dot = x_dot.tolist()
x_dot.insert(arg_max[i], 0) #put back column m
X_dot[:, i] = x_dot #go back to initial norm
x_dot.insert(arg_max[i], 0) # put back column m
X_dot[:, i] = x_dot # go back to initial norm
return np.asarray(lmb_dot), X_dot
def derivative_lhs(E, H_s, H_leads, V_leads, transf, sparse=False):
"""
Compute the derivative of the lhs of Eq.(14) (H_eff in its hermitian
......@@ -115,7 +158,8 @@ def derivative_lhs(E, H_s, H_leads, V_leads, transf, sparse=False):
lmb_inv = 1 / np.diag(lmb)
lmb_dot_inv, X_dot_full = direct_method(A, B, A_dot, lmb_inv, X_full)
lmb_dot = -1 / lmb_inv**2 * lmb_dot_inv #because original problem formulated in 1/lmb
# because original problem formulated in 1/lmb
lmb_dot = -1 / lmb_inv**2 * lmb_dot_inv
X_dot.append(X_dot_full[:N_l, :])
X.append(X_full[:N_l, :])
......@@ -156,24 +200,27 @@ def derivative_lhs(E, H_s, H_leads, V_leads, transf, sparse=False):
H_prime[N_s:, N_s:] -= L.conj().dot(X.conj().T).dot(V).dot(X_dot)
return H_prime
def eig_val_derivative(E, H_s, H_leads, V_leads, transf, sparse=False,
uv=True, need_orth=True, eig_val=2):
"""
Compute the derivative of the smallest eigenvalue of H_eff wrt to the
energy, Eq.20., de_alpha / dE = psi_dag dH_eff / dE psi.
"""
if not np.isscalar(E): # because fmin_tnc uses array and not scalars
if not np.isscalar(E): # because fmin_tnc uses array and not scalars
E = E[0]
L_out_leads, X_out_leads = leads_modes(H_leads, V_leads, E, need_orth=False)
L_out_leads, X_out_leads = leads_modes(H_leads, V_leads, E,
need_orth=False)
H = setup_lhs_H(H_s, V_leads, transf, L_out_leads, X_out_leads, E,
sparse=sparse)
sparse=sparse)
# if eig_val_derivative is called from a root finder, then setup_lhs_H()
# has already been called in eig_values() at the same energy. Should be fixed...
# has already been called in eig_values() at the same energy.
# Should be fixed...
if sparse:
vals, vecs = sp.linalg.eigsh(H, k=3, sigma=0)
min_val = np.argmin(abs(vals)) # find the closest eig value to 0
min_val = np.argmin(abs(vals)) # find the closest eig value to 0
vec = vecs[:, min_val]
H_prime = derivative_lhs(E, H_s, H_leads, V_leads,
transf, sparse=sparse)
......@@ -182,37 +229,56 @@ def eig_val_derivative(E, H_s, H_leads, V_leads, transf, sparse=False,
else:
# compute all eigenvalues
vals, vecs = la.eigh(H)
min_val = np.argmin(abs(vals)) # find the closest eig value to 0
min_val = np.argmin(abs(vals)) # find the closest eig value to 0
vec = vecs[:, min_val]
#~ vec /= np.sqrt(np.dot(vec.conj().T, vec)) # already normalized to 1
# ~ vec /= np.sqrt(np.dot(vec.conj().T, vec)) # already normalized to 1
H_prime = derivative_lhs(E, H_s, H_leads, V_leads,
transf, sparse=sparse)
return vec.conj().T.dot(H_prime).dot(vec) #simple expression because of hermitian formulation
# Simple expression because of hermitian formulation
return vec.conj().T.dot(H_prime).dot(vec)
def BS_solver(H_s, H_leads, V_leads, transf, E, tol=1e-4, sparse=False):
def eig_values(syst, E, args=(), sparse=False, uv=True,
need_orth=True, sing_values=4, tol=1e-5):
"""
Returns the elements to compute the wavefunction of a bound state
Raise an error if there is no bound state (singular value > tol)
Solver that returns the singular values of the SVD of
the left hand side of the bound state equation (L in the notes)
Parameters:
transf: numpy array or sparse csr matrix. V.dot(transf) gives the
hopping matrix between the scattering region and the leads.
H_s: Hamiltonian of the central system
H_leads: tuple,
countains the onsite matrices of the unit cells of the leads
V_leads: tuple, countains the hopping matrices between the unit
cells of every lead
uv: Whether to compute the the full SVD of the matrix L
or only the singular values
eig_val: the number of eigenvalues to be computed. Only if sparse
is True
returns:
if uv is false: the singular values of the matrix L
else:
S: The singular values of L
Wh: The right matrix of the svd decomposition
L: array that countains the evanescent outgoing lambdas
X_out: Columns are the corresponding evanescent wavefunctions
"""
H_s, H_leads, V_leads, transf = extract_kwant_matrices(syst, args=args,
sparse=sparse)
return _eig_values(H_s, H_leads, V_leads, transf, E, sparse, uv,
need_orth, sing_values, tol)
vals, vecs, L_leads, X_out_leads = eig_values(H_s, H_leads, V_leads, transf, E,
sparse=sparse, uv=True)
N = H_s.shape[0]
return compute_wf(vals, vecs, L_leads, X_out_leads, N, eps=tol)
def eig_values(H_s, H_leads, V_leads, transf, E, sparse=False, uv=True,
need_orth=True, sing_values=4, tol=1e-5):
def _eig_values(H_s, H_leads, V_leads, transf, E, sparse=False, uv=True,
need_orth=True, sing_values=4, tol=1e-5):
"""
Solver that returns the singular values of the SVD of
the left hand side of the bound state equation (L in the notes)
Parameters:
H_s: Hamiltonian of the central system
H_leads: tuple, countains the onsite matrices of the unit cells of the leads
H_leads: tuple,
countains the onsite matrices of the unit cells of the leads
V_leads: tuple, countains the hopping matrices between the unit
cells of every lead
uv: Whether to compute the the full SVD of the matrix L
......@@ -231,7 +297,7 @@ def eig_values(H_s, H_leads, V_leads, transf, E, sparse=False, uv=True,
L_out_leads, X_out_leads = leads_modes(H_leads, V_leads, E, tol=tol)
lhs = setup_lhs_H(H_s, V_leads, transf, L_out_leads, X_out_leads, E,
sparse=sparse)
sparse=sparse)
if sparse:
if not uv:
......@@ -249,7 +315,9 @@ def eig_values(H_s, H_leads, V_leads, transf, E, sparse=False, uv=True,
vals, vecs = la.eigh(lhs)
return vals, vecs, L_out_leads, X_out_leads
def setup_lhs_H(H_s, V_leads, transf, L_out_leads, X_out_leads, E, sparse=False):
def setup_lhs_H(H_s, V_leads, transf, L_out_leads, X_out_leads, E,
sparse=False):
"""
Setup the left hand side of the equation to find bound states.
Parameters:
......@@ -264,7 +332,7 @@ def setup_lhs_H(H_s, V_leads, transf, L_out_leads, X_out_leads, E, sparse=False)
X_out = block_diag(*X_out_leads)
V = block_diag(*V_leads)
if sparse: #Uses either scipy.sparse if sparse of numpy if dense
if sparse: # Uses either scipy.sparse if sparse of numpy if dense
d = sp.csr_matrix.dot
lin = sp
V = sp.csr_matrix(V)
......@@ -282,59 +350,50 @@ def setup_lhs_H(H_s, V_leads, transf, L_out_leads, X_out_leads, E, sparse=False)
lhs = lin.hstack([H_s - E * Id_N, d(V_ls.conj().T.tocsr(), X_L)])
else:
L_X = X_L.conj().T
lhs = lin.hstack([H_s - E * Id_N, d(V_ls.conj().T, X_L)])
lhs = lin.vstack([lhs, lin.hstack([d(L_X, V_ls), -d(L_X, d(V, X_out))])])
return lhs
def plot_eig_values(sys, e_min, e_max, N=500, args=(), sparse=False, k=5):
energies = np.linspace(e_min, e_max, N)
H_s, H_leads, V_leads, transf = extract_kwant_matrices(sys, args=args, sparse=sparse)
if sparse:
dumb = np.zeros(shape=(N, k))
else:
# at maximum the number of singular value is the dimension of H_s
# + the dimension of all the leads together
dumb = np.zeros(shape=(N, H_s.shape[0] + sum(V.shape[0] for V in V_leads)))
# the zeros are later replaced with NaN so they are not plotted.
for i, e in enumerate(energies):
s = eig_values(H_s, H_leads, V_leads, transf, e, sparse,
uv=False, sing_values=k)
def finder(syst, E_0, args=(), sparse=False, tol=1e-6, sing_values=4,
uv=False, fprime=eig_val_derivative):
"""
Function that search for a minimum of the eigenvalues of the lhs of Eq.13
Uses the newton algorithm.
dumb[i, :len(s)] = s
dumb = fill_zero_with_nan(dumb)
Parameters:
E_0 : float. An initial estimate of the bound stat energy.
fprime: None or function (same arguments for eig_val and
eig_val_derivative) that returns the derivative of an
eigenvalue of Eq.13
return (energies, dumb)
Returns:
E: The energy of the bound state if it found one or print a message
def false_positive_BS(H_s, H_leads, V_leads, transf, E, eps=1e-8,
sparse=False):
"""
Check the False positive computing the singular values of H_eff in
it non-hermitian form.
It would be also possible to check Eq.19, something like:
np.amax(block_diag(V_leads) @ (P @ psi_alpha_0 - phi_e @ q_e)) < eps
Important note:
If there are degenerated eigenvalues, then fprime should be equal to
None, as I have no clue how to compute the derivative of fully degenerated
eigenvalues (or the derivative of an invariant subspace).
Ideas welcome...
"""
sing = singular_values(H_s, H_leads, V_leads, transf, E, uv=False,
sparse=False)
if np.amin(sing) < eps:
return False # True bound state
else:
return True # True means a true "false positive" (i.e., no bound states)
H_s, H_leads, V_leads, transf = extract_kwant_matrices(syst, args=args,
sparse=sparse)
return _finder(H_s, H_leads, V_leads, transf, E_0, sparse,
tol, sing_values, uv, fprime)
def BS_finder_derivative(H_s, H_leads, V_leads, transf, E_0, sparse=False,
tol=1e-6, sing_values=4, uv=False,
fprime=eig_val_derivative):
def _finder(H_s, H_leads, V_leads, transf, E_0, sparse=False,
tol=1e-6, sing_values=4, uv=False, fprime=eig_val_derivative):
"""
Function that search for a minimum of the eigenvalues of the lhs of Eq.13
Uses the newton algorithm.
Parameters:
fprime: None or function (same arguments for eig_val and eig_val_derivative)
that returns the derivative of an eigenvalue of Eq.13
fprime: None or function (same arguments for eig_val and
eig_val_derivative) that returns the derivative of an
eigenvalue of Eq.13
Returns:
E: The energy of the bound state if it found one or print a message
......@@ -347,15 +406,15 @@ def BS_finder_derivative(H_s, H_leads, V_leads, transf, E_0, sparse=False,
"""
def dumb(e, H_s, H_leads, V_leads, transf, sparse=sparse):
eig = eig_values(H_s, H_leads, V_leads, transf, e, sparse=sparse,
uv=False, sing_values=sing_values)
eig = _eig_values(H_s, H_leads, V_leads, transf, e, sparse=sparse,
uv=False, sing_values=sing_values)
idx = np.argmin(abs(eig))
return eig[idx]
try:
E = newton(dumb, E_0, fprime=fprime, tol=tol,
args=(H_s, H_leads, V_leads, transf, sparse))
if false_positive_BS(H_s, H_leads, V_leads, transf, E, tol=tol):
if false_positive_BS(H_s, H_leads, V_leads, transf, E, eps=tol):
print('False positive, no bound state')
return False
else:
......@@ -363,3 +422,56 @@ def BS_finder_derivative(H_s, H_leads, V_leads, transf, E_0, sparse=False,
return E
except RuntimeError:
print('no bound states found')
def scan_eig_values(syst, e_min, e_max, args=(), N=500, sparse=False, k=5):
energies = np.linspace(e_min, e_max, N)
"""Compute the eigenvalues on a energy range.
See 'eig_values' for details.
"""
H_s, H_leads, V_leads, transf = extract_kwant_matrices(syst, args=args,
sparse=sparse)
return _scan_eig_values(H_s, H_leads, V_leads, transf, e_min, e_max, N,
sparse, k)
def _scan_eig_values(H_s, H_leads, V_leads, transf, e_min, e_max, N=500,
sparse=False, k=5):
energies = np.linspace(e_min, e_max, N)
"""Compute the eigenvalues on a energy range.
See 'eig_values' for details.
"""
if sparse:
dumb = np.zeros(shape=(N, k))
else:
# at maximum the number of singular value is the dimension of H_s
# + the dimension of all the leads together
dumb = np.zeros(shape=(N, H_s.shape[0] + sum(V.shape[0]
for V in V_leads)))
# the zeros are later replaced with NaN so they are not plotted.
for i, e in enumerate(energies):
s = _eig_values(H_s, H_leads, V_leads, transf, e, sparse,
uv=False, sing_values=k)
dumb[i, :len(s)] = s
dumb = fill_zero_with_nan(dumb)