leads.py 27.2 KB
Newer Older
1
# Copyright 2011-2013 Kwant authors.
2
#
3
# This file is part of Kwant.  It is subject to the license terms in the
4
# LICENSE file found in the top-level directory of this distribution and at
5
# http://kwant-project.org/license.  A list of Kwant authors can be found in
6
7
8
# the AUTHORS file at the top-level directory of this distribution and at
# http://kwant-project.org/authors.

9
10
from __future__ import division
from math import sin, cos, sqrt, pi, copysign
11
from collections import namedtuple
12
from itertools import izip
13
14
15
import numpy as np
import numpy.linalg as npl
import scipy.linalg as la
16
from .. import linalg as kla
17
18
19

dot = np.dot

20
__all__ = ['selfenergy', 'modes', 'PropagatingModes', 'StabilizedModes']
21

Anton Akhmerov's avatar
Anton Akhmerov committed
22

23
24
25
26
27
# Container classes
Linsys = namedtuple('Linsys', ['eigenproblem', 'v', 'extract'])


class PropagatingModes(object):
28
    """The calculated propagating modes of a lead.
29

30
31
    Attributes
    ----------
32
    wave_functions : numpy array
33
        The wave functions of the propagating modes.
34
35
36
37
38
39
40
41
42
43
44
    momenta : numpy array
        Momenta of the modes.
    velocities : numpy array
        Velocities of the modes.

    Notes
    =====
    The sort order of all the three arrays is identical. The first half of the
    modes have negative velocity, the second half have positive velocity. The
    modes with negative velocity are ordered from larger to lower momenta, the
    modes with positive velocity vice versa.
45

Christoph Groth's avatar
Christoph Groth committed
46
47
48
49
50
    The first dimension of `wave_functions` corresponds to the orbitals of all
    the sites in a unit cell, the second one to the number of the mode.  Each
    mode is normalized to carry unit current. If several modes have the same
    momentum and velocity, an arbitrary orthonormal basis in the subspace of
    these modes is chosen.
51
52
53
54
55
56
57
58
59
60
61
    """
    def __init__(self, wave_functions, velocities, momenta):
        kwargs = locals()
        kwargs.pop('self')
        self.__dict__.update(kwargs)


class StabilizedModes(object):
    """Stabilized eigendecomposition of the translation operator.

    Due to the lack of Hermiticity of the translation operator, its
62
    eigendecomposition is frequently poorly conditioned. Solvers in Kwant use
63
64
65
66
67
68
69
70
71
72
73
74
    this stabilized decomposition of the propagating and evanescent modes in
    the leads. If the hopping between the unit cells of an infinite system is
    invertible, the translation eigenproblem is written in the basis `psi_n,
    h_hop^+ * psi_(n+1)`, with ``h_hop`` the hopping between unit cells.  If
    `h_hop` is not invertible, and has the singular value decomposition `u s
    v^+`, then the eigenproblem is written in the basis `sqrt(s) v^+ psi_n,
    sqrt(s) u^+ psi_(n+1)`. In this basis we calculate the eigenvectors of the
    propagating modes, and the Schur vectors (an orthogonal basis) in the space
    of evanescent modes.

    `vecs` and `vecslmbdainv` are the first and the second halves of the wave
    functions.  The first `nmodes` are eigenmodes moving in the negative
75
    direction (hence they are incoming into the system in Kwant convention),
76
77
78
79
80
81
    the second `nmodes` are eigenmodes moving in the positive direction. The
    remaining modes are the Schur vectors of the modes evanescent in the
    positive direction. Propagating modes with the same eigenvalue are
    orthogonalized, and all the propagating modes are normalized to carry unit
    current. Finally the `sqrt_hop` attribute is `v sqrt(s)`.

82
83
    Attributes
    ----------
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
    vecs : numpy array
        Translation eigenvectors.
    vecslmbdainv : numpy array
        Translation eigenvectors divided by the corresponding eigenvalues.
    nmodes : int
        Number of left-moving (or right-moving) modes.
    sqrt_hop : numpy array or None
        Part of the SVD of `h_hop`, or None if the latter is invertible.
    """

    def __init__(self, vecs, vecslmbdainv, nmodes, sqrt_hop=None):
        kwargs = locals()
        kwargs.pop('self')
        self.__dict__.update(kwargs)

99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
    def selfenergy(self):
        """
        Compute the self-energy generated by lead modes.

        Returns
        -------
        Sigma : numpy array, real or complex, shape (M,M)
            The computed self-energy. Note that even if `h_cell` and `h_hop` are
            both real, `Sigma` will typically be complex. (More precisely, if
            there is a propagating mode, `Sigma` will definitely be complex.)
        """
        v = self.sqrt_hop
        vecs = self.vecs[:, self.nmodes:]
        vecslmbdainv = self.vecslmbdainv[:, self.nmodes:]
        if v is not None:
            return dot(v, dot(vecs, la.solve(vecslmbdainv, v.T.conj())))
        else:
            return la.solve(vecslmbdainv.T, vecs.T).T

118
119

# Auxiliary functions that perform different parts of the calculation.
120
121
def setup_linsys(h_cell, h_hop, tol=1e6, stabilization=None):
    """Make an eigenvalue problem for eigenvectors of translation operator.
122
123
124

    Parameters
    ----------
125
126
    h_cell : numpy array with shape (n, n)
        Hamiltonian of a single lead unit cell.
127
    h_hop : numpy array with shape (n, m), m <= n
128
        Hopping Hamiltonian from a cell to the next one.
129
130
131
    tol : float
        Numbers are considered zero when they are smaller than `tol` times
        the machine precision.
132
    stabilization : list of 2 booleans or None
Christoph Groth's avatar
Christoph Groth committed
133
        Which steps of the eigenvalue problem stabilization to perform. If the
134
        value is `None`, then Kwant chooses the fastest (and least stable)
135
        algorithm that is expected to be sufficient.  For any other value,
136
137
        Kwant forms the eigenvalue problem in the basis of the hopping singular
        values.  The first element set to `True` forces Kwant to add an
138
139
140
        anti-Hermitian term to the cell Hamiltonian before inverting. If it is
        set to `False`, the extra term will only be added if the cell
        Hamiltonian isn't invertible. The second element set to `True` forces
141
        Kwant to solve a generalized eigenvalue problem, and not to reduce it
142
        to the regular one.  If it is `False`, reduction to a regular problem
143
        is performed if possible.
144
145
146

    Returns
    -------
147
    linsys : namedtuple
148
149
150
151
152
        A named tuple containing `matrices` a matrix pencil defining
        the eigenproblem, `v` a hermitian conjugate of the last matrix in
        the hopping singular value decomposition, and functions for
        extracting the wave function in the unit cell from the wave function
        in the basis of the nonzero singular exponents of the hopping.
153
154
155

    Notes
    -----
156
157
    The lead problem with degenerate hopping is rather complicated, and the
    details of the algorithm will be published elsewhere.
158

159
    """
160
    n = h_cell.shape[0]
161
    m = h_hop.shape[1]
162
163
    if stabilization is not None:
        stabilization = list(stabilization)
164

Anton Akhmerov's avatar
Anton Akhmerov committed
165
    if not (np.any(h_hop.real) or np.any(h_hop.imag)):
166
        # Inter-cell hopping is zero.  The current algorithm is not suited to
167
168
169
        # treat this extremely singular case.
        # Note: np.any(h_hop) returns (at least from numpy 1.6.*)
        #       False if h_hop is purely imaginary
170
        raise ValueError("Inter-cell hopping is exactly zero.")
171

172
    # If both h and t are real, it may be possible to use the real eigenproblem.
173
    if (not np.any(h_hop.imag)) and (not np.any(h_cell.imag)):
174
        h_hop = h_hop.real
175
        h_cell = h_cell.real
176

177
    eps = np.finfo(np.common_type(h_cell, h_hop)).eps * tol
178

179
180
181
    # First check if the hopping matrix has singular values close to 0.
    # (Close to zero is defined here as |x| < eps * tol * s[0] , where
    # s[0] is the largest singular value.)
182

183
    u, s, vh = la.svd(h_hop)
184
    assert m == vh.shape[1], "Corrupt output of svd."
185
    n_nonsing = np.sum(s > eps * s[0])
186

187
    if (n_nonsing == n and stabilization is None):
188
        # The hopping matrix is well-conditioned and can be safely inverted.
189
        # Hence the regular transfer matrix may be used.
190
        hop_inv = la.inv(h_hop)
191

192
193
        A = np.zeros((2*n, 2*n), dtype=np.common_type(h_cell, h_hop))
        A[:n, :n] = dot(hop_inv, -h_cell)
194
195
        A[:n, n:] = -hop_inv
        A[n:, :n] = h_hop.T.conj()
196

197
        # The function that can extract the full wave function psi from the
198
199
200
201
        # projected one. Here it is almost trivial, but used for simplifying
        # the logic.

        def extract_wf(psi, lmbdainv):
202
            return np.copy(psi[:n])
203

204
        matrices = (A, None)
205
        v_out = None
206
    else:
207
208
        if stabilization is None:
            stabilization = [None, False]
209

210
211
212
213
214
        # The hopping matrix has eigenvalues close to 0 - those
        # need to be eliminated.

        # Recast the svd of h_hop = u s v^dagger such that
        # u, v are matrices with shape n x n_nonsing.
Anton Akhmerov's avatar
Anton Akhmerov committed
215
216
        u = u[:, :n_nonsing]
        s = s[:n_nonsing]
217
        u = u * np.sqrt(s)
218
219
        # pad v with zeros if necessary
        v = np.zeros((n, n_nonsing), dtype=vh.dtype)
220
        v[:vh.shape[1]] = vh[:n_nonsing].T.conj()
221
        v = v * np.sqrt(s)
222

223
224
225
226
227
228
229
230
231
232
233
        # Eliminating the zero eigenvalues requires inverting the on-site
        # Hamiltonian, possibly including a self-energy-like term.  The
        # self-energy-like term stabilizes the inversion, but the most stable
        # choice is inherently complex. This can be disadvantageous if the
        # Hamiltonian is real, since staying in real arithmetics can be
        # significantly faster.  The strategy here is to add a complex
        # self-energy-like term always if the original Hamiltonian is complex,
        # and check for invertibility first if it is real

        matrices_real = issubclass(np.common_type(h_cell, h_hop), np.floating)
        add_imaginary = stabilization[0] or ((stabilization[0] is None) and
234
                                             not matrices_real)
235
236
237
238
        # Check if there is a chance we will not need to add an imaginary term.
        if not add_imaginary:
            h = h_cell
            sol = kla.lu_factor(h)
239
240
            rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))

241
242
243
            if rcond < eps:
                need_to_stabilize = True
            else:
244
                need_to_stabilize = False
245

246
        if add_imaginary or need_to_stabilize:
247
            need_to_stabilize = True
248
            # Matrices are complex or need self-energy-like term to be
249
            # stabilized.
250
            temp = dot(u, u.T.conj()) + dot(v, v.T.conj())
251
            h = h_cell + 1j * temp
252
253
254
255
256
257

            sol = kla.lu_factor(h)
            rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))

            # If the condition number of the stabilized h is
            # still bad, there is nothing we can do.
258
            if rcond < eps:
259
260
261
                raise RuntimeError("Flat band encountered at the requested "
                                   "energy, result is badly defined.")

262
263
        # The function that can extract the full wave function psi from
        # the projected one (v^dagger psi lambda^-1, s u^dagger psi).
264
265

        def extract_wf(psi, lmbdainv):
266
            wf = - dot(u, psi[: n_nonsing] * lmbdainv) - \
267
268
269
270
271
                 dot(v, psi[n_nonsing:])
            if need_to_stabilize:
                wf += 1j * (dot(v, psi[: n_nonsing]) +
                            dot(u, psi[n_nonsing:] * lmbdainv))
            return kla.lu_solve(sol, wf)
272
273
274

        # Setup the generalized eigenvalue problem.

275
276
        A = np.zeros((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
        B = np.zeros((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
277

278
279
        begin, end = slice(n_nonsing), slice(n_nonsing, None)

Anton Akhmerov's avatar
Anton Akhmerov committed
280
        A[end, begin] = np.identity(n_nonsing)
281
        temp = kla.lu_solve(sol, v)
282
        temp2 = dot(u.T.conj(), temp)
283
        if need_to_stabilize:
Anton Akhmerov's avatar
Anton Akhmerov committed
284
285
            A[begin, begin] = -1j * temp2
        A[begin, end] = temp2
286
        temp2 = dot(v.T.conj(), temp)
287
        if need_to_stabilize:
Anton Akhmerov's avatar
Anton Akhmerov committed
288
289
            A[end, begin] -= 1j *temp2
        A[end, end] = temp2
290

Anton Akhmerov's avatar
Anton Akhmerov committed
291
        B[begin, end] = -np.identity(n_nonsing)
292
293
        temp = kla.lu_solve(sol, u)
        temp2 = dot(u.T.conj(), temp)
Anton Akhmerov's avatar
Anton Akhmerov committed
294
        B[begin, begin] = -temp2
295
        if need_to_stabilize:
Anton Akhmerov's avatar
Anton Akhmerov committed
296
            B[begin, end] += 1j * temp2
297
        temp2 = dot(v.T.conj(), temp)
Anton Akhmerov's avatar
Anton Akhmerov committed
298
        B[end, begin] = -temp2
299
        if need_to_stabilize:
Anton Akhmerov's avatar
Anton Akhmerov committed
300
            B[end, end] = 1j * temp2
301
302

        v_out = v[:m]
303
304
305
306
307
308
309
310
311
312

        # Solving a generalized eigenproblem is about twice as expensive
        # as solving a regular eigenvalue problem.
        # Computing the LU factorization is negligible compared to both
        # (approximately 1/30th of a regular eigenvalue problem).
        # Because of this, it makes sense to try to reduce
        # the generalized eigenvalue problem to a regular one, provided
        # the matrix B can be safely inverted.

        lu_b = kla.lu_factor(B)
313
        if not stabilization[1]:
314
315
316
            rcond = kla.rcond_from_lu(lu_b, npl.norm(B, 1))
            # A more stringent condition is used here since errors can
            # accumulate from here to the eigenvalue calculation later.
317
            stabilization[1] = rcond > eps * tol
318

319
        if stabilization[1]:
320
            matrices = (kla.lu_solve(lu_b, A), None)
321
        else:
322
            matrices = (A, B)
323
324
325
326
327
328
329
330
    return Linsys(matrices, v_out, extract_wf)


def unified_eigenproblem(a, b=None, tol=1e6):
    """A helper routine for modes(), that wraps eigenproblems.

    This routine wraps the regular and general eigenproblems that can arise
    in a unfied way.
331

332
333
334
335
336
337
338
339
340
341
    Parameters
    ----------
    a : numpy array
        The matrix on the left hand side of a regular or generalized eigenvalue
        problem.
    b : numpy array or None
        The matrix on the right hand side of the generalized eigenvalue problem.
    tol : float
        The tolerance for separating eigenvalues with absolute value 1 from the
        rest.
342

343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
    Returns
    -------
    ev : numpy array
        An array of eigenvalues (can contain NaNs and Infs, but those
        are not accessed in `modes()`) The number of eigenvalues equals
        twice the number of nonzero singular values of
        `h_hop` (so `2*h_cell.shape[0]` if `h_hop` is invertible).
    evanselect : numpy integer array
        Index array of right-decaying modes.
    propselect : numpy integer array
        Index array of propagating modes (both left and right).
    vec_gen(select) : function
        A function that computes the eigenvectors chosen by the array select.
    ord_schur(select) : function
        A function that computes the unitary matrix (corresponding to the right
        eigenvector space) of the (general) Schur decomposition reordered such
        that the eigenvalues chosen by the array select are in the top left
        block.
361
    """
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
    if b is None:
        eps = np.finfo(a.dtype).eps * tol
        t, z, ev = kla.schur(a)

        # Right-decaying modes.
        select = abs(ev) > 1 + eps
        # Propagating modes.
        propselect = abs(abs(ev) - 1) < eps

        vec_gen = lambda x: kla.evecs_from_schur(t, z, select=x)
        ord_schur = lambda x: kla.order_schur(x, t, z, calc_ev=False)[1]

    else:
        eps = np.finfo(np.common_type(a, b)).eps * tol
        s, t, z, alpha, beta = kla.gen_schur(a, b, calc_q=False)

        # Right-decaying modes.
        select = abs(alpha) > (1 + eps) * abs(beta)
        # Propagating modes.
        propselect = (abs(abs(alpha) - abs(beta)) < eps * abs(beta))

        warning_settings = np.seterr(divide='ignore', invalid='ignore')
        ev = alpha / beta
        np.seterr(**warning_settings)
        # Note: the division is OK here, since we later only access
        #       eigenvalues close to the unit circle

        vec_gen = lambda x: kla.evecs_from_gen_schur(s, t, z=z, select=x)
        ord_schur = lambda x: kla.order_gen_schur(x, s, t, z=z,
                                                  calc_ev=False)[2]

    return ev, select, propselect, vec_gen, ord_schur


def make_proper_modes(lmbdainv, psi, extract, tol=1e6):
    """
    Find, normalize and sort the propagating eigenmodes.
399
400
401
402
403

    Special care is taken of the case of degenerate k-values, where the
    numerically computed modes are typically a superposition of the real
    modes. In this case, also the proper (orthogonal) modes are computed.
    """
404
    vel_eps = np.finfo(psi.dtype).eps * tol
405
406

    nmodes = psi.shape[1]
407
    n = len(psi) // 2
408
409

    # Array for the velocities.
410
    velocities = np.empty(nmodes, dtype=float)
411

412
413
    # Calculate the full wave function in real space.
    full_psi = extract(psi, lmbdainv)
414

415
416
417
    # Find clusters of nearby eigenvalues. Since the eigenvalues occupy the
    # unit circle, special care has to be taken to not introduce a cut at
    # lambda = -1.
418
419
420
421
422
    eps = np.finfo(lmbdainv.dtype).eps * tol
    angles = np.angle(lmbdainv)
    sort_order = np.resize(np.argsort(angles), (2 * len(angles,)))
    boundaries = np.argwhere(np.abs(np.diff(lmbdainv[sort_order]))
                             > eps).flatten() + 1
423

424
425
426
427
    # Detect the singular case of all eigenvalues equal.
    if boundaries.shape == (0,) and len(angles):
        boundaries = np.array([0, len(angles)])

428
429
430
    for interval in izip(boundaries[:-1], boundaries[1:]):
        if interval[1] > boundaries[0] + len(angles):
            break
431

432
        indx = sort_order[interval[0] : interval[1]]
433

434
435
436
437
438
439
        # If there is a degenerate eigenvalue with several different
        # eigenvectors, the numerical routines return some arbitrary
        # overlap of the real, physical solutions. In order
        # to figure out the correct wave function, we need to
        # have the full, not the projected wave functions
        # (at least to our current knowledge).
440

441
        # Finding the true modes is done in two steps:
442

443
444
445
446
447
448
449
450
        # 1. The true transversal modes should be orthogonal to each other, as
        # they share the same Bloch momentum (note that transversal modes with
        # different Bloch momenta k1 and k2 need not be orthogonal, the full
        # modes are orthogonal because of the longitudinal dependence e^{i k1
        # x} and e^{i k2 x}).  The modes with the same k are therefore
        # orthogonalized. Moreover for the velocity to have a proper value the
        # modes should also be normalized.

451
        q, r = la.qr(full_psi[:, indx], mode='economic')
452

453
454
455
456
457
458
        # If the eigenvectors were purely real up to this stage,
        # they will typically become complex after the rotation.
        if psi.dtype != np.common_type(psi, r):
            psi = psi.astype(np.common_type(psi, r))
        if full_psi.dtype != np.common_type(full_psi, q):
            full_psi = full_psi.astype(np.common_type(psi, q))
459

460
461
        full_psi[:, indx] = q
        psi[:, indx] = la.solve(r.T, psi[:, indx].T).T
462

463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
        # 2. Moving infinitesimally away from the degeneracy
        # point, the modes should diagonalize the velocity
        # operator (i.e. when they are non-degenerate any more)
        # The modes are therefore rotated properly such that they
        # diagonalize the velocity operator.
        # Note that step 2. does not give a unique result if there are
        # two modes with the same velocity, or if the modes stay
        # degenerate even for a range of Bloch momenta (and hence
        # must have the same velocity). However, this does not matter,
        # since we are happy with any superposition in this case.

        vel_op = -1j * dot(psi[n:, indx].T.conj(), psi[:n, indx])
        vel_op = vel_op + vel_op.T.conj()
        vel_vals, rot = la.eigh(vel_op)

        # If the eigenvectors were purely real up to this stage,
        # they will typically become complex after the rotation.

        if psi.dtype != np.common_type(psi, rot):
            psi = psi.astype(np.common_type(psi, rot))
483
484
485
        if full_psi.dtype != np.common_type(full_psi, rot):
            full_psi = full_psi.astype(np.common_type(psi, rot))

486
        psi[:, indx] = dot(psi[:, indx], rot)
487
        full_psi[:, indx] = dot(full_psi[:, indx], rot)
488
489
490
491
        velocities[indx] = vel_vals

    if np.any(abs(velocities) < vel_eps):
        raise RuntimeError("Found a mode with zero or close to zero velocity.")
492
    if 2 * np.sum(velocities < 0) != len(velocities):
493
        raise RuntimeError("Numbers of left- and right-propagating "
494
495
                           "modes differ, possibly due to a numerical "
                           "instability.")
496
497
498
499
    momenta = -np.angle(lmbdainv)
    order = np.lexsort([velocities, -np.sign(velocities) * momenta,
                        np.sign(velocities)])

500
501
502
503
504
505
506
507
508
509
510
    # The following is necessary due to wrong numpy handling of zero length
    # arrays, which is going to be fixed in numpy 1.8.
    if not len(order):
        order = slice(None)
    velocities = velocities[order]
    norm = np.sqrt(abs(velocities))
    full_psi = full_psi[:, order] / norm
    psi = psi[:, order] / norm
    momenta = momenta[order]

    return psi, PropagatingModes(full_psi, velocities, momenta)
511

Anton Akhmerov's avatar
Anton Akhmerov committed
512

513
def modes(h_cell, h_hop, tol=1e6, stabilization=None):
514
    """Compute the eigendecomposition of a translation operator of a lead.
515
516
517

    Parameters
    ----------
518
519
    h_cell : numpy array, real or complex, shape (N,N) The unit cell
        Hamiltonian of the lead unit cell.
520
    h_hop : numpy array, real or complex, shape (N,M)
521
        The hopping matrix from a lead cell to the one on which self-energy
522
        has to be calculated (and any other hopping in the same direction).
523
524
525
    tol : float
        Numbers and differences are considered zero when they are smaller
        than `tol` times the machine precision.
526
527
    stabilization : list of 2 booleans or None
        Which steps of the eigenvalue prolem stabilization to perform. If the
528
        value is `None`, then Kwant chooses the fastest (and least stable)
529
        algorithm that is expected to be sufficient.  For any other value,
530
531
        Kwant forms the eigenvalue problem in the basis of the hopping singular
        values.  The first element set to `True` forces Kwant to add an
532
533
534
        anti-Hermitian term to the cell Hamiltonian before inverting. If it is
        set to `False`, the extra term will only be added if the cell
        Hamiltonian isn't invertible. The second element set to `True` forces
535
        Kwant to solve a generalized eigenvalue problem, and not to reduce it
536
        to the regular one.  If it is `False`, reduction to a regular problem
537
        is performed if possible.  Selecting the stabilization manually is
538
        mostly necessary for testing purposes.
539
540
541

    Returns
    -------
542
    propagating : `~kwant.physics.PropagatingModes`
543
544
        Contains the array of the wave functions of propagating modes, their
        momenta, and their velocities. It can be used to identify the gauge in
545
546
        which the scattering problem is solved.
    stabilized : `~kwant.physics.StabilizedModes`
547
        A basis of propagating and evanescent modes used by the solvers.
548
549
550

    Notes
    -----
551
552
553
554
555
556
    The propagating modes are sorted according to the longitudinal component of
    their k-vector, with incoming modes having k sorted in descending order,
    and outgoing modes having k sorted in ascending order.  In simple cases
    where bands do not cross, this ordering corresponds to "lowest modes
    first". In general, however, it is necessary to examine the band structure
    -- something this function is not doing by design.
557

558
    Propagating modes with the same momentum are orthogonalized. All the
559
560
561
    propagating modes are normalized by current.

    This function uses the most stable and efficient algorithm for calculating
562
    the mode decomposition that the Kwant authors are aware about. Its details
563
    are to be published.
564
565
    """
    m = h_hop.shape[1]
566
    n = h_cell.shape[0]
567

568
569
570
    if (h_cell.shape[0] != h_cell.shape[1] or
        h_cell.shape[0] != h_hop.shape[0]):
        raise ValueError("Incompatible matrix sizes for h_cell and h_hop.")
571

572
    # Note: np.any(h_hop) returns (at least from numpy 1.6.1 - 1.8-devel)
573
574
    #       False if h_hop is purely imaginary
    if not (np.any(h_hop.real) or np.any(h_hop.imag)):
575
576
577
578
        v = np.zeros((0, m))
        return (PropagatingModes(np.zeros((0, n)), np.zeros((0,)),
                                 np.zeros((0,))),
                StabilizedModes(np.zeros((0, 0)), np.zeros((0, 0)), 0, v))
579

580
    # Defer most of the calculation to helper routines.
581
    matrices, v, extract = setup_linsys(h_cell, h_hop, tol, stabilization)
582
583
    ev, evanselect, propselect, vec_gen, ord_schur =\
         unified_eigenproblem(*(matrices + (tol,)))
584

585
586
    if v is not None:
        n = v.shape[1]
587

588
589
    nrightmovers = np.sum(propselect) // 2
    nevan = n - nrightmovers
Anton Akhmerov's avatar
Anton Akhmerov committed
590
    evan_vecs = ord_schur(evanselect)[:, :nevan]
591

592
593
594
    # Compute the propagating eigenvectors.
    prop_vecs = vec_gen(propselect)
    # Compute their velocity, and, if necessary, rotate them
595
596
    prop_vecs, real_space_data = \
            make_proper_modes(ev[propselect], prop_vecs, extract, tol)
597

598
599
    vecs = np.c_[prop_vecs[n:], evan_vecs[n:]]
    vecslmbdainv = np.c_[prop_vecs[:n], evan_vecs[:n]]
600

601
    return real_space_data, StabilizedModes(vecs, vecslmbdainv, nrightmovers, v)
602
603


604
def selfenergy(h_cell, h_hop, tol=1e6):
605
606
607
608
609
    """
    Compute the self-energy generated by the lead.

    Parameters
    ----------
610
611
    h_cell : numpy array, real or complex, shape (N,N) The unit cell Hamiltonian
        of the lead unit cell.
612
    h_hop : numpy array, real or complex, shape (N,M)
613
        The hopping matrix from a lead cell to the one on which self-energy
614
615
        has to be calculated (and any other hopping in the same direction).
    tol : float
616
617
        Numbers are considered zero when they are smaller than `tol` times
        the machine precision.
618
619
620
621

    Returns
    -------
    Sigma : numpy array, real or complex, shape (M,M)
622
623
624
        The computed self-energy. Note that even if `h_cell` and `h_hop` are
        both real, `Sigma` will typically be complex. (More precisely, if there
        is a propagating mode, `Sigma` will definitely be complex.)
625
626
627
628
629
630

    Notes
    -----
    For simplicity this function internally calculates the modes first.
    This may cause a small slowdown, and can be improved if necessary.
    """
631
632
    propagating, stabilized = modes(h_cell, h_hop, tol)
    return stabilized.selfenergy()
633
634


635
def square_selfenergy(width, hopping, fermi_energy):
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
    """
    Calculate analytically the self energy for a square lattice.

    The lattice is assumed to have a single orbital per site and
    nearest-neighbor hopping.

    Parameters
    ----------
    width : integer
        width of the lattice
    """

    # Following appendix C of M. Wimmer's diploma thesis:
    # http://www.physik.uni-regensburg.de/forschung/\
    # richter/richter/media/research/publications2004/wimmer-Diplomarbeit.pdf

652
    # p labels transversal modes.  i and j label the sites of a cell.
653
654
655
656
657
658
659
660
661
662
663
664

    # Precalculate the transverse wave function.
    psi_p_i = np.empty((width, width))
    factor = pi / (width + 1)
    prefactor = sqrt(2 / (width + 1))
    for p in xrange(width):
        for i in xrange(width):
            psi_p_i[p, i] = prefactor * sin(factor * (p + 1) * (i + 1))

    # Precalculate the integrals of the longitudinal wave functions.
    def f(q):
        if abs(q) <= 2:
Anton Akhmerov's avatar
Anton Akhmerov committed
665
            return q/2 - 1j * sqrt(1 - (q / 2) ** 2)
666
        else:
Anton Akhmerov's avatar
Anton Akhmerov committed
667
            return q/2 - copysign(sqrt((q / 2) ** 2 - 1), q)
668
669
670
    f_p = np.empty((width,), dtype=complex)
    for p in xrange(width):
        e = 2 * hopping * (1 - cos(factor * (p + 1)))
Anton Akhmerov's avatar
Anton Akhmerov committed
671
        q = (fermi_energy - e) / hopping - 2
672
673
674
675
676
677
678
679
680
        f_p[p] = f(q)

    # Put everything together into the self energy and return it.
    result = np.empty((width, width), dtype=complex)
    for i in xrange(width):
        for j in xrange(width):
            result[i, j] = hopping * \
                (psi_p_i[:, i] * psi_p_i[:, j].conj() * f_p).sum()
    return result