Coverage for kwant/physics/leads.py : 97%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Copyright 2011-2016 Kwant authors.
2#
3# This file is part of Kwant. It is subject to the license terms in the file
4# LICENSE.rst found in the top-level directory of this distribution and at
5# https://kwant-project.org/license. A list of Kwant authors can be found in
6# the file AUTHORS.rst at the top-level directory of this distribution and at
7# https://kwant-project.org/authors.
10from math import sin, cos, sqrt, pi, copysign
11from collections import namedtuple
13from itertools import combinations_with_replacement
14import numpy as np
15import numpy.linalg as npl
16import scipy.linalg as la
17from .. import linalg as kla
18from scipy.linalg import block_diag
19from scipy.sparse import (identity as sp_identity, hstack as sp_hstack,
20 csr_matrix)
22dot = np.dot
24__all__ = ['selfenergy', 'modes', 'PropagatingModes', 'StabilizedModes']
27def nonzero_symm_projection(matrix):
28 """Check whether a discrete symmetry relation between two blocks of the
29 Hamiltonian vanishes or not.
31 For a discrete symmetry S, the projection that maps from the block i to
32 block j of the Hamiltonian with projectors P_i and P_j is
33 S_ji = P_j^+ S P_i.
34 This function determines whether S_ji vanishes or not, i. e. whether a
35 symmetry relation exists between blocks i and j.
37 Here, discrete symmetries and projectors are in canonical form, such that
38 S maps each block at most either to itself or to a single other block.
39 In other words, for each j, the block S_ji is nonzero at most for one i,
40 while all other blocks vanish. If a nonzero block exists, it is unitary
41 and hence has norm 1. To identify nonzero blocks without worrying about
42 numerical errors, we thus check that its norm is larger than 0.5.
43 """
44 if not isinstance(matrix, np.ndarray):
45 matrix = matrix.data
46 return np.linalg.norm(matrix) > 0.5
49def group_halves(arr_list):
50 """Split and rearrange a list of arrays.
52 Combine a list of arrays into a single array where the first halves
53 of each array appear first:
54 `[a b], [c d], [e f] -> [a c e b d f]`
55 """
56 list_ = [np.split(arr, 2) for arr in arr_list]
57 lefts, rights = zip(*list_)
58 return np.r_[tuple(lefts + rights)]
61# Container classes
62Linsys = namedtuple('Linsys', ['eigenproblem', 'v', 'extract'])
65class PropagatingModes:
66 """The calculated propagating modes of a lead.
68 Attributes
69 ----------
70 wave_functions : numpy array
71 The wave functions of the propagating modes.
72 momenta : numpy array
73 Momenta of the modes.
74 velocities : numpy array
75 Velocities of the modes.
76 block_nmodes: list of integers
77 Number of left or right moving propagating modes
78 per conservation law block of the Hamiltonian.
80 Notes
81 =====
82 The sort order of all the three arrays is identical. The first half of the
83 modes have negative velocity, the second half have positive velocity.
84 Within these halves the modes are ordered by the eigenvalues of any
85 declared conservation laws. Within blocks with the same conservation law
86 eigenvalue the modes with negative velocity are ordered by increasing
87 momentum, and the modes with positive velocity are ordered by decreasing
88 momentum. Finally, modes are ordered by the magnitude of their velocity.
89 To summarize, the modes are ordered according to the key
90 `(sign(v), conserved_quantity, sign(v) * k , abs(v))` where `v` is
91 velocity, `k` is momentum and `conserved_quantity` is the conservation
92 law eigenvalue.
94 In the above, the positive velocity and momentum directions are defined
95 with respect to the translational symmetry direction of the system.
97 The first dimension of `wave_functions` corresponds to the orbitals of all
98 the sites in a unit cell, the second one to the number of the mode. Each
99 mode is normalized to carry unit current. If several modes have the same
100 momentum and velocity, an arbitrary orthonormal basis in the subspace of
101 these modes is chosen.
103 If a conservation law is specified to block diagonalize the Hamiltonian,
104 then `block_nmodes[i]` is the number of left or right moving propagating
105 modes in conservation law block `i`. The ordering of blocks is the same as
106 the ordering of the projectors used to block diagonalize the Hamiltonian.
107 """
108 def __init__(self, wave_functions, velocities, momenta):
109 kwargs = locals()
110 kwargs.pop('self')
111 self.__dict__.update(kwargs)
114class StabilizedModes:
115 """Stabilized eigendecomposition of the translation operator.
117 Due to the lack of Hermiticity of the translation operator, its
118 eigendecomposition is frequently poorly conditioned. Solvers in Kwant use
119 this stabilized decomposition of the propagating and evanescent modes in
120 the leads. If the hopping between the unit cells of an infinite system is
121 invertible, the translation eigenproblem is written in the basis `psi_n,
122 h_hop^+ * psi_(n+1)`, with ``h_hop`` the hopping between unit cells. If
123 `h_hop` is not invertible, and has the singular value decomposition `u s
124 v^+`, then the eigenproblem is written in the basis `sqrt(s) v^+ psi_n,
125 sqrt(s) u^+ psi_(n+1)`. In this basis we calculate the eigenvectors of the
126 propagating modes, and the Schur vectors (an orthogonal basis) in the space
127 of evanescent modes.
129 `vecs` and `vecslmbdainv` are the first and the second halves of the wave
130 functions. The first `nmodes` are eigenmodes moving in the negative
131 direction (hence they are incoming into the system in Kwant convention),
132 the second `nmodes` are eigenmodes moving in the positive direction. The
133 remaining modes are the Schur vectors of the modes evanescent in the
134 positive direction. Propagating modes with the same eigenvalue are
135 orthogonalized, and all the propagating modes are normalized to carry unit
136 current. Finally the `sqrt_hop` attribute is `v sqrt(s)`.
138 Attributes
139 ----------
140 vecs : numpy array
141 Translation eigenvectors.
142 vecslmbdainv : numpy array
143 Translation eigenvectors divided by the corresponding eigenvalues.
144 nmodes : int
145 Number of left-moving (or right-moving) modes.
146 sqrt_hop : numpy array
147 Part of the SVD of `h_hop`.
148 """
150 def __init__(self, vecs, vecslmbdainv, nmodes, sqrt_hop):
151 kwargs = locals()
152 kwargs.pop('self')
153 self.__dict__.update(kwargs)
155 def selfenergy(self):
156 """
157 Compute the self-energy generated by lead modes.
159 Returns
160 -------
161 Sigma : numpy array, real or complex, shape (M,M)
162 The computed self-energy. Note that even if `h_cell` and `h_hop` are
163 both real, `Sigma` will typically be complex. (More precisely, if
164 there is a propagating mode, `Sigma` will definitely be complex.)
165 """
166 v = self.sqrt_hop
167 # Only include the *outgoing* modes (propagating + evanescent),
168 # meaning this is the *retarded* self-energy.
169 outgoing = slice(self.nmodes, None)
170 vecs = self.vecs[:, outgoing]
171 vecslmbdainv = self.vecslmbdainv[:, outgoing]
172 return dot(v, dot(vecs, la.solve(vecslmbdainv, v.T.conj())))
175# Auxiliary functions that perform different parts of the calculation.
176def setup_linsys(h_cell, h_hop, tol=1e6, stabilization=None):
177 """Make an eigenvalue problem for eigenvectors of translation operator.
179 Parameters
180 ----------
181 h_cell : numpy array with shape (n, n)
182 Hamiltonian of a single lead unit cell.
183 h_hop : numpy array with shape (n, m), m <= n
184 Hopping Hamiltonian from a cell to the next one.
185 tol : float
186 Numbers are considered zero when they are smaller than `tol` times
187 the machine precision.
188 stabilization : sequence of 2 booleans or None
189 Which steps of the eigenvalue problem stabilization to perform. If the
190 value is `None`, then Kwant chooses the fastest (and least stable)
191 algorithm that is expected to be sufficient. For any other value,
192 Kwant forms the eigenvalue problem in the basis of the hopping singular
193 values. The first element set to `True` forces Kwant to add an
194 anti-Hermitian term to the cell Hamiltonian before inverting. If it is
195 set to `False`, the extra term will only be added if the cell
196 Hamiltonian isn't invertible. The second element set to `True` forces
197 Kwant to solve a generalized eigenvalue problem, and not to reduce it
198 to the regular one. If it is `False`, reduction to a regular problem
199 is performed if possible.
201 Returns
202 -------
203 linsys : namedtuple
204 A named tuple containing `matrices` a matrix pencil defining
205 the eigenproblem, `v` a hermitian conjugate of the last matrix in
206 the hopping singular value decomposition, and functions for
207 extracting the wave function in the unit cell from the wave function
208 in the basis of the nonzero singular exponents of the hopping.
210 Notes
211 -----
212 The lead problem with degenerate hopping is rather complicated, and the
213 details of the algorithm will be published elsewhere.
215 """
216 n = h_cell.shape[0]
217 m = h_hop.shape[1]
218 if stabilization is not None:
219 stabilization = list(stabilization)
221 if not np.any(h_hop): # skip coverage
222 # Inter-cell hopping is zero. The current algorithm is not suited to
223 # treat this extremely singular case.
224 raise ValueError("Inter-cell hopping is exactly zero.")
226 # If both h and t are real, it may be possible to use the real eigenproblem.
227 if (not np.any(h_hop.imag)) and (not np.any(h_cell.imag)):
228 h_hop = h_hop.real
229 h_cell = h_cell.real
231 eps = np.finfo(np.common_type(h_cell, h_hop)).eps * tol
233 # First check if the hopping matrix has singular values close to 0.
234 # (Close to zero is defined here as |x| < eps * tol * s[0] , where
235 # s[0] is the largest singular value.)
237 u, s, vh = la.svd(h_hop)
238 assert m == vh.shape[1], "Corrupt output of svd."
239 n_nonsing = np.sum(s > eps * s[0])
241 if (n_nonsing == n and stabilization is None):
242 # The hopping matrix is well-conditioned and can be safely inverted.
243 # Hence the regular transfer matrix may be used.
244 h_hop_sqrt = sqrt(np.linalg.norm(h_hop))
245 A = h_hop / h_hop_sqrt
246 B = h_hop_sqrt
247 B_H_inv = 1.0 / B # just a real scalar here
248 A_inv = la.inv(A)
250 lhs = np.zeros((2*n, 2*n), dtype=np.common_type(h_cell, h_hop))
251 lhs[:n, :n] = -dot(A_inv, h_cell) * B_H_inv
252 lhs[:n, n:] = -A_inv * B
253 lhs[n:, :n] = A.T.conj() * B_H_inv
255 def extract_wf(psi, lmbdainv):
256 return B_H_inv * np.copy(psi[:n])
258 matrices = (lhs, None)
259 v_out = h_hop_sqrt * np.eye(n)
260 else:
261 if stabilization is None:
262 stabilization = [None, False]
264 # The hopping matrix has eigenvalues close to 0 - those
265 # need to be eliminated.
267 # Recast the svd of h_hop = u s v^dagger such that
268 # u, v are matrices with shape n x n_nonsing.
269 u = u[:, :n_nonsing]
270 s = s[:n_nonsing]
271 u = u * np.sqrt(s)
272 # pad v with zeros if necessary
273 v = np.zeros((n, n_nonsing), dtype=vh.dtype)
274 v[:vh.shape[1]] = vh[:n_nonsing].T.conj()
275 v = v * np.sqrt(s)
277 # Eliminating the zero eigenvalues requires inverting the on-site
278 # Hamiltonian, possibly including a self-energy-like term. The
279 # self-energy-like term stabilizes the inversion, but the most stable
280 # choice is inherently complex. This can be disadvantageous if the
281 # Hamiltonian is real, since staying in real arithmetics can be
282 # significantly faster. The strategy here is to add a complex
283 # self-energy-like term always if the original Hamiltonian is complex,
284 # and check for invertibility first if it is real
286 matrices_real = issubclass(np.common_type(h_cell, h_hop), np.floating)
287 add_imaginary = stabilization[0] or ((stabilization[0] is None) and
288 not matrices_real)
289 # Check if there is a chance we will not need to add an imaginary term.
290 if not add_imaginary:
291 h = h_cell
292 sol = kla.lu_factor(h)
293 rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))
295 if rcond < eps:
296 need_to_stabilize = True
297 else:
298 need_to_stabilize = False
300 if add_imaginary or need_to_stabilize:
301 need_to_stabilize = True
302 # Matrices are complex or need self-energy-like term to be
303 # stabilized.
304 temp = dot(u, u.T.conj()) + dot(v, v.T.conj())
305 h = h_cell + 1j * temp
307 sol = kla.lu_factor(h)
308 rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))
310 # If the condition number of the stabilized h is
311 # still bad, there is nothing we can do.
312 if rcond < eps: 312 ↛ 313line 312 didn't jump to line 313, because the condition on line 312 was never true
313 raise RuntimeError("Flat band encountered at the requested "
314 "energy, result is badly defined.")
316 # The function that can extract the full wave function psi from
317 # the projected one (v^dagger psi lambda^-1, s u^dagger psi).
319 def extract_wf(psi, lmbdainv):
320 wf = -dot(u, psi[: n_nonsing] * lmbdainv) - dot(v, psi[n_nonsing:])
321 if need_to_stabilize:
322 wf += 1j * (dot(v, psi[: n_nonsing]) +
323 dot(u, psi[n_nonsing:] * lmbdainv))
324 return kla.lu_solve(sol, wf)
326 # Setup the generalized eigenvalue problem.
328 A = np.zeros((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
329 B = np.zeros((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
331 begin, end = slice(n_nonsing), slice(n_nonsing, None)
333 A[end, begin] = np.identity(n_nonsing)
334 temp = kla.lu_solve(sol, v)
335 temp2 = dot(u.T.conj(), temp)
336 if need_to_stabilize:
337 A[begin, begin] = -1j * temp2
338 A[begin, end] = temp2
339 temp2 = dot(v.T.conj(), temp)
340 if need_to_stabilize:
341 A[end, begin] -= 1j *temp2
342 A[end, end] = temp2
344 B[begin, end] = -np.identity(n_nonsing)
345 temp = kla.lu_solve(sol, u)
346 temp2 = dot(u.T.conj(), temp)
347 B[begin, begin] = -temp2
348 if need_to_stabilize:
349 B[begin, end] += 1j * temp2
350 temp2 = dot(v.T.conj(), temp)
351 B[end, begin] = -temp2
352 if need_to_stabilize:
353 B[end, end] = 1j * temp2
355 v_out = v[:m]
357 # Solving a generalized eigenproblem is about twice as expensive
358 # as solving a regular eigenvalue problem.
359 # Computing the LU factorization is negligible compared to both
360 # (approximately 1/30th of a regular eigenvalue problem).
361 # Because of this, it makes sense to try to reduce
362 # the generalized eigenvalue problem to a regular one, provided
363 # the matrix B can be safely inverted.
365 lu_b = kla.lu_factor(B)
366 if not stabilization[1]:
367 rcond = kla.rcond_from_lu(lu_b, npl.norm(B, 1))
368 # A more stringent condition is used here since errors can
369 # accumulate from here to the eigenvalue calculation later.
370 stabilization[1] = rcond > eps * tol
372 if stabilization[1]:
373 matrices = (kla.lu_solve(lu_b, A), None)
374 else:
375 matrices = (A, B)
376 return Linsys(matrices, v_out, extract_wf)
379def unified_eigenproblem(a, b=None, tol=1e6):
380 """A helper routine for modes(), that wraps eigenproblems.
382 This routine wraps the regular and general eigenproblems that can arise
383 in a unfied way.
385 Parameters
386 ----------
387 a : numpy array
388 The matrix on the left hand side of a regular or generalized eigenvalue
389 problem.
390 b : numpy array or None
391 The matrix on the right hand side of the generalized eigenvalue problem.
392 tol : float
393 The tolerance for separating eigenvalues with absolute value 1 from the
394 rest.
396 Returns
397 -------
398 ev : numpy array
399 An array of eigenvalues (can contain NaNs and Infs, but those
400 are not accessed in `modes()`) The number of eigenvalues equals
401 twice the number of nonzero singular values of
402 `h_hop` (so `2*h_cell.shape[0]` if `h_hop` is invertible).
403 evanselect : numpy integer array
404 Index array of right-decaying modes.
405 propselect : numpy integer array
406 Index array of propagating modes (both left and right).
407 vec_gen(select) : function
408 A function that computes the eigenvectors chosen by the array select.
409 ord_schur(select) : function
410 A function that computes the unitary matrix (corresponding to the right
411 eigenvector space) of the (general) Schur decomposition reordered such
412 that the eigenvalues chosen by the array select are in the top left
413 block.
414 """
415 if b is None:
416 eps = np.finfo(a.dtype).eps * tol
417 t, z, ev = kla.schur(a)
419 # Right-decaying modes.
420 select = abs(ev) > 1 + eps
421 # Propagating modes.
422 propselect = abs(abs(ev) - 1) < eps
424 vec_gen = lambda x: kla.evecs_from_schur(t, z, select=x)
425 ord_schur = lambda x: kla.order_schur(x, t, z, calc_ev=False)[1]
427 else:
428 eps = np.finfo(np.common_type(a, b)).eps * tol
429 s, t, z, alpha, beta = kla.gen_schur(a, b, calc_q=False)
431 # Right-decaying modes.
432 select = abs(alpha) > (1 + eps) * abs(beta)
433 # Propagating modes.
434 propselect = (abs(abs(alpha) - abs(beta)) < eps * abs(beta))
436 with np.errstate(divide='ignore', invalid='ignore'):
437 ev = alpha / beta
438 # Note: the division is OK here, since we later only access
439 # eigenvalues close to the unit circle
441 vec_gen = lambda x: kla.evecs_from_gen_schur(s, t, z=z, select=x)
442 ord_schur = lambda x: kla.order_gen_schur(x, s, t, z=z,
443 calc_ev=False)[2]
445 return ev, select, propselect, vec_gen, ord_schur
448def phs_symmetrization(wfs, particle_hole):
449 """Makes the wave functions that have the same velocity at a time-reversal
450 invariant momentum (TRIM) particle-hole symmetric.
452 If P is the particle-hole operator and P^2 = 1, then a particle-hole
453 symmetric wave function at a TRIM is an eigenstate of P with eigenvalue 1.
454 If P^2 = -1, wave functions with the same velocity at a TRIM come in pairs.
455 Such a pair is particle-hole symmetric if the wave functions are related by
456 P, i. e. the pair can be expressed as [psi_n, P psi_n] where psi_n is a wave
457 function.
459 To ensure proper ordering of modes, this function also returns an array
460 of indices which ensures that particle-hole partners are properly ordered
461 in this subspace of modes. These are later used with np.lexsort to ensure
462 proper ordering.
464 Parameters
465 ----------
466 wfs : numpy array
467 A matrix of propagating wave functions at a TRIM that all have the same
468 velocity. The orthonormal wave functions form the columns of this matrix.
469 particle_hole : numpy array
470 The matrix representation of the unitary part of the particle-hole
471 operator, expressed in the tight binding basis.
473 Returns
474 -------
475 new_wfs : numpy array
476 The matrix of particle-hole symmetric wave functions.
477 TRIM_sort: numpy integer array
478 Index array that stores the proper sort order of particle-hole symmetric
479 wave functions in this subspace.
480 """
482 def Pdot(mat):
483 """Apply the particle-hole operator to an array. """
484 return particle_hole.dot(mat.conj())
486 # Take P in the subspace of W = wfs: U = W^+ @ P @ W^*.
487 U = wfs.T.conj().dot(Pdot(wfs))
488 # Check that wfs are orthonormal and the space spanned
489 # by them is closed under ph, meaning U is unitary.
490 if not np.allclose(U.dot(U.T.conj()), np.eye(U.shape[0])): 490 ↛ 491line 490 didn't jump to line 491, because the condition on line 490 was never true
491 raise ValueError('wfs are not orthonormal or not closed under particle_hole.')
492 P_squared = U.dot(U.conj())
493 if np.allclose(P_squared, np.eye(U.shape[0])):
494 P_squared = 1
495 elif np.allclose(P_squared, -np.eye(U.shape[0])): 495 ↛ 498line 495 didn't jump to line 498, because the condition on line 495 was never false
496 P_squared = -1
497 else:
498 raise ValueError('particle_hole must square to +1 or -1. P_squared = {}'.format(P_squared))
500 if P_squared == 1:
501 # Use the matrix square root method from
502 # Applied Mathematics and Computation 234 (2014) 380-384.
503 assert np.allclose(U, U.T)
504 # Schur decomposition guarantees that vecs are orthonormal.
505 vals, vecs = la.schur(U)
506 # U should be normal, so vals is diagonal.
507 assert np.allclose(np.diag(np.diag(vals)), vals)
508 vals = np.diag(vals)
509 # Need to take safe square root of U, the branch cut should not go
510 # through any eigenvalues. Otherwise the square root may not be symmetric.
511 # Find largest gap between eigenvalues
512 phases = np.sort(np.angle(vals))
513 dph = np.append(np.diff(phases), phases[0] + 2*np.pi - phases[-1])
514 i = np.argmax(dph)
515 shift = -np.pi - (phases[i] + dph[i]/2)
516 # Take matrix square root with branch cut in largest gap
517 vals = np.sqrt(vals * np.exp(1j * shift)) * np.exp(-0.5j * shift)
518 sqrtU = vecs.dot(np.diag(vals)).dot(vecs.T.conj())
519 # For symmetric U sqrt(U) is also symmetric.
520 assert np.allclose(sqrtU, sqrtU.T)
521 # We want a new basis W_new such that W_new^+ @ P @ W_new^* = 1.
522 # This is achieved by W_new = W @ sqrt(U).
523 new_wfs = wfs.dot(sqrtU)
524 # If P^2 = 1, there is no need to sort the modes further.
525 TRIM_sort = np.zeros((wfs.shape[1],), dtype=int)
526 else:
527 # P^2 = -1.
528 # Iterate over wave functions to construct
529 # particle-hole partners.
530 new_wfs = []
531 # The number of modes. This is always an even number >=2.
532 N_modes = wfs.shape[1]
533 # If there are only two modes in this subspace, they are orthogonal
534 # so we replace the second one with the P applied to the first one.
535 if N_modes == 2:
536 wf = wfs[:, 0]
537 # Store psi_n and P psi_n.
538 new_wfs.append(wf)
539 new_wfs.append(Pdot(wf))
540 # If there are more than two modes, iterate over wave functions
541 # and construct their particle-hole partners one by one.
542 else:
543 # We construct pairs of modes that are particle-hole partners.
544 # Need to iterate over all pairs except the final one.
545 iterations = range((N_modes-2)//2)
546 for i in iterations:
547 # Take a mode psi_n from the basis - the first column
548 # of the matrix of remaining modes.
549 wf = wfs[:, 0]
550 # Store psi_n and P psi_n.
551 new_wfs.append(wf)
552 P_wf = Pdot(wf)
553 new_wfs.append(P_wf)
554 # Remove psi_n and P psi_n from the basis matrix of modes.
555 # First remove psi_n.
556 wfs = wfs[:, 1:]
557 # Now we project the remaining modes onto the orthogonal
558 # complement of P psi_n. projector:
559 projector = wfs.dot(wfs.T.conj()) - \
560 np.outer(P_wf, P_wf.T.conj())
561 # After the projection, the mode matrix is rank deficient -
562 # the span of the column space has dimension one less than
563 # the number of columns.
564 wfs = projector.dot(wfs)
565 wfs = la.qr(wfs, mode='economic', pivoting=True)[0]
566 # Remove the redundant column.
567 wfs = wfs[:, :-1]
568 # If this is the final iteration, we only have two modes
569 # left and can construct particle-hole partners without
570 # the projection.
571 if i == iterations[-1]:
572 assert wfs.shape[1] == 2
573 wf = wfs[:, 0]
574 # Store psi_n and P psi_n.
575 new_wfs.append(wf)
576 new_wfs.append(Pdot(wf))
577 assert np.allclose(wfs.T.conj().dot(wfs),
578 np.eye(wfs.shape[1]))
579 new_wfs = np.hstack([col.reshape(len(col), 1)/npl.norm(col) for
580 col in new_wfs])
581 assert np.allclose(new_wfs[:, 1::2], Pdot(new_wfs[:, ::2]))
582 # Store sort ordering in this subspace of modes
583 TRIM_sort = np.arange(new_wfs.shape[1])
584 assert np.allclose(new_wfs.T.conj().dot(new_wfs), np.eye(new_wfs.shape[1]))
585 return new_wfs, TRIM_sort
588def make_proper_modes(lmbdainv, psi, extract, tol, particle_hole,
589 time_reversal, chiral):
590 """
591 Find, normalize and sort the propagating eigenmodes.
593 Special care is taken of the case of degenerate k-values, where the
594 numerically computed modes are typically a superposition of the real
595 modes. In this case, also the proper (orthogonal) modes are computed.
596 """
597 vel_eps = np.finfo(psi.dtype).eps * tol
599 nmodes = psi.shape[1]
600 n = len(psi) // 2
602 # Array for the velocities.
603 velocities = np.empty(nmodes, dtype=float)
605 # Array of indices to sort modes at a TRIM by PHS.
606 TRIM_PHS_sort = np.zeros(nmodes, dtype=int)
608 # Calculate the full wave function in real space.
609 full_psi = extract(psi, lmbdainv)
611 # Cast the types if any of the symmetry operators is complex
612 for symm in time_reversal, particle_hole, chiral:
613 if symm is None:
614 continue
615 full_psi = full_psi.astype(np.common_type(symm, full_psi))
616 psi = psi.astype(np.common_type(symm, psi))
618 # Find clusters of nearby eigenvalues. Since the eigenvalues occupy the
619 # unit circle, special care has to be taken to not introduce a cut at
620 # lambda = -1.
621 eps = np.finfo(lmbdainv.dtype).eps * tol
622 angles = np.angle(lmbdainv)
623 sort_order = np.resize(np.argsort(angles), (2 * len(angles,)))
624 boundaries = np.argwhere(np.abs(np.diff(lmbdainv[sort_order]))
625 > eps).flatten() + 1
627 # Detect the singular case of all eigenvalues equal.
628 if boundaries.shape == (0,) and len(angles):
629 boundaries = np.array([0, len(angles)])
631 for interval in zip(boundaries[:-1], boundaries[1:]):
632 if interval[1] > boundaries[0] + len(angles):
633 break
635 indx = sort_order[interval[0] : interval[1]]
637 # If there is a degenerate eigenvalue with several different
638 # eigenvectors, the numerical routines return some arbitrary
639 # overlap of the real, physical solutions. In order
640 # to figure out the correct wave function, we need to
641 # have the full, not the projected wave functions
642 # (at least to our current knowledge).
644 # Finding the true modes is done in two steps:
646 # 1. The true transversal modes should be orthogonal to each other, as
647 # they share the same Bloch momentum (note that transversal modes with
648 # different Bloch momenta k1 and k2 need not be orthogonal, the full
649 # modes are orthogonal because of the longitudinal dependence e^{i k1
650 # x} and e^{i k2 x}). The modes with the same k are therefore
651 # orthogonalized. Moreover for the velocity to have a proper value the
652 # modes should also be normalized.
654 q, r = la.qr(full_psi[:, indx], mode='economic')
656 # If the eigenvectors were purely real up to this stage,
657 # they will typically become complex after the rotation.
658 if psi.dtype != np.common_type(psi, r): 658 ↛ 659line 658 didn't jump to line 659, because the condition on line 658 was never true
659 psi = psi.astype(np.common_type(psi, r))
660 if full_psi.dtype != np.common_type(full_psi, q): 660 ↛ 661line 660 didn't jump to line 661, because the condition on line 660 was never true
661 full_psi = full_psi.astype(np.common_type(psi, q))
663 full_psi[:, indx] = q
664 psi[:, indx] = la.solve(r.T, psi[:, indx].T).T
666 # 2. Moving infinitesimally away from the degeneracy
667 # point, the modes should diagonalize the velocity
668 # operator (i.e. when they are non-degenerate any more)
669 # The modes are therefore rotated properly such that they
670 # diagonalize the velocity operator.
671 # Note that step 2. does not give a unique result if there are
672 # two modes with the same velocity, or if the modes stay
673 # degenerate even for a range of Bloch momenta (and hence
674 # must have the same velocity). However, this does not matter,
675 # since we are happy with any superposition in this case.
677 vel_op = -1j * dot(psi[n:, indx].T.conj(), psi[:n, indx])
678 vel_op = vel_op + vel_op.T.conj()
679 vel_vals, rot = la.eigh(vel_op)
681 # If the eigenvectors were purely real up to this stage,
682 # they will typically become complex after the rotation.
684 if psi.dtype != np.common_type(psi, rot):
685 psi = psi.astype(np.common_type(psi, rot))
686 if full_psi.dtype != np.common_type(full_psi, rot):
687 full_psi = full_psi.astype(np.common_type(psi, rot))
689 psi[:, indx] = dot(psi[:, indx], rot)
690 full_psi[:, indx] = dot(full_psi[:, indx], rot)
691 velocities[indx] = vel_vals
693 # With particle-hole symmetry, treat TRIMs individually.
694 # Particle-hole conserves velocity.
695 # If P^2 = 1, we can pick modes at a TRIM as particle-hole eigenstates.
696 # If P^2 = -1, a mode at a TRIM and its particle-hole partner are
697 # orthogonal, and we pick modes such that they are related by
698 # particle-hole symmetry.
700 # At a TRIM, propagating translation eigenvalues are +1 or -1.
701 if (particle_hole is not None and
702 (np.abs(np.abs(lmbdainv[indx].real) - 1) < eps).all()):
703 assert not len(indx) % 2
704 # Set the eigenvalues to the exact TRIM values.
705 if (np.abs(lmbdainv[indx].real - 1) < eps).all():
706 lmbdainv[indx] = 1
707 else:
708 # Momenta are the negative arguments of the translation
709 # eigenvalues, as computed below using np.angle. np.angle of -1
710 # is pi, so this assigns k = -pi to modes with translation
711 # eigenvalue -1.
712 lmbdainv[indx] = -1
714 # Original wave functions
715 orig_wf = full_psi[:, indx]
717 # Modes are already sorted by velocity in ascending order, as
718 # returned by la.eigh above. The first half is thus incident,
719 # and the second half outgoing.
720 # Here we work within a subspace of modes with a fixed velocity.
721 # Mostly, this is done to ensure that modes of different velocities
722 # are not mixed when particle-hole partners are constructed for
723 # P^2 = -1. First, we identify which modes have the same velocity.
724 # In each such subspace of modes, we construct wave functions that
725 # are particle-hole partners.
726 vels = velocities[indx]
727 # Velocities are sorted in ascending order. Find the indices of the
728 # last instance of each unique velocity.
729 inds = [ind+1 for ind, vel in enumerate(vels[:-1])
730 if np.abs(vel-vels[ind+1])>vel_eps]
731 inds = [0] + inds + [len(vels)]
732 inds = zip(inds[:-1], inds[1:])
733 # Now possess an iterator over tuples, where each tuple (i,j)
734 # contains the starting and final indices i and j of a submatrix
735 # of the modes matrix, such that all modes in the submatrix
736 # have the same velocity.
738 # Iterate over all submatrices of modes with the same velocity.
739 new_wf = []
740 TRIM_sorts = []
741 for ind_tuple in inds:
742 # Pick out wave functions that have a given velocity
743 wfs = orig_wf[:, slice(*ind_tuple)]
744 # Make particle-hole symmetric modes
745 new_modes, TRIM_sort = phs_symmetrization(wfs, particle_hole)
746 new_wf.append(new_modes)
747 # Store sorting indices of the TRIM modes with the given
748 # velocity.
749 TRIM_sorts.append(TRIM_sort)
750 # Gather into a matrix of modes
751 new_wf = np.hstack(new_wf)
752 # Store the sort order of all modes at the TRIM.
753 # Used later with np.lexsort when the ordering
754 # of modes is done.
755 TRIM_PHS_sort[indx] = np.hstack(TRIM_sorts)
756 # Replace the old modes.
757 full_psi[:, indx] = new_wf
758 # For both cases P^2 = +-1, must rotate wave functions in the
759 # singular value basis. Find the rotation from new basis to old.
760 rot = new_wf.T.conj().dot(orig_wf)
761 # Rotate the wave functions in the singular value basis
762 psi[:, indx] = psi[:, indx].dot(rot.T.conj())
764 # Ensure proper usage of chiral symmetry.
765 if chiral is not None and time_reversal is None:
766 out_orig = full_psi[:, indx[len(indx)//2:]]
767 out = chiral.dot(full_psi[:, indx[:len(indx)//2]])
768 # No least squares below because the modes should be orthogonal.
769 rot = out_orig.T.conj().dot(out)
770 full_psi[:, indx[len(indx)//2:]] = out
771 psi[:, indx[len(indx)//2:]] = psi[:, indx[len(indx)//2:]].dot(rot)
773 if np.any(abs(velocities) < vel_eps): 773 ↛ 774line 773 didn't jump to line 774, because the condition on line 773 was never true
774 raise RuntimeError("Found a mode with zero or close to zero velocity.")
775 if 2 * np.sum(velocities < 0) != len(velocities): 775 ↛ 776line 775 didn't jump to line 776, because the condition on line 775 was never true
776 raise RuntimeError("Numbers of left- and right-propagating "
777 "modes differ, possibly due to a numerical "
778 "instability.")
780 momenta = -np.angle(lmbdainv)
781 # Sort the modes. The modes are sorted first by velocity and momentum,
782 # and finally TRIM modes are properly ordered.
783 order = np.lexsort([TRIM_PHS_sort, velocities,
784 -np.sign(velocities) * momenta, np.sign(velocities)])
786 velocities = velocities[order]
787 momenta = momenta[order]
788 full_psi = full_psi[:, order]
789 psi = psi[:, order]
791 # Use particle-hole symmetry to relate modes that propagate in the
792 # same direction but at opposite momenta.
793 # Modes are sorted by velocity (first incident, then outgoing).
794 # Incident modes are then sorted by momentum in ascending order,
795 # and outgoing modes in descending order.
796 # Adopted convention is to get modes with negative k (both in and out)
797 # by applying particle-hole operator to modes with positive k.
798 if particle_hole is not None:
799 N = nmodes//2 # Number of incident or outgoing modes.
800 # With particle-hole symmetry, N must be an even number.
801 # Incident modes
802 positive_k = (np.pi - eps > momenta[:N]) * (momenta[:N] > eps)
803 # Original wave functions with negative values of k
804 orig_neg_k = full_psi[:, :N][:, positive_k[::-1]]
805 # For incident modes, ordering of wfs by momentum as returned by kwant
806 # is [-k2, -k1, k1, k2], if k2, k1 > 0 and k2 > k1.
807 # To maintain this ordering with ki and -ki as particle-hole partners,
808 # reverse the order of the product at the end.
809 wf_neg_k = particle_hole.dot(
810 (full_psi[:, :N][:, positive_k]).conj())[:, ::-1]
811 rot = la.lstsq(orig_neg_k, wf_neg_k)[0]
812 full_psi[:, :N][:, positive_k[::-1]] = wf_neg_k
813 psi[:, :N][:, positive_k[::-1]] = \
814 psi[:, :N][:, positive_k[::-1]].dot(rot)
816 # Outgoing modes
817 positive_k = (np.pi - eps > momenta[N:]) * (momenta[N:] > eps)
818 # Original wave functions with negative values of k
819 orig_neg_k = full_psi[:, N:][:, positive_k[::-1]]
820 # For outgoing modes, ordering of wfs by momentum as returned by kwant
821 # is like [k2, k1, -k1, -k2], if k2, k1 > 0 and k2 > k1.
823 # Reverse order at the end to match momenta of opposite sign.
824 wf_neg_k = particle_hole.dot(
825 full_psi[:, N:][:, positive_k].conj())[:, ::-1]
826 rot = la.lstsq(orig_neg_k, wf_neg_k)[0]
827 full_psi[:, N:][:, positive_k[::-1]] = wf_neg_k
828 psi[:, N:][:, positive_k[::-1]] = \
829 psi[:, N:][:, positive_k[::-1]].dot(rot)
831 # Modes are ordered by velocity.
832 # Use time-reversal symmetry to relate modes of opposite velocity.
833 if time_reversal is not None:
834 # Note: within this function, nmodes refers to the total number
835 # of propagating modes, not either left or right movers.
836 out_orig = full_psi[:, nmodes//2:]
837 out = time_reversal.dot(full_psi[:, :nmodes//2].conj())
838 rot = la.lstsq(out_orig, out)[0]
839 full_psi[:, nmodes//2:] = out
840 psi[:, nmodes//2:] = psi[:, nmodes//2:].dot(rot)
842 norm = np.sqrt(abs(velocities))
843 full_psi = full_psi / norm
844 psi = psi / norm
846 return psi, PropagatingModes(full_psi, velocities, momenta)
849def compute_block_modes(h_cell, h_hop, tol, stabilization,
850 time_reversal, particle_hole, chiral):
851 """Calculate modes corresponding to a single projector. """
852 n, m = h_hop.shape
854 # Defer most of the calculation to helper routines.
855 matrices, v, extract = setup_linsys(h_cell, h_hop, tol, stabilization)
856 ev, evanselect, propselect, vec_gen, ord_schur = unified_eigenproblem(
857 *(matrices + (tol,)))
859 # v is never None.
860 # h_hop.shape[0] and v.shape[1] not always the same.
861 # Adding this makes the difference for some tests
862 n = v.shape[1]
864 nrightmovers = np.sum(propselect) // 2
865 nevan = n - nrightmovers
866 evan_vecs = ord_schur(evanselect)[:, :nevan]
868 # Compute the propagating eigenvectors.
869 prop_vecs = vec_gen(propselect)
870 # Compute their velocity, and, if necessary, rotate them.
872 # prop_vecs here is 'psi' in make_proper_modes, i.e. the wf in the SVD
873 # basis. It is in turn used to construct vecs and vecslmbdainv (the
874 # propagating parts). The evanescent parts of vecs and vecslmbdainv come
875 # from evan_vecs above.
876 prop_vecs, real_space_data = make_proper_modes(ev[propselect], prop_vecs,
877 extract, tol, particle_hole,
878 time_reversal, chiral)
880 vecs = np.c_[prop_vecs[n:], evan_vecs[n:]]
881 vecslmbdainv = np.c_[prop_vecs[:n], evan_vecs[:n]]
883 # Prepare output for a single block
884 wave_functions = real_space_data.wave_functions
885 momenta = real_space_data.momenta
886 velocities = real_space_data.velocities
888 return (wave_functions, momenta, velocities, vecs, vecslmbdainv, v)
891def transform_modes(modes_data, unitary=None, time_reversal=None,
892 particle_hole=None, chiral=None):
893 """Transform the modes data for a given block of the Hamiltonian using a
894 discrete symmetry (see arguments). The symmetry operator can also be
895 specified as can also be identity, for the case when blocks are identical.
897 Assume that modes_data has the form returned by compute_block_modes, i.e. a
898 tuple (wave_functions, momenta, velocities, nmodes, vecs, vecslmbdainv, v)
899 containing the block modes data.
901 Assume that the symmetry operator is written in the proper basis (block
902 basis, not full TB).
904 """
905 wave_functions, momenta, velocities, vecs, vecslmbdainv, v = modes_data
907 # Copy to not overwrite modes from previous blocks
908 wave_functions = wave_functions.copy()
909 momenta = momenta.copy()
910 velocities = velocities.copy()
911 vecs = vecs.copy()
912 vecslmbdainv = vecslmbdainv.copy()
913 v = v.copy()
915 nmodes = wave_functions.shape[1] // 2
917 if unitary is not None:
918 perm = np.arange(2*nmodes)
919 conj = False
920 flip_energy = False
921 elif time_reversal is not None:
922 unitary = time_reversal
923 perm = np.arange(2*nmodes)[::-1]
924 conj = True
925 flip_energy = False
926 elif particle_hole is not None:
927 unitary = particle_hole
928 perm = ((-1-np.arange(2*nmodes)) % nmodes +
929 nmodes * (np.arange(2*nmodes) // nmodes))
930 conj = True
931 flip_energy = True
932 elif chiral is not None:
933 unitary = chiral
934 perm = (np.arange(2*nmodes) % nmodes +
935 nmodes * (np.arange(2*nmodes) < nmodes))
936 conj = False
937 flip_energy = True
938 else: # skip coverage
939 raise ValueError("No relation between blocks was provided.")
941 if conj:
942 momenta *= -1
943 vecs = vecs.conj()
944 vecslmbdainv = vecslmbdainv.conj()
945 wave_functions = wave_functions.conj()
946 v = v.conj()
948 if flip_energy:
949 vecslmbdainv *= -1
951 if flip_energy != conj:
952 velocities *= -1
954 wave_functions = unitary.dot(wave_functions)[:, perm]
955 v = unitary.dot(v)
956 vecs[:, :2*nmodes] = vecs[:, perm]
957 vecslmbdainv[:, :2*nmodes] = vecslmbdainv[:, perm]
958 velocities = velocities[perm]
959 momenta = momenta[perm]
960 return (wave_functions, momenta, velocities, vecs, vecslmbdainv, v)
963def modes(h_cell, h_hop, tol=1e6, stabilization=None, *,
964 discrete_symmetry=None, projectors=None, time_reversal=None,
965 particle_hole=None, chiral=None):
966 """Compute the eigendecomposition of a translation operator of a lead.
968 Parameters
969 ----------
970 h_cell : numpy array, real or complex, shape (N,N) The unit cell
971 Hamiltonian of the lead unit cell.
972 h_hop : numpy array, real or complex, shape (N,M)
973 The hopping matrix from a lead cell to the one on which self-energy
974 has to be calculated (and any other hopping in the same direction).
975 tol : float
976 Numbers and differences are considered zero when they are smaller
977 than `tol` times the machine precision.
978 stabilization : sequence of 2 booleans or None
979 Which steps of the eigenvalue prolem stabilization to perform. If the
980 value is `None`, then Kwant chooses the fastest (and least stable)
981 algorithm that is expected to be sufficient. For any other value,
982 Kwant forms the eigenvalue problem in the basis of the hopping singular
983 values. The first element set to `True` forces Kwant to add an
984 anti-Hermitian term to the cell Hamiltonian before inverting. If it is
985 set to `False`, the extra term will only be added if the cell
986 Hamiltonian isn't invertible. The second element set to `True` forces
987 Kwant to solve a generalized eigenvalue problem, and not to reduce it
988 to the regular one. If it is `False`, reduction to a regular problem
989 is performed if possible. Selecting the stabilization manually is
990 mostly necessary for testing purposes.
991 particle_hole : sparse or dense square matrix
992 The unitary part of the particle-hole symmetry operator.
993 time_reversal : sparse or dense square matrix
994 The unitary part of the time-reversal symmetry operator.
995 chiral : sparse or dense square matrix
996 The chiral symmetry operator.
997 projectors : an iterable of sparse or dense matrices
998 Projectors that block diagonalize the Hamiltonian in accordance
999 with a conservation law.
1001 Returns
1002 -------
1003 propagating : `~kwant.physics.PropagatingModes`
1004 Contains the array of the wave functions of propagating modes, their
1005 momenta, and their velocities. It can be used to identify the gauge in
1006 which the scattering problem is solved.
1007 stabilized : `~kwant.physics.StabilizedModes`
1008 A basis of propagating and evanescent modes used by the solvers.
1010 Notes
1011 -----
1012 The sorting of the propagating modes is fully described in the
1013 documentation for `~kwant.physics.PropagatingModes`. In simple cases where
1014 bands do not cross, this ordering corresponds to "lowest modes first". In
1015 general, however, it is necessary to examine the band structure --
1016 something this function is not doing by design.
1018 Propagating modes with the same momentum are orthogonalized. All the
1019 propagating modes are normalized by current.
1021 `projectors`, `time_reversal`, `particle_hole`, and `chiral` affect the
1022 basis in which the scattering modes are expressed - see
1023 `~kwant.physics.DiscreteSymmetry` for details.
1025 This function uses the most stable and efficient algorithm for calculating
1026 the mode decomposition that the Kwant authors are aware about. Its details
1027 are to be published.
1028 """
1029 if discrete_symmetry is not None:
1030 projectors, time_reversal, particle_hole, chiral = discrete_symmetry
1031 n, m = h_hop.shape
1033 if h_cell.shape != (n, n): 1033 ↛ 1034line 1033 didn't jump to line 1034, because the condition on line 1033 was never true
1034 raise ValueError("Incompatible matrix sizes for h_cell and h_hop.")
1036 if not np.any(h_hop):
1037 wf = np.zeros((n, 0))
1038 v = np.zeros((m, 0))
1039 m = np.zeros((0, 0))
1040 vec = np.zeros((0,))
1041 return (PropagatingModes(wf, vec, vec), StabilizedModes(m, m, 0, v))
1043 ham = h_cell
1044 # Avoid the trouble of dealing with non-square hopping matrices.
1045 # TODO: How to avoid this while not doing a lot of book-keeping?
1046 hop = np.empty_like(ham, dtype=h_hop.dtype)
1047 hop[:, :m] = h_hop
1048 hop[:, m:] = 0
1050 # Provide default values to not deal with special cases.
1051 if projectors is None:
1052 projectors = [sp_identity(ham.shape[0], format='csr')]
1054 def to_dense_or_csr(matrix):
1055 if matrix is None:
1056 return csr_matrix(ham.shape)
1057 try:
1058 return matrix.tocsr()
1059 except AttributeError:
1060 return matrix
1062 symmetries = map(to_dense_or_csr,
1063 (time_reversal, particle_hole, chiral))
1064 time_reversal, particle_hole, chiral = symmetries
1066 offsets = np.cumsum([0] + [projector.shape[1] for projector in projectors])
1067 indices = [slice(*i) for i in np.vstack([offsets[:-1], offsets[1:]]).T]
1068 projection_op = sp_hstack(projectors)
1070 def basis_change(a, antiunitary=False):
1071 b = projection_op
1072 # We need the extra transposes to ensure that sparse dot is used.
1073 if antiunitary:
1074 # b.T.conj() @ a @ b.conj()
1075 return (b.T.conj().dot((b.T.conj().dot(a)).T)).T
1076 else:
1077 # b.T.conj() @ a @ b
1078 return (b.T.dot((b.T.conj().dot(a)).T)).T
1080 # Conservation law basis
1081 ham_cons = basis_change(ham)
1082 hop_cons = basis_change(hop)
1083 trs = basis_change(time_reversal, True)
1084 phs = basis_change(particle_hole, True)
1085 sls = basis_change(chiral)
1087 # Check that the Hamiltonian has the conservation law
1088 block_modes = len(projectors) * [None]
1089 numbers_coords = combinations_with_replacement(enumerate(indices), 2)
1090 for (i, x), (j, y) in numbers_coords:
1091 h = ham_cons[x, y]
1092 t = hop_cons[x, y]
1093 # Symmetries that project from block x to block y
1094 symmetries = [symm[y, x] for symm in (trs, phs, sls)]
1095 symmetries = [(symm if nonzero_symm_projection(symm) else None) for
1096 symm in symmetries]
1097 if i == j:
1098 if block_modes[i] is not None:
1099 continue
1100 # We did not compute this block yet.
1101 block_modes[i] = compute_block_modes(h, t, tol, stabilization,
1102 *symmetries)
1103 else:
1104 if block_modes[j] is not None: 1104 ↛ 1106line 1104 didn't jump to line 1106, because the condition on line 1104 was never true
1105 # Modes in the block already computed.
1106 continue
1107 x, y = tuple(np.mgrid[x, x]), tuple(np.mgrid[y, y])
1109 if ham_cons[x].shape != ham_cons[y].shape:
1110 continue
1111 if (np.allclose(ham_cons[x], ham_cons[y]) and
1112 np.allclose(hop_cons[x], hop_cons[y])):
1113 unitary = sp_identity(h.shape[0])
1114 else:
1115 unitary = None
1116 if any(op is not None for op in symmetries + [unitary]):
1117 block_modes[j] = transform_modes(block_modes[i], unitary,
1118 *symmetries)
1119 (wave_functions, momenta, velocities,
1120 vecs, vecslmbdainv, sqrt_hops) = zip(*block_modes)
1122 # Reorder by direction of propagation
1123 wave_functions = group_halves([(projector.dot(wf)).T for wf, projector in
1124 zip(wave_functions, projectors)]).T
1125 # Propagating modes object to return
1126 prop_modes = PropagatingModes(wave_functions, group_halves(velocities),
1127 group_halves(momenta))
1128 nmodes = [len(v) // 2 for v in velocities]
1129 # Store the number of modes per block as an attribute.
1130 # nmodes[i] is the number of left or right moving modes in block i.
1131 # In the module that makes leads with conservation laws, this is necessary
1132 # to keep track of the block structure of the scattering matrix.
1133 prop_modes.block_nmodes = nmodes
1135 parts = zip(*([v[:, :n], v[:, n:2*n], v[:, 2*n:]]
1136 for n, v in zip(nmodes, vecs)))
1137 vecs = np.hstack([block_diag(*part) for part in parts])
1139 parts = zip(*([v[:, :n], v[:, n:2*n], v[:, 2*n:]]
1140 for n, v in zip(nmodes, vecslmbdainv)))
1141 vecslmbdainv = np.hstack([block_diag(*part) for part in parts])
1143 sqrt_hops = np.hstack([projector.dot(hop) for projector, hop in
1144 zip(projectors, sqrt_hops)])[:h_hop.shape[1]]
1146 stab_modes = StabilizedModes(vecs, vecslmbdainv, sum(nmodes), sqrt_hops)
1148 return prop_modes, stab_modes
1151def selfenergy(h_cell, h_hop, tol=1e6):
1152 """
1153 Compute the self-energy generated by the lead.
1155 Parameters
1156 ----------
1157 h_cell : numpy array, real or complex, shape (N,N) The unit cell Hamiltonian
1158 of the lead unit cell.
1159 h_hop : numpy array, real or complex, shape (N,M)
1160 The hopping matrix from a lead cell to the one on which self-energy
1161 has to be calculated (and any other hopping in the same direction).
1162 tol : float
1163 Numbers are considered zero when they are smaller than `tol` times
1164 the machine precision.
1166 Returns
1167 -------
1168 Sigma : numpy array, real or complex, shape (M,M)
1169 The computed self-energy. Note that even if `h_cell` and `h_hop` are
1170 both real, `Sigma` will typically be complex. (More precisely, if there
1171 is a propagating mode, `Sigma` will definitely be complex.)
1173 Notes
1174 -----
1175 For simplicity this function internally calculates the modes first.
1176 This may cause a small slowdown, and can be improved if necessary.
1177 """
1178 stabilized = modes(h_cell, h_hop, tol)[1]
1179 return stabilized.selfenergy()
1182def square_selfenergy(width, hopping, fermi_energy):
1183 """
1184 Calculate analytically the self energy for a square lattice.
1186 The lattice is assumed to have a single orbital per site and
1187 nearest-neighbor hopping.
1189 Parameters
1190 ----------
1191 width : integer
1192 width of the lattice
1193 """
1195 # Following appendix C of M. Wimmer's diploma thesis:
1196 # http://www.physik.uni-regensburg.de/forschung/\
1197 # richter/richter/media/research/publications2004/wimmer-Diplomarbeit.pdf
1199 # p labels transversal modes. i and j label the sites of a cell.
1201 # Precalculate the transverse wave function.
1202 psi_p_i = np.empty((width, width))
1203 factor = pi / (width + 1)
1204 prefactor = sqrt(2 / (width + 1))
1205 for p in range(width):
1206 for i in range(width):
1207 psi_p_i[p, i] = prefactor * sin(factor * (p + 1) * (i + 1))
1209 # Precalculate the integrals of the longitudinal wave functions.
1210 def f(q):
1211 if abs(q) <= 2:
1212 return q/2 - 1j * sqrt(1 - (q / 2) ** 2)
1213 else:
1214 return q/2 - copysign(sqrt((q / 2) ** 2 - 1), q)
1215 f_p = np.empty((width,), dtype=complex)
1216 for p in range(width):
1217 e = 2 * hopping * (1 - cos(factor * (p + 1)))
1218 q = (fermi_energy - e) / hopping - 2
1219 f_p[p] = f(q)
1221 # Put everything together into the self energy and return it.
1222 result = np.empty((width, width), dtype=complex)
1223 for i in range(width):
1224 for j in range(width):
1225 result[i, j] = hopping * (
1226 psi_p_i[:, i] * psi_p_i[:, j].conj() * f_p).sum()
1227 return result