Commit aae2de5e authored by Mathieu Istas's avatar Mathieu Istas Committed by Christoph Groth
Browse files

add bound state code

parents
==================
Boundstate license
==================
Copyright 2016-2019 M. Istas (CEA), B. Rossignol (CEA), C. Groth (CEA),
X. Waintal (CEA), A. Akhmerov (TUD), and others.
All rights reserved.
(CEA = Commissariat à l'énergie atomique et aux énergies alternatives,
TUD = Technische Universiteit Delft)
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
This package implements an algorithm for computing bound states of infinite
tight-binding systems that are made up of a finite scattering region connected
to semi-infinite leads. The details of the algorithms are explained in an
`accompanying article <dx.doi.org/10.21468/SciPostPhys.4.5.026>`_.
This package is free software under the `BSD License <LICENSE.rst>`_. If you
use it for work that leads to a scientific publication, please cite the above
article.
from __future__ import division
import numpy as np
import scipy.linalg as la
from kwant.physics.leads import unified_eigenproblem, setup_linsys
from scipy.optimize import minimize_scalar
np.set_printoptions(precision=16, suppress=1)
import scipy.sparse as sp
import scipy
def extract_out_modes(H_0, V, E, eps=1e-6):
"""
parameters:
H_0: hamiltonian inside the unite cell of the lead
V: hoppings between the cells
E: Energy at which the mode is calculated
returns:
modes: a list of evanescent right going modes
"""
n = H_0.shape[0]
Id = np.eye(n)
H = H_0 - E * Id
V_dag = V.conj().T
matrices, v, extract = setup_linsys(H, V)
all_lmb, evanselect, propselect, vec_gen, ord_schur = unified_eigenproblem(*(matrices))
lmb_inv = all_lmb[evanselect]
full_rank = len(all_lmb) // 2 == n
all_vecs = vec_gen(evanselect)
X_out = np.empty(shape=(n, len(lmb_inv)), dtype=complex)
for i, l in enumerate(lmb_inv):
phi = extract(all_vecs[:, i], l)
assert(np.allclose(np.dot(H + V * l + V_dag / l, phi), np.zeros(H.shape[0])))
X_out[:, i] = phi
return 1 / lmb_inv, X_out, full_rank
def singular_values(H_s, H, V, V_ls, DBh, DBh_ls, E, sparse,
uv=True, sing_values=4, eps=1e-4, sigma=None):
"""
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
D: Singular values of the hopping
Bh: right matrix of the SVD of the hopping, such that
A D Bh = V
V: hopping matrix between the cells of the lead
V_ls: hopping matrix between the scattering region and the lead(s)
uv: Whether to compute the the full SVD of the matrix L
or only the singular values
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
"""
small_num = np.finfo(float).eps * 1e2
if sparse: #Uses either scipy.sparse if sparse of numpy if dense
d = sp.csr_matrix.dot
lin = sp
else:
d = np.dot
lin = np
Id_N = lin.eye(H_s.shape[0])
X_out = []
L = []
for i, lead in enumerate(sys.leads):
L_lead, X_out_lead, full_rank = extract_out_modes(H, V, E, eps=eps)
X_out.append(X_out_lead)
L.append(L_lead)
l = len(L)
L_out = np.diag(L)
V_ls_dag = V_ls.conj().T
if sparse:
V_ls_dag = V_ls_dag.tocsr()
V = sp.csr_matrix(V)
L_out = sp.csr_matrix(L_out)
X_out = sp.csr_matrix(X_out)
DBh = sp.csr_matrix(DBh)
DBh_ls = sp.csr_matrix(DBh_ls)
# l is the number of out evanescent states, while r correpond to
# the rank of V, so half of the non trivial solutions to the mode equation
lhs = lin.hstack([H_s - E * Id_N, d(V_ls_dag, d(X_out, L_out))])
if full_rank:
lhs = lin.vstack([lhs, lin.hstack([V_ls, -d(V, X_out)])])
else:
lhs = lin.vstack([lhs, lin.hstack([DBh_ls, -d(DBh, X_out)])])
if sparse:
B = sp.coo_matrix.dot(lhs.conj().T, lhs)
if sigma is None:
# First compute the biggest eigenvalue
max_eval = sp.linalg.eigsh(B, k=1, which='LM', return_eigenvectors=False)[0]
# shift the matrix so that the lowest eigenvalue become the biggest in magnitude
B = -B + 2 * max_eval * sp.eye(B.shape[0])
if not uv:
try: # if the singular value is exactly zero
evals = sp.linalg.eigsh(B, return_eigenvectors=False,
k=sing_values, sigma=sigma)
except RuntimeError:
evals, evecs = sp.linalg.eigsh(B + eps*sp.eye(B.shape[0]),
return_eigenvectors=True,
k=sing_values, sigma=sigma)
evals -= eps
if sigma is None: # that is, use the shift method
return np.sqrt(2 * max_eval - evals)
else: # use of the shift invert method
return np.sqrt(evals)
else:
try:
evals, evecs = sp.linalg.eigsh(B, return_eigenvectors=True,
k=sing_values, sigma=sigma)
except RuntimeError:
evals, evecs = sp.linalg.eigsh(B + eps*sp.eye(B.shape[0]),
return_eigenvectors=True,
k=sing_values, sigma=sigma)
evals -= eps
if sigma is None: # shift method
return np.sqrt(2 * max_eval - evals), evecs.T, L, X_out
else: # shift invert method
# If an eigenvalue is negative due to numerical
# precision, then it is replaced by 0
evals = np.clip(evals, 0, np.inf)
return np.sqrt(evals), evecs.T, L, X_out
else:
if not uv:
return la.svd(lhs, compute_uv=False)
else:
U, S, Wh = la.svd(lhs, compute_uv=True)
return S, Wh, L, X_out
def BS_solver(H_S, H_0, V, transf, E, eps=1e-4, sparse=False, sigma=None):
"""
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.
"""
H_S, DBh = setup_matrices(H_S, V, sparse)
if sparse:
V_ls = sp.csr_matrix(V).dot(transf)
DBh = sp.csr_matrix(DBh)
DBh_ls = DBh.dot(transf)
else:
if type(transf).__module__ is scipy:
transf = sp.spmatrix.toarray(transf)
V_ls = V.dot(transf)
DBh_ls = DBh.dot(transf)
S, Wh, L, X_out = singular_values(H_S, H_0, V, V_ls, DBh,
DBh_ls, E, sparse, uv=True, sigma=sigma)
assert(min(S) < eps), 'The energy does not match the bound state energy, {0}'.format(min(S))
N = H_S.shape[0]
return compute_wf(S, Wh, L, X_out, N, sparse=sparse, eps=eps)
def compute_wf(S, Wh, L, X_out, N, sparse=False, eps=1e-4):
W = Wh.conj().T
zero_schur = S < eps
# the wave function is given by the columns of W
# Psi_alpha_0 is the wavefunction of the bound states in the
# system.
psi_alpha_0, a_alpha = W[:N, zero_schur], W[N:, zero_schur]
#### normalization
if sparse:
X_out = sp.csr_matrix.toarray(X_out)
M = np.dot(X_out.conj().T, X_out)
for i, l1 in enumerate(L):
for j, l2 in enumerate(L):
ll = np.conj(l1) * l2
M[i, j] *= ll / (1 - ll)
for i, psi in enumerate(psi_alpha_0.T): #loop if degeneracy
a = a_alpha[:, i]
norm = np.dot(psi.conj().T, psi)
norm += np.dot(a.conj().T, np.dot(M, a))
norm = np.sqrt(norm)
psi_alpha_0[:, i] /= norm
a_alpha[:, i] /= norm
return psi_alpha_0, a_alpha, L, X_out
def BS_finder(H_S, H_0, V, transf, e_range, sparse=False, eps=1e-6, sigma=None, tol=1e-6):
"""
Finds the minimum of the singular values of L in e_range.
If the singular value is smaller than eps, then it returns
the energy of the bound state.
Parameters:
e_range: tuple countaining the minimum and the maximum of the
energy range
returns:
E: energy of the bound state if minimal singular value is smaller
than eps
None otherwise
"""
H_S, DBh = setup_matrices(H_S, V, sparse, eps=eps)
if sparse:
V_ls = sp.csr_matrix(V).dot(transf)
DBh = sp.csr_matrix(DBh)
DBh_ls = DBh.dot(transf)
else:
try:
transf = sp.spmatrix.toarray(transf)
except TypeError:
pass
V_ls = V.dot(transf)
DBh_ls = DBh.dot(transf)
def dumb(e):
s = min(singular_values(H_S, H_0, V, V_ls, DBh,
DBh_ls, e, sparse, uv=False, sigma=sigma))
return s
minimum = minimize_scalar(dumb, method='Bounded', bounds=e_range, tol=tol)
if minimum.fun < eps:
# if the last singular value is smaller than eps,
# then there is a bound state.
E = minimum.x
return E
else:
return None
def setup_matrices(H_S, V, sparse, eps=1e-4):
_, D, Bh = la.svd(V)
# keep only the vectors assiocated with non zero singular value
non_zeros = sum(D > eps)
D = np.diag(D[:non_zeros])
Bh = Bh[:non_zeros, :]
DBh = np.dot(D, Bh)
if sparse:
H_S = H_S.tocsr()
return H_S, DBh
def extract_kwant_matrices(sys, sparse=False):
"""
Simple function that extract the matrices from a kwant system
"""
H_S = sys.hamiltonian_submatrix(sparse=sparse)
H_0 = sys.leads[0].cell_hamiltonian()
V = sys.leads[0].inter_cell_hopping()
def P_matrix(sys):
lhs, norb = sys.hamiltonian_submatrix(sparse=sparse, return_norb=True)[:2]
offsets = np.empty(norb.shape[0] + 1, int)
offsets[0] = 0
offsets[1 :] = np.cumsum(norb)
interface = sys.lead_interfaces[0]
iface_orbs = np.r_[tuple(slice(offsets[i], offsets[i + 1])
for i in interface)]
coords = np.r_[[np.arange(len(iface_orbs))], [iface_orbs]]
transf = sp.csr_matrix((np.ones(len(iface_orbs)), coords),
shape=(iface_orbs.size, lhs.shape[0]))
return transf
transf = P_matrix(sys)
return H_S, H_0, V, transf
def extract_kwant_matrices_v_2(sys, sparse=False):
#Simple function that extract the matrices from a kwant system
H_S = sys.hamiltonian_submatrix(sparse=sparse)
H_0s = []
Vs = []
coords = []
for i, lead in enumerate(sys.leads):
H_0s.append(sys.leads[i].cell_hamiltonian())
Vs.append(sys.leads[i].inter_cell_hopping())
coords = np.concatenate((coords, sys.lead_interfaces[i]))
H_0 = la.block_diag(*H_0s)
V = la.block_diag(*Vs)
norb = sys.hamiltonian_submatrix(sparse=sparse, return_norb=True)[1]
offsets = np.empty(norb.shape[0] + 1, int)
offsets[0] = 0
offsets[1 :] = np.cumsum(norb)
ones = []
for coord in coords:
coord = int(coord)
for i in range(offsets[coord], offsets[coord] + np.array(norb)[coord]):
ones.append(i)
transf = sp.csr_matrix((np.ones(len(ones)), (range(len(ones)),ones)),
shape=(len(ones),sum(norb)))
return H_S, H_0, V, transf
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment