diff --git a/TODO b/TODO
index 72de50c890af3cbe6a8abc1ba67910b7b3aad14a..c3b3e7ad622d68585941e3cd739499b80a915b8c 100644
--- a/TODO
+++ b/TODO
@@ -4,8 +4,6 @@ Roughly in order of importance. -*-org-*-
* Make kwant run on windows
-* Implement "lead freezing"
-
* Add calculation of current density
* Consider making the b parameter of _solve_linear_sys a matrix instead of a
diff --git a/doc/other/linear_system.pdf b/doc/other/linear_system.pdf
deleted file mode 100644
index 384562b861868e7d3b961721c4a74698f6837c7f..0000000000000000000000000000000000000000
Binary files a/doc/other/linear_system.pdf and /dev/null differ
diff --git a/doc/other/reduce_modes.pdf b/doc/other/reduce_modes.pdf
deleted file mode 100644
index a01aad6df8305badad3a88190589bfb5608c26a2..0000000000000000000000000000000000000000
Binary files a/doc/other/reduce_modes.pdf and /dev/null differ
diff --git a/doc/source/reference/kwant.builder.rst b/doc/source/reference/kwant.builder.rst
index 5d0d063c1ff23b958f0fd57fdd2698e9e0cef601..973b8f5a3f06a4fc3dfa5f2e5df794eb75a3e92e 100644
--- a/doc/source/reference/kwant.builder.rst
+++ b/doc/source/reference/kwant.builder.rst
@@ -14,6 +14,7 @@ Types
SimpleSiteFamily
BuilderLead
SelfEnergy
+ ModesLead
Abstract base classes
---------------------
diff --git a/doc/source/reference/kwant.physics.rst b/doc/source/reference/kwant.physics.rst
index 90128888cb7aa734a978a45e3f6ed0d099d48ff5..7b202da06aba3f9d78badfc4ba20796c34a52320 100644
--- a/doc/source/reference/kwant.physics.rst
+++ b/doc/source/reference/kwant.physics.rst
@@ -10,4 +10,6 @@ Leads
Bands
modes
- self_energy
+ selfenergy
+ selfenergy_from_modes
+ ModesTuple
diff --git a/doc/source/whatsnew/0.3.rst b/doc/source/whatsnew/0.3.rst
index f6dad072149c807207fda73d4ca9450eb6047e73..84b0da4b49d9fdd7daa39c7da11980e7aa0f52b9 100644
--- a/doc/source/whatsnew/0.3.rst
+++ b/doc/source/whatsnew/0.3.rst
@@ -47,6 +47,8 @@ Some renames
------------
* site groups are now called site families. This affects all the names that
used to contain "group" or "groups".
+* ``self_energy`` has been renamed to ``selfenergy`` in all cases, most notably
+ in `kwant.physics.selfenergy`.
* ``wave_func`` has been renamed to `~kwant.solvers.default.wave_function`,
* ``MonatomicLattice`` has been renamed to `~kwant.lattice.Monatomic`,
* ``PolyatomicLattice`` has been renamed to `~kwant.lattice.Polyatomic`.
@@ -138,3 +140,19 @@ collection could be a dictionary, or a class instance, for example::
Arguments can be passed in an equivalent way to
`~kwant.solvers.default.wave_function`,
`~kwant.system.System.hamiltonian_submatrix`, etc.
+
+Calculation of modes separated from solving
+-------------------------------------------
+The interface that solvers expect from leads attached to a
+`~kwant.system.FiniteSystem` has been simplified and codified (see there).
+Similar to self-energy, calculation of modes is now the lead's own
+responsibility.
+
+The new class `~kwant.builder.ModesLead` allows to attach leads that have a
+custom way of calculating their modes (e.g. ideal leads) directly to a
+`~kwant.builder.Builder`.
+
+Modes or self-energies can now be precomputed before passing the system to a
+solver, using the method `~kwant.system.FiniteSystem.precalculate`. This may
+save time, when the linear system has to be solved many times with the same
+lead parameters.
diff --git a/examples/square.py b/examples/square.py
index d390d9c5f28e630906280ab3329dfd9c3054148c..19224461d04ff3b997613c0d6e8bff3485910577 100644
--- a/examples/square.py
+++ b/examples/square.py
@@ -6,7 +6,7 @@ from __future__ import division
import numpy as np
from matplotlib import pyplot
import kwant
-from kwant.physics.selfenergy import square_self_energy
+from kwant.physics.selfenergy import square_selfenergy
__all__ = ['System']
@@ -16,8 +16,13 @@ class Lead(object):
self.t = t
self.potential = potential
+<<<<<<< HEAD
def self_energy(self, fermi_energy, args=()):
return square_self_energy(self.width, self.t, self.potential,
+=======
+ def selfenergy(self, fermi_energy):
+ return square_selfenergy(self.width, self.t, self.potential,
+>>>>>>> modes
fermi_energy)
class System(kwant.system.FiniteSystem):
diff --git a/examples/tests/test_square.py b/examples/tests/test_square.py
index e49d67fffcd208e99c58b5c468c349f4d956a9e6..9abfbbb1aad05377e33f9d0bbeb491369faf7fa3 100644
--- a/examples/tests/test_square.py
+++ b/examples/tests/test_square.py
@@ -26,12 +26,12 @@ def test_hamiltonian():
assert_almost_equal(m, m_herm)
assert_almost_equal(m_herm, sys.hamiltonian(j, i))
-def test_self_energy():
+def test_selfenergy():
sys = square.System((2, 4), 1)
for lead in xrange(len(sys.lead_interfaces)):
n_orb = sum(
sys.num_orbitals(site) for site in sys.lead_interfaces[lead])
- se = sys.self_energy(lead, 0)
+ se = sys.selfenergy(lead, 0)
assert_equal(len(se.shape), 2)
assert_equal(se.shape[0], se.shape[1])
assert_equal(se.shape[0], n_orb)
diff --git a/kwant/builder.py b/kwant/builder.py
index 82462dc9dd53251914cbc8d09dfb15167ddbc7e7..b6346e644e7635d9ae55ccc1b32f4af2eb86bc57 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -9,7 +9,7 @@
from __future__ import division
__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
- 'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergy']
+ 'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergy', 'ModesLead']
import abc
import sys
@@ -18,7 +18,7 @@ import operator
from itertools import izip, islice, chain
import tinyarray as ta
import numpy as np
-from . import system, graph
+from . import system, graph, physics
################ Sites and site families
@@ -387,15 +387,12 @@ class Lead(object):
Notes
-----
- The finalized lead must at least have a single method
- ``self_energy(energy)`` but it can be a full
- `kwant.system.InfiniteSystem` as well.
+ The finalized lead must be an object that can be used as a lead in a
+ `kwant.system.FiniteSystem`. It could be an instance of
+ `kwant.system.InfiniteSystem` for example.
- The method ``self_energy`` of the finalized lead must return a square
- matrix of appropriate size.
-
- The order of interface sites is assumed to be preserved during
- finalization.
+ The order of sites for the finalized lead must be the one specified in
+ `interface`.
"""
pass
@@ -444,21 +441,49 @@ class SelfEnergy(Lead):
Parameters
----------
- self_energy_func : function
+ selfenergy_func : function
Function which returns the self energy matrix for the interface sites
given the energy and optionally a list of extra arguments.
interface : sequence of `Site` instances
"""
- def __init__(self, self_energy_func, interface):
- self.self_energy_func = self_energy_func
+ def __init__(self, selfenergy_func, interface):
+ self._selfenergy_func = selfenergy_func
+ self.interface = tuple(interface)
+
+ def finalized(self):
+ """Trivial finalization: the object is returned itself."""
+ return self
+
+ def selfenergy(self, energy, args=()):
+ return self._selfenergy_func(energy, args)
+
+
+class ModesLead(Lead):
+ """A general lead defined by its modes wave functions.
+
+ Parameters
+ ----------
+ modes_func : function
+ Function which returns the modes of the lead in the format specified
+ in `~kwant.physics.ModesTuple` given the energy and optionally a list of
+ extra arguments.
+ interface : sequence of `Site` instances
+ """
+ def __init__(self, modes_func, interface):
+ self._modes_func = modes_func
self.interface = tuple(interface)
def finalized(self):
"""Trivial finalization: the object is returned itself."""
return self
- def self_energy(self, energy, args=()):
- return self.self_energy_func(energy, args)
+ def modes(self, energy, args=()):
+ return self._modes_func(energy, args)
+
+ def selfenergy(self, energy, args=()):
+ modes = self.modes(energy, args)
+ return physics.selfenergy_from_modes(modes=modes)
+
################ Builder class
diff --git a/kwant/physics/__init__.py b/kwant/physics/__init__.py
index 2b081aa8ddde4a367c3f7a1344c0865a360e77fe..cf6a2c5b010caf97bcd58e45450c6c389dd047b0 100644
--- a/kwant/physics/__init__.py
+++ b/kwant/physics/__init__.py
@@ -10,7 +10,7 @@
# Merge the public interface of all submodules.
__all__ = []
-for module in ['selfenergy', 'dispersion']:
+for module in ['leads', 'dispersion']:
exec 'from . import {0}'.format(module)
exec 'from .{0} import *'.format(module)
exec '__all__.extend({0}.__all__)'.format(module)
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
new file mode 100644
index 0000000000000000000000000000000000000000..d417fc4cd5922b3ef25b7cc7322069fa9972128f
--- /dev/null
+++ b/kwant/physics/leads.py
@@ -0,0 +1,670 @@
+# Copyright 2011-2013 kwant authors.
+#
+# This file is part of kwant. It is subject to the license terms in the
+# LICENSE file found in the top-level directory of this distribution and at
+# http://kwant-project.org/license. A list of kwant authors can be found in
+# the AUTHORS file at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from __future__ import division
+from math import sin, cos, sqrt, pi, copysign
+from collections import namedtuple
+from itertools import izip
+import numpy as np
+import numpy.linalg as npl
+import scipy.linalg as la
+from .. import linalg as kla
+
+dot = np.dot
+
+__all__ = ['selfenergy', 'modes', 'ModesTuple', 'selfenergy_from_modes']
+
+Linsys = namedtuple('Linsys', ['eigenproblem', 'v', 'extract', 'project'])
+
+def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
+ """
+ Make an eigenvalue problem for eigenvectors of translation operator.
+
+ Parameters
+ ----------
+ h_onslice : numpy array with shape (n, n)
+ Hamiltonian of a single lead slice.
+ h_hop : numpy array with shape (n, m), m <= n
+ Hopping Hamiltonian from the slice to the next one.
+ tol : float
+ Numbers are considered zero when they are smaller than `tol` times
+ the machine precision.
+ algorithm : tuple of 3 booleans or None
+ Which steps of the eigenvalue prolem stabilization to perform.
+ The first value selects whether to work in the basis of the
+ hopping svd, or lattice basis. If the real space basis is chosen, the
+ following two options do not apply.
+ The second value selects whether to add an anti-Hermitian term to the
+ slice Hamiltonian before inverting. Finally the third value selects
+ whether to reduce a generalized eigenvalue problem to the regular one.
+ The default value, `None`, results in kwant selecting the algorithm
+ that is the most efficient without sacrificing stability. Manual
+ selection may result in either slower performance, or large numerical
+ errors, and is mostly required for testing purposes.
+
+ Returns
+ -------
+ linsys : namedtuple
+ 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.
+
+ Notes
+ -----
+ The lead problem with degenerate hopping is rather complicated, and the
+ details of the algorithm will be published elsewhere.
+ """
+ n = h_onslice.shape[0]
+ m = h_hop.shape[1]
+
+ if not (np.any(h_hop.real) or np.any(h_hop.imag)):
+ # Inter-slice hopping is zero. The current algorithm is not suited to
+ # treat this extremely singular case.
+ # Note: np.any(h_hop) returns (at least from numpy 1.6.*)
+ # False if h_hop is purely imaginary
+ raise ValueError("Inter-slice hopping is exactly zero.")
+
+ eps = np.finfo(np.common_type(h_onslice, h_hop)).eps * tol
+
+ # 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.)
+
+ u, s, vh = la.svd(h_hop)
+ assert m == vh.shape[1], "Corrupt output of svd."
+ n_nonsing = np.sum(s > eps * s[0])
+
+ if (n_nonsing == n and algorithm is None) or (algorithm is not None and
+ algorithm[0]):
+ # The hopping matrix is well-conditioned and can be safely inverted.
+ # Hence the regular transfer matrix may be used.
+ hop_inv = la.inv(h_hop)
+
+ A = np.zeros((2*n, 2*n), dtype=np.common_type(h_onslice, h_hop))
+ A[:n, :n] = dot(hop_inv, -h_onslice)
+ A[:n, n:] = -hop_inv
+ A[n:, :n] = h_hop.T.conj()
+
+ # The function that can extract the full wave function psi from the
+ # projected one. Here it is almost trivial, but used for simplifying
+ # the logic.
+
+ def extract_wf(psi, lmbdainv):
+ return psi[:n]
+
+ # Project the full wave function back (also trivial).
+
+ def project_wf(psi, lmbdainv):
+ return np.r_[psi * lmbdainv, dot(h_hop.T.conj(), psi)]
+
+ matrices = (A, None)
+ v_out = None
+ else:
+ if algorithm is not None:
+ need_to_stabilize, divide = algorithm[1:]
+ else:
+ need_to_stabilize = None
+
+ # 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.
+ u = u[:, :n_nonsing]
+ s = s[:n_nonsing]
+ u = u * np.sqrt(s)
+ # pad v with zeros if necessary
+ v = np.zeros((n, n_nonsing), dtype=vh.dtype)
+ v[:vh.shape[1]] = vh[:n_nonsing].T.conj()
+ v = v * np.sqrt(s)
+
+ # 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 - as 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
+
+ h = h_onslice
+ sol = kla.lu_factor(h)
+ if issubclass(np.common_type(h_onslice, h_hop), np.floating) \
+ and need_to_stabilize is None:
+ # Check if stabilization is needed.
+ rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))
+
+ if rcond > eps:
+ need_to_stabilize = False
+
+ if need_to_stabilize is None:
+ need_to_stabilize = True
+
+ if need_to_stabilize:
+ # Matrices are complex or need self-energy-like term to be
+ # stabilized.
+ temp = dot(u, u.T.conj()) + dot(v, v.T.conj())
+ h = h_onslice + 1j * temp
+
+ 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.
+ if rcond < eps:
+ raise RuntimeError("Flat band encountered at the requested "
+ "energy, result is badly defined.")
+
+ # The function that can extract the full wave function psi from
+ # the projected one (v^dagger psi lambda^-1, s u^dagger psi).
+
+ def extract_wf(psi, lmbdainv):
+ wf = - dot(u, psi[: n_nonsing] * lmbdainv) - \
+ 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)
+
+ # Project the full wave function back.
+
+ def project_wf(psi, lmbdainv):
+ return np.r_[dot(v.T.conj(), psi * lmbdainv),
+ dot(u.T.conj(), psi)]
+
+ # Setup the generalized eigenvalue problem.
+
+ 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))
+
+ begin, end = slice(n_nonsing), slice(n_nonsing, None)
+
+ A[end, begin] = np.identity(n_nonsing)
+ temp = kla.lu_solve(sol, v)
+ temp2 = dot(u.T.conj(), temp)
+ if need_to_stabilize:
+ A[begin, begin] = -1j * temp2
+ A[begin, end] = temp2
+ temp2 = dot(v.T.conj(), temp)
+ if need_to_stabilize:
+ A[end, begin] -= 1j *temp2
+ A[end, end] = temp2
+
+ B[begin, end] = -np.identity(n_nonsing)
+ temp = kla.lu_solve(sol, u)
+ temp2 = dot(u.T.conj(), temp)
+ B[begin, begin] = -temp2
+ if need_to_stabilize:
+ B[begin, end] += 1j * temp2
+ temp2 = dot(v.T.conj(), temp)
+ B[end, begin] = -temp2
+ if need_to_stabilize:
+ B[end, end] = 1j * temp2
+
+ v_out = v[:m]
+
+ # 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)
+ if algorithm is None:
+ 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.
+ divide = rcond > eps * tol
+
+ if divide:
+ matrices = (kla.lu_solve(lu_b, A), None)
+ else:
+ matrices = (A, B)
+ return Linsys(matrices, v_out, extract_wf, project_wf)
+
+
+def make_proper_modes(lmbdainv, psi, extract, project, tol=1e6):
+ """
+ Determine the velocities and direction of the propagating eigenmodes.
+
+ 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.
+ """
+ vel_eps = np.finfo(psi.dtype).eps * tol
+
+ # h_hop is either the full hopping matrix, or the singular
+ # values vector of the svd.
+ nmodes = psi.shape[1]
+ n = len(psi) // 2
+
+ # Array for the velocities.
+ velocities = np.empty(nmodes, dtype=float)
+
+ # Mark the right-going modes.
+ rightselect = np.zeros(nmodes, dtype=bool)
+
+ # Find clusters of nearby points.
+ 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
+
+ for interval in izip(boundaries[:-1], boundaries[1:]):
+ if interval[1] > boundaries[0] + len(angles):
+ break
+
+ indx = sort_order[interval[0] : interval[1]]
+
+ # 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).
+
+ # Finding the true modes is done in two steps:
+
+ # 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:
+
+ if len(indx) > 1:
+ full_psi = extract(psi[:, indx], lmbdainv[indx])
+
+ # Note: Here's a workaround for the fact that the interface
+ # to qr changed from SciPy 0.8.0 to 0.9.0
+ try:
+ full_psi = la.qr(full_psi, econ=True, mode='qr')[0]
+ except TypeError:
+ full_psi = la.qr(full_psi, mode='economic')[0]
+
+ psi[:, indx] = project(full_psi, lmbdainv[indx])
+
+ # 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))
+ psi[:, indx] = dot(psi[:, indx], rot)
+ velocities[indx] = vel_vals
+
+ rightselect = velocities > vel_eps
+ if np.any(abs(velocities) < vel_eps):
+ raise RuntimeError("Found a mode with zero or close to zero velocity.")
+ if 2 * np.sum(rightselect) != len(velocities):
+ raise RuntimeError("Numbers of left- and right-propagating "
+ "modes differ, possibly due to a numerical "
+ "instability.")
+
+ return psi, velocities, rightselect
+
+
+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.
+
+ 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.
+
+ 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_onslice.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.
+ """
+ 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
+
+
+_Modes = namedtuple('ModesTuple', ['vecs', 'vecslmbdainv', 'nmodes',
+ 'sqrt_hop'])
+modes_docstring = """Eigendecomposition of a translation operator.
+
+A named tuple with attributes `(vecs, vecslmbdainv, nmodes, sqrt_hop)`. 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
+with 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)`. `vecs` and
+`vecslmbdainv` are the first and the second halves of the wave functions. The
+first `nmodes` are eigenmodes moving in the negative direction (hence they are
+incoming into the system in kwant convention), the second `nmodes` are
+eigenmodes moving in the positive direction. The remaining modes are 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 he `sqrt_hop` attribute is `v
+sqrt(s)`.
+"""
+d = dict(_Modes.__dict__)
+d.update(__doc__=modes_docstring)
+ModesTuple = type('ModesTuple', _Modes.__bases__, d)
+del _Modes, d, modes_docstring
+
+def modes(h_onslice, h_hop, tol=1e6, algorithm=None):
+ """
+ Compute the eigendecomposition of a translation operator of a lead.
+
+ Parameters
+ ----------
+ h_onslice : numpy array, real or complex, shape (N,N) The unit cell
+ Hamiltonian of the lead slice.
+ h_hop : numpy array, real or complex, shape (N,M)
+ The hopping matrix from a lead slice to the one on which self-energy
+ has to be calculated (and any other hopping in the same direction).
+ tol : float
+ Numbers and differences are considered zero when they are smaller
+ than `tol` times the machine precision.
+ algorithm : tuple or 3 booleans or None
+ Which steps of the eigenvalue prolem stabilization to perform. The
+ default value, `None`, results in kwant selecting the algorithm that is
+ the most efficient without sacrificing stability. Manual selection may
+ result in either slower performance, or large numerical errors, and is
+ mostly required for testing purposes. The first value selects whether
+ to work in the basis of the hopping svd, or lattice basis. If the real
+ space basis is chosen, the following two options do not apply. The
+ second value selects whether to add an anti-Hermitian term to the slice
+ Hamiltonian before inverting. Finally the third value selects whether
+ to reduce a generalized eigenvalue problem to the regular one.
+
+
+ Returns
+ -------
+ ModesTuple(vecs, vecslmbdainv, nmodes, v) : a named tuple
+ See notes below and `~kwant.physics.ModesTuple`
+ documentation for details.
+
+ Notes
+ -----
+ Only the propagating modes and the modes decaying away from the system are
+ returned: The first `nmodes` columns in `vecs` correspond to incoming modes
+ (coming from the lead into the system), the following `nmodes` columns
+ correspond to outgoing modes (going into the lead, away from the system),
+ the remaining columns are evanescent modes, decaying away from the system.
+
+ 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.
+
+ If `h_hop` is invertible, the full transverse wave functions are returned.
+ If it is singular, the projections (s u^dagger psi, v^dagger psi lambda^-1)
+ are returned.
+
+ In order for the linear system to be well-defined, instead of the
+ evanescent modes, an orthogonal basis in the space of evanescent modes is
+ returned.
+
+ Propagating modes with the same lambda are orthogonalized. All the
+ propagating modes are normalized by current.
+
+ This function uses the most stable and efficient algorithm for calculating
+ the mode decomposition. Its details will be published elsewhere.
+
+ """
+ m = h_hop.shape[1]
+
+ if (h_onslice.shape[0] != h_onslice.shape[1] or
+ h_onslice.shape[0] != h_hop.shape[0]):
+ raise ValueError("Incompatible matrix sizes for h_onslice and h_hop.")
+
+ # Note: np.any(h_hop) returns (at least from numpy 1.6.1 - 1.8-devel)
+ # False if h_hop is purely imaginary
+ if not (np.any(h_hop.real) or np.any(h_hop.imag)):
+ n = h_hop.shape[0]
+ v = np.empty((0, m))
+ return ModesTuple(np.empty((0, 0)), np.empty((0, 0)), 0, v)
+
+ # Defer most of the calculation to helper routines.
+ matrices, v, extract, project = setup_linsys(h_onslice, h_hop, tol,
+ algorithm)
+ ev, evanselect, propselect, vec_gen, ord_schur =\
+ unified_eigenproblem(*(matrices + (tol,)))
+
+ if v is not None:
+ n = v.shape[1]
+ else:
+ n = h_onslice.shape[0]
+
+ nprop = np.sum(propselect)
+ nevan = n - nprop // 2
+ evanselect_bool = np.zeros((2*n), dtype='bool')
+ evanselect_bool[evanselect] = True
+ evan_vecs = ord_schur(evanselect)[:, :nevan]
+
+ # Compute the propagating eigenvectors.
+ prop_vecs = vec_gen(propselect)
+ # Compute their velocity, and, if necessary, rotate them
+ prop_vecs, vel, rprop = \
+ make_proper_modes(ev[propselect], prop_vecs, extract, project, tol)
+
+ # Normalize propagating eigenvectors by velocities.
+ prop_vecs /= np.sqrt(abs(vel))
+
+ # Fix the phase factor - make maximum of the transverse wave function real
+ # TODO (Anton): Take care of multiple maxima when normalizing.
+ maxnode = prop_vecs[n + np.argmax(abs(prop_vecs[n:, :]), axis=0),
+ np.arange(prop_vecs.shape[1])]
+ maxnode /= abs(maxnode)
+ prop_vecs /= maxnode
+
+ lprop = np.logical_not(rprop)
+ nmodes = np.sum(rprop)
+
+ # Sort modes according to their k-vector (1/lambda = e^{-ik}):
+ # - right-going modes: sort that k is in ascending order
+ # - left-going modes: sort that k is in descending order
+ # (note that k can be positive or negative). With this convention,
+ # the modes of a simple square lattice strip are ordered as
+ # expected (lowest subband first, etc.)
+
+ prop_ev = ev[propselect]
+ rsort = np.argsort((-1j * np.log(prop_ev[rprop])).real)
+ lsort = np.argsort((1j * np.log(prop_ev[lprop])).real)
+
+ # The following is necessary due to how numpy deals with indexing of empty.
+ # arrays.
+ if nmodes == 0:
+ lprop = rprop = rsort = lsort = slice(None)
+
+ vecs = np.c_[prop_vecs[n:, lprop][:, lsort],
+ prop_vecs[n:, rprop][:, rsort],
+ evan_vecs[n:]]
+ vecslmbdainv = np.c_[prop_vecs[:n, lprop][:, lsort],
+ prop_vecs[:n, rprop][:, rsort],
+ evan_vecs[:n]]
+
+ return ModesTuple(vecs, vecslmbdainv, nmodes, v)
+
+
+def selfenergy_from_modes(lead_modes):
+ """
+ Compute the self-energy generated by lead modes.
+
+ Parameters
+ ----------
+ lead_modes : ModesTuple(vecs, vecslmbdainv, nmodes, v) a named tuple
+ The modes in the lead, with format defined in
+ `~kwant.physics.ModesTuple`.
+
+ Returns
+ -------
+ Sigma : numpy array, real or complex, shape (M,M)
+ The computed self-energy. Note that even if `h_onslice` and `h_hop`
+ are both real, `Sigma` will typically be complex. (More precisely, if
+ there is a propagating mode, `Sigma` will definitely be complex.)
+
+ Notes
+ -----
+ For simplicity this function relies on the calculation of modes as input.
+ This may cause a small slowdown, and can be improved if necessary.
+ """
+ vecs, vecslmbdainv, nmodes, v = lead_modes
+ vecs = vecs[:, nmodes:]
+ vecslmbdainv = vecslmbdainv[:, 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
+
+
+def selfenergy(h_onslice, h_hop, tol=1e6):
+ """
+ Compute the self-energy generated by the lead.
+
+ Parameters
+ ----------
+ h_onslice : numpy array, real or complex, shape (N,N) The unit cell
+ Hamiltonian of the lead slice.
+ h_hop : numpy array, real or complex, shape (N,M)
+ The hopping matrix from a lead slice to the one on which self-energy
+ has to be calculated (and any other hopping in the same direction).
+ tol : float
+ Numbers are considered zero when they are smaller than `tol` times
+ the machine precision.
+
+ Returns
+ -------
+ Sigma : numpy array, real or complex, shape (M,M)
+ The computed self-energy. Note that even if `h_onslice` and `h_hop`
+ are both real, `Sigma` will typically be complex. (More precisely, if
+ there is a propagating mode, `Sigma` will definitely be complex.)
+
+ Notes
+ -----
+ For simplicity this function internally calculates the modes first.
+ This may cause a small slowdown, and can be improved if necessary.
+ """
+ return selfenergy_from_modes(modes(h_onslice, h_hop, tol))
+
+
+def square_selfenergy(width, hopping, fermi_energy):
+ """
+ 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
+
+ # p labels transversal modes. i and j label the sites of a slice.
+
+ # 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:
+ return q/2 - 1j * sqrt(1 - (q / 2) ** 2)
+ else:
+ return q/2 - copysign(sqrt((q / 2) ** 2 - 1), q)
+ f_p = np.empty((width,), dtype=complex)
+ for p in xrange(width):
+ e = 2 * hopping * (1 - cos(factor * (p + 1)))
+ q = (fermi_energy - e) / hopping - 2
+ 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
diff --git a/kwant/physics/selfenergy.py b/kwant/physics/selfenergy.py
deleted file mode 100644
index e741a66193ba99cf654c02001d99a4cee336178a..0000000000000000000000000000000000000000
--- a/kwant/physics/selfenergy.py
+++ /dev/null
@@ -1,795 +0,0 @@
-# Copyright 2011-2013 kwant authors.
-#
-# This file is part of kwant. It is subject to the license terms in the
-# LICENSE file found in the top-level directory of this distribution and at
-# http://kwant-project.org/license. A list of kwant authors can be found in
-# the AUTHORS file at the top-level directory of this distribution and at
-# http://kwant-project.org/authors.
-
-from __future__ import division
-from math import sin, cos, sqrt, pi, copysign
-from collections import namedtuple
-import numpy as np
-import numpy.linalg as npl
-import scipy.linalg as la
-from .. import linalg as kla
-
-dot = np.dot
-
-__all__ = ['self_energy', 'modes', 'Modes']
-
-
-def setup_linsys(h_onslice, h_hop, tol=1e6):
- """
- Make an eigenvalue problem for eigenvectors of translation operator.
-
- Parameters
- ----------
- h_onslice : NumPy array with shape (n, n)
- Hamiltonian of a single lead slice.
- h_hop : NumPy array with shape (n, m), m <= n
- Hopping Hamiltonian from the slice to the next one.
-
- Returns
- -------
- linsys : matrix or tuple
- if the hopping is nonsingular, a single matrix defining an eigenvalue
- problem is returned, othewise a tuple of two matrices defining a
- generalized eigenvalue problem together with additional information is
- returned.
-
- Notes
- -----
- The lead problem with degenerate hopping is rather complicated, and it is
- described in kwant/doc/other/lead_modes.pdf.
- """
- n = h_onslice.shape[0]
- m = h_hop.shape[1]
-
- # Inter-slice hopping is zero. The current algorithm is not suited to
- # treat this extremely singular case.
- # Note: np.any(h_hop) returns (at least from NumPy 1.6.1 - 1.8-devel)
- # False if h_hop is purely imaginary
- assert np.any(h_hop.real) or np.any(h_hop.imag)
-
- eps = np.finfo(np.common_type(h_onslice, h_hop)).eps
-
- # First check if the hopping matrix has eigenvalues close to 0.
- u, s, vh = la.svd(h_hop)
-
- assert m == vh.shape[1], "Corrupt output of svd."
-
- # Count the number of singular values close to zero.
- # (Close to zero is defined here as |x| < eps * tol * s[0] , where
- # s[0] is the largest singular value.)
- n_nonsing = np.sum(s > eps * tol * s[0])
-
- if n_nonsing == n:
- # The hopping matrix is well-conditioned and can be safely inverted.
- sol = kla.lu_factor(h_hop)
-
- A = np.empty((2*n, 2*n), dtype=np.common_type(h_onslice, h_hop))
-
- A[0:n, 0:n] = kla.lu_solve(sol, -h_onslice)
- A[0:n, n:2*n] = kla.lu_solve(sol, -h_hop.T.conj())
- A[n:2*n, 0:n] = np.identity(n)
- A[n:2*n, n:2*n] = 0
-
- return A
- else:
- # 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.
- u = u[:, :n_nonsing]
- s = s[:n_nonsing]
- # pad v with zeros if necessary
- v = np.zeros((n, n_nonsing), dtype=vh.dtype)
- v[:vh.shape[1], :] = vh[:n_nonsing, :].T.conj()
-
- # 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 - as 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
-
- gamma = None
-
- if issubclass(np.common_type(h_onslice, h_hop), np.floating):
-
- # Check if stabilization is needed.
- h = h_onslice
-
- sol = kla.lu_factor(h)
- rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))
-
- if rcond > eps * tol:
- gamma = 0
-
- if gamma is None:
- # Matrices are complex or need self-energy-like term to be
- # stabilized.
-
- # Normalize such that the maximum entry in the
- # self-energy-like term has a value comparable to the
- # maximum entry in h_onslice.
-
- temp = dot(u, u.T.conj()) + dot(v, v.T.conj())
-
- max_h = np.amax(abs(h_onslice))
- max_temp = np.amax(abs(temp))
-
- gamma = max_h / max_temp * 1j
-
- h = h_onslice + gamma * temp
-
- 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.
- if rcond < eps * tol:
- raise RuntimeError("Flat band encountered at the requested "
- "energy, result is badly defined.")
-
- # Function that can extract the full wave function psi from
- # the projected one (v^dagger psi lambda^-1, u^dagger psi).
-
- def extract_wf(psi, lmbdainv):
- return kla.lu_solve(sol,
- gamma * dot(v, psi[: n_nonsing]) +
- gamma * dot(u, psi[n_nonsing:] * lmbdainv) -
- dot(u * s, psi[: n_nonsing] * lmbdainv) -
- dot(v * s, psi[n_nonsing:]))
-
- # Project a full wave function back.
-
- def project_wf(psi, lmbdainv):
- return np.asarray(np.bmat([[dot(v.T.conj(), psi * lmbdainv)],
- [dot(u.T.conj(), psi)]]))
-
- # Setup the generalized eigenvalue problem.
-
- A = np.empty((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
- B = np.empty((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
-
- A[:n_nonsing, :n_nonsing] = -np.eye(n_nonsing)
-
- B[n_nonsing: 2 * n_nonsing,
- n_nonsing: 2 * n_nonsing] = np.eye(n_nonsing)
-
- temp = kla.lu_solve(sol, v)
- temp2 = dot(u.T.conj(), temp)
- A[n_nonsing : 2 * n_nonsing, :n_nonsing] = gamma * temp2
- A[n_nonsing : 2 * n_nonsing, n_nonsing: 2 * n_nonsing] = - temp2 * s
- temp2 = dot(v.T.conj(), temp)
- A[:n_nonsing, :n_nonsing] += gamma * temp2
- A[:n_nonsing, n_nonsing : 2 * n_nonsing] = - temp2 * s
-
- temp = kla.lu_solve(sol, u)
- temp2 = dot(u.T.conj(), temp)
- B[n_nonsing : 2 * n_nonsing, :n_nonsing] = temp2 * s
- B[n_nonsing : 2 * n_nonsing, n_nonsing : 2 * n_nonsing] -= \
- gamma * temp2
- temp2 = dot(v.T.conj(), temp)
- B[:n_nonsing, :n_nonsing] = temp2 * s
- B[:n_nonsing, n_nonsing : 2 * n_nonsing] = - gamma * temp2
-
- # 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)
- rcond = kla.rcond_from_lu(lu_b, npl.norm(B, 1))
-
- # I put a more stringent condition here - errors can accumulate
- # from here to the eigenvalue calculation later.
- if rcond > eps * tol**2:
- return (kla.lu_solve(lu_b, A), (u, s, v[:m, :]),
- (extract_wf, project_wf))
- else:
- return (A, B, (u, s, v[:m]), (extract_wf, project_wf))
-
-
-def split_degenerate(evs, tol=1e6):
- """
- Find sets of approximately degenerate list elements on a unit circle.
-
- Given a list of eigenvalues on the unit circle, return a list containing
- tuples of indices of eigenvalues that are numerically degenerate. Two
- eigenvalues ev[i] and ev[j] are considered to be numerically degenerate if
- abs(ev[i] - ev[j]) < eps * tol, where eps is the machine precision.
-
- Example
- -------
- >>> split_degenerate(np.array([1,-1,1,1], dtype=complex))
- [(1,), (0, 2, 3)].
- """
- eps = np.finfo(evs.dtype).eps
-
- n = evs.size
- evlist = []
-
- # Figure out if there are degenerate eigenvalues.
- # For this, sort according to k, which is i*log(ev) (ev is exp(-ik)).
- k = np.log(evs).imag
- sortindx = np.argsort(k)
- evs_sorted = evs[sortindx]
-
- # Note that we sorted eigenvalues on the unit circle, cutting
- # the unit circle at -1. We thus must search for degeneracies also
- # across this cut
-
- start = 0
- while (start - 1 > -n and
- abs(evs_sorted[start - 1] - evs_sorted[start]) < eps * tol):
- start = start - 1
-
- stop = n + start
-
- while start < stop:
- deglist = [sortindx[start]]
- while (start + 1 < stop and
- abs(evs_sorted[start] - evs_sorted[start + 1]) < eps * tol):
- start += 1
- deglist.append(sortindx[start])
-
- evlist.append(tuple(deglist))
- start += 1
-
- return evlist
-
-
-def make_proper_modes(lmbdainv, psi, h_hop, extract=None,
- project=None, tol=1e6):
- """
- Determine the velocities and direction of the propagating eigenmodes.
-
- 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.
- """
-
- vel_eps = np.finfo(np.common_type(psi, h_hop)).eps * tol
-
- # h_hop is either the full hopping matrix, or the singular
- # values vector of the svd.
- if h_hop.ndim == 2:
- n = h_hop.shape[0]
- m = h_hop.shape[1]
- else:
- n = h_hop.size
-
- nmodes = psi.shape[1]
-
- if nmodes == 0:
- raise ValueError('Empty mode array.')
-
- # Array for the velocities.
- v = np.empty(nmodes, dtype=float)
-
- # Mark the right-going modes.
- rightselect = np.zeros(nmodes, dtype=bool)
-
- n_left = n_right = 0
- crossing = False
-
- indxclust = split_degenerate(lmbdainv)
-
- for indx in indxclust:
- if len(indx) > 1:
- # Several degenerate propagating modes. In this case, the computed
- # eigenvectors do not orthogonalize the velocity
- # operator, i.e. they do not have a proper velocity.
-
- indx = np.array(indx)
-
- # 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).
-
- if extract is not None:
- full_psi = extract(psi[:, indx], lmbdainv[indx])
- else:
- full_psi = psi[:n, indx]
-
- # Finding the true modes is done in two steps:
-
- # 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 are therefore orthogonalized:
-
- # Note: Here's a workaround for the fact that the interface
- # to qr changed from SciPy 0.8.0 to 0.9.0
- try:
- full_psi = la.qr(full_psi, econ=True, mode='qr')[0]
- except TypeError:
- full_psi = la.qr(full_psi, mode='economic')[0]
-
- if project:
- psi[:, indx] = project(full_psi, lmbdainv[indx])
- else:
- psi[:n, indx] = full_psi * lmbdainv[indx]
- psi[n:2*n, indx] = full_psi
-
- # 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,
- # as we are happy with any superposition in this case.
-
- if h_hop.ndim == 2:
- vel_op = -1j * dot(psi[n:, indx].T.conj(),
- dot(h_hop, psi[:m, indx]))
- else:
- vel_op = -1j * dot(psi[n:, indx].T.conj() * h_hop,
- 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,
- # will typically become complex after the rotation.
- if psi.dtype != np.common_type(psi, rot):
- psi = psi.astype(np.common_type(psi, rot))
-
- psi[:, indx] = dot(psi[:, indx], rot)
-
- v[indx] = vel_vals
-
- # For some of the self-energy methods it matters
- # whether the degeneracy is a crossing with velocities
- # of different sign
- if not ((vel_vals > 0).all() or (vel_vals < 0).all()):
- crossing = True
-
- for (vel, k) in zip(vel_vals, indx):
- if vel > vel_eps:
- n_right += 1
- rightselect[k] = True
- elif vel < -vel_eps:
- n_left += 1
- else:
- raise RuntimeError("Found a mode with zero or close to "
- "zero velocity.")
- else:
- # A single, unique propagating mode
- k = indx[0]
-
- if h_hop.ndim == 2:
- v[k] = 2 * dot(dot(psi[n:2*n, k:k+1].T.conj(), h_hop),
- psi[:m, k:k+1]).imag
- else:
- v[k] = 2 * dot(psi[n:2*n, k:k+1].T.conj() * h_hop,
- psi[0:n, k:k+1]).imag
-
- if v[k] > vel_eps:
- rightselect[k] = True
- n_right += 1
- elif v[k] < -vel_eps:
- n_left += 1
- else:
- raise RuntimeError("Found a mode with zero or close to "
- "zero velocity.")
-
- if n_left != n_right:
- raise RuntimeError("Numbers of left- and right-propagating "
- "modes differ.")
-
- return psi, v, rightselect, crossing
-
-
-def unified_eigenproblem(h_onslice, h_hop, tol):
- """A helper routine for general() and modes(), that wraps eigenproblems.
-
- This routine wraps the different types of eigenproblems that can arise
- in a unfied way.
-
- Returns
- -------
- ev : NumPy array
- an array of eigenvalues (can contain NaNs and Infs, but those
- are not accessed in `general()` and `modes()`) The number of
- eigenvalues is given by twice the number of nonzero singular values of
- `h_hop` (i.e. `2*h_onslice.shape[0]` if `h_hop` is invertible).
- evanselect : NumPy array
- index array of right-decaying modes.
- propselect : NumPy 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.
- u, w, v :
- if the hopping is singular, the svd of the hopping matrix, otherwise
- they all the three values are None.
- extract, project : functions
- functions to extract the full wave function from the projected wave
- functions, and project it back. Both are equal to None if the hopping
- is invertible.
- """
-
- eps = np.finfo(np.common_type(h_onslice, h_hop)).eps
-
- linsys = setup_linsys(h_onslice, h_hop)
-
- if isinstance(linsys, tuple):
- # In the singular case, it depends on the details of the system
- # whether one needs to solve a regular or a generalized
- # eigenproblem.
-
- assert len(linsys) == 3 or len(linsys) == 4, \
- "Corrupt lead eigenproblem data."
-
- if len(linsys) == 3:
- t, z, ev = kla.schur(linsys[0])
-
- # Right-decaying modes.
- select = abs(ev) > 1 + eps * tol
- # Propagating modes.
- propselect = abs(abs(ev) - 1) < eps * tol
-
- u, w, v = linsys[1]
- extract, project = linsys[2]
-
- 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:
- s, t, z, alpha, beta = kla.gen_schur(linsys[0], linsys[1],
- calc_q=False)
-
- # Right-decaying modes.
- select = abs(alpha) > (1 + eps * tol) * abs(beta)
- # Propagating modes.
- propselect = (abs(abs(alpha) - abs(beta)) <
- eps * tol * abs(beta))
-
- warning_settings = np.seterr(divide='ignore', invalid='ignore')
- ev = alpha / beta
- np.seterr(**warning_settings)
- # Note: the division is OK here, as we later only access
- # eigenvalues close to the unit circle
-
- u, w, v = linsys[2]
- extract, project = linsys[3]
-
- 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]
- else:
- # Hopping matrix can be safely inverted -> regular eigenproblem can be
- # used. This also means, that the hopping matrix is n x n square.
-
- t, z, ev = kla.schur(linsys)
-
- # Right-decaying modes.
- select = abs(ev) > 1 + eps * tol
- # Propagating modes.
- propselect = abs(abs(ev) - 1) < eps * tol
-
- # Signal that we are in the regular case.
- u = v = w = None
- extract = project = None
-
- 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]
-
- return ev, select, propselect, vec_gen, ord_schur,\
- u, w, v, extract, project
-
-
-def self_energy(h_onslice, h_hop, tol=1e6):
- """
- Compute the self-energy generated by a lead.
-
- The lead is described by the unit-cell
- Hamiltonian h_onslice and the hopping matrix h_hop.
-
- Parameters
- ----------
- h_onslice : NumPy array, real or complex, shape (N,N) The unit cell
- Hamiltonian of the lead slice.
- h_hop : NumPy array, real or complex, shape (N,M)
- the hopping matrix from a lead slice to the one on which self-energy
- has to be calculated (and any other hopping in the same direction).
-
- Returns
- -------
- Sigma : NumPy array, real or complex, shape (M,M)
- The computed self-energy. Note that even if `h_onslice` and `h_hop`
- are both real, `Sigma` will typically be complex. (More precisely, if
- there is a propagating mode, `Sigma` will definitely be complex.)
-
- Notes
- -----
- This function uses the most stable and efficient algorithm for calculating
- self-energy, described in kwant/doc/other/lead_modes.pdf
- """
-
- m = h_hop.shape[1]
-
- # Note: np.any(h_hop) returns (at least from NumPy 1.6.1 - 1.8-devel)
- # False if h_hop is purely imaginary
- if not (np.any(h_hop.real) or np.any(h_hop.imag)):
- return np.zeros((m, m))
-
- if (h_onslice.shape[0] != h_onslice.shape[1] or
- h_onslice.shape[0] != h_hop.shape[0]):
- raise ValueError("Incompatible matrix sizes for h_onslice and h_hop.")
-
- #defer most of the calculation to a helper routine (also used by modes)
- ev, select, propselect, vec_gen, ord_schur,\
- u, w, v, extract, project = unified_eigenproblem(h_onslice, h_hop, tol)
-
- if w is not None:
- n = w.size
- h_hop = w
- else:
- n = h_onslice.shape[0]
-
- # Compute the propagating eigenvectors, if they are present.
- nprop = np.sum(propselect)
-
- if nprop > 0:
- prop_vecs = vec_gen(propselect)
-
- prop_vecs, vel, rselect, crossing = \
- make_proper_modes(ev[propselect], prop_vecs, h_hop,
- extract, project)
- else:
- # Without propagating modes, the Schur methods certainly work.
- crossing = False
-
- if crossing:
- # Schur decomposition method does not work in this case, we need to
- # compute all the eigenvectors.
-
- # We already have the propagating ones, now we just need the
- # evanescent ones in addition.
-
- if nprop > 0:
- vecs = np.empty((2*n, n),
- dtype=np.common_type(ev, prop_vecs))
- else:
- vecs = np.empty((2*n, n), dtype=ev.dtype)
- # Note: rationale for the dtype: only if all the eigenvalues are real,
- # (which can only happen if the original eigenproblem was
- # real) and all propagating modes are real, the matrix of
- # eigenvectors will be real, too.
-
- if nprop > 0:
- nmodes = np.sum(rselect)
- vecs[:, :nmodes] = prop_vecs[:, rselect]
- else:
- nmodes = 0
-
- vecs[:, nmodes:] = vec_gen(select)
-
- if v is not None:
- return dot(v * w, dot(vecs[n:], dot(npl.inv(vecs[:n]),
- v.T.conj())))
- else:
- return dot(h_hop.T.conj(), dot(vecs[n:], npl.inv(vecs[:n])))
- else:
- # Reorder all the right-going eigenmodes to the top left part of
- # the Schur decomposition.
-
- if nprop > 0:
- select[propselect] = rselect
-
- z = ord_schur(select)
-
- if v is not None:
- return dot(v * w, dot(z[n:, :n], dot(npl.inv(z[:n, :n]),
- v.T.conj())))
- else:
- return dot(h_hop.T.conj(), dot(z[n:, :n], npl.inv(z[:n, :n])))
-
-Modes = namedtuple('Modes', ['vecs', 'vecslmbdainv', 'nmodes', 'svd'])
-
-
-def modes(h_onslice, h_hop, tol=1e6):
- """
- Compute the eigendecomposition of a translation operator of a lead.
-
- Parameters
- ----------
- h_onslice : NumPy array, real or complex, shape (N,N) The unit cell
- Hamiltonian of the lead slice.
- h_hop : NumPy array, real or complex, shape (N,M)
- the hopping matrix from a lead slice to the one on which self-energy
- has to be calculated (and any other hopping in the same direction).
-
- Returns
- -------
- Modes(vecs, vecslmbdainv, nmodes, svd) : a named tuple
- `vecs` is the matrix of eigenvectors of the translation operator.
- `vecslmbdainv` is the matrix of eigenvectors multiplied with their
- corresponding inverse eigenvalue. `nmodes` is the number of
- propagating modes in either direction. `svd` is a tuple (u, s, v)
- holding the singular value decomposition of the hopping matrix, or a
- single None if `h_hop` is invertible.
-
- Notes
- -----
- Only propagating modes and modes decaying away from the system are
- returned: The first `nmodes` columns in `vecs` correspond to incoming modes
- (coming from the lead into the system), the following `nmodes` columns
- correspond to outgoing modes (going into the lead, away from the system),
- the remaining columns are evanescent modes, decaying away from the system.
-
- 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.
-
- If `h_hop` is invertible, the full transverse wave functions are returned.
- If it is singular, the projections (u^dagger psi, v^dagger psi lambda^-1)
- are returned.
-
- In order for the linear system to be well-defined, instead of the
- evanescent modes, an orthogonal basis in the space of evanescent modes is
- returned.
-
- Propagating modes with the same lambda are orthogonalized. All the
- propagating modes are normalized by current.
-
- This function uses the most stable and efficient algorithm for calculating
- self-energy, described in kwant/doc/other/lead_modes.pdf
- """
-
- m = h_hop.shape[1]
-
- if (h_onslice.shape[0] != h_onslice.shape[1] or
- h_onslice.shape[0] != h_hop.shape[0]):
- raise ValueError("Incompatible matrix sizes for h_onslice and h_hop.")
-
- # Note: np.any(h_hop) returns (at least from NumPy 1.6.1 - 1.8-devel)
- # False if h_hop is purely imaginary
- if not (np.any(h_hop.real) or np.any(h_hop.imag)):
- n = h_hop.shape[0]
- svd = (np.empty((n, 0)), np.empty((0, 0)), np.empty((0, m)))
- return Modes(np.empty((0, 0)), np.empty((0, 0)), 0, svd)
-
- # Defer most of the calculation to a helper routine.
- ev, evanselect, propselect, vec_gen, ord_schur, \
- u, s, v, extract, project = unified_eigenproblem(h_onslice, h_hop, tol)
-
- if s is not None:
- n = s.size
- else:
- n = h_onslice.shape[0]
-
- nprop = np.sum(propselect)
- nevan = n - nprop // 2
- evanselect_bool = np.zeros((2*n), dtype='bool')
- evanselect_bool[evanselect] = True
- evan_vecs = ord_schur(evanselect)[:, :nevan]
-
- if nprop > 0:
- # Compute the propagating eigenvectors.
- prop_vecs = vec_gen(propselect)
-
- # Compute their velocity, and, if necessary, rotate them
- if s is not None:
- prop_vecs, vel, rprop, crossing = \
- make_proper_modes(ev[propselect], prop_vecs, s,
- extract, project)
- else:
- prop_vecs, vel, rprop, crossing = \
- make_proper_modes(ev[propselect], prop_vecs, h_hop)
-
- # Normalize propagating eigenvectors by velocities.
- prop_vecs /= np.sqrt(abs(vel))
-
- # Fix phase factor - make maximum of transverse wave function real
- # TODO (Anton): Take care of multiple maxima when normalizing.
- maxnode = prop_vecs[n + np.argmax(abs(prop_vecs[n:, :]), axis=0),
- np.arange(prop_vecs.shape[1])]
- maxnode /= abs(maxnode)
- prop_vecs /= maxnode
-
- lprop = np.logical_not(rprop)
- nmodes = np.sum(rprop)
-
- # Sort modes according to their k-vector (1/lambda = e^{-ik}):
- # - right-going modes: sort that k is in ascending order
- # - left-going modes: sort that k is in descending order
- # (note that k can be positive or negative). With this convention,
- # the modes of a simple square lattice strip are ordered as
- # expected (lowest subband first, etc.)
-
- prop_ev = ev[propselect]
- rsort = np.argsort((-1j * np.log(prop_ev[rprop])).real)
- lsort = np.argsort((1j * np.log(prop_ev[lprop])).real)
-
- vecs = np.c_[prop_vecs[n:, lprop][:, lsort],
- prop_vecs[n:, rprop][:, rsort],
- evan_vecs[n:]]
- vecslmbdainv = np.c_[prop_vecs[:n, lprop][:, lsort],
- prop_vecs[:n, rprop][:, rsort],
- evan_vecs[:n]]
-
- else:
- vecs = evan_vecs[n:]
- vecslmbdainv = evan_vecs[:n]
- nmodes = 0
-
- svd = None if s is None else (u, s, v)
- return Modes(vecs, vecslmbdainv, nmodes, svd)
-
-
-def square_self_energy(width, hopping, potential, fermi_energy):
- """
- 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
-
- # p labels transversal modes. i and j label the sites of a slice.
-
- # 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:
- return q/2 - 1j * sqrt(1 - (q / 2) ** 2)
- else:
- return q/2 - copysign(sqrt((q / 2) ** 2 - 1), q)
- f_p = np.empty((width,), dtype=complex)
- for p in xrange(width):
- e = 2 * hopping * (1 - cos(factor * (p + 1)))
- q = (fermi_energy - potential - e) / hopping - 2
- 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
diff --git a/kwant/physics/tests/test_leads.py b/kwant/physics/tests/test_leads.py
new file mode 100644
index 0000000000000000000000000000000000000000..88e0ea809f90ad9ecb2889ac4d449297d8a4d1b7
--- /dev/null
+++ b/kwant/physics/tests/test_leads.py
@@ -0,0 +1,280 @@
+# Copyright 2011-2013 kwant authors.
+#
+# This file is part of kwant. It is subject to the license terms in the
+# LICENSE file found in the top-level directory of this distribution and at
+# http://kwant-project.org/license. A list of kwant authors can be found in
+# the AUTHORS file at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from __future__ import division
+import numpy as np
+from itertools import product, izip
+from numpy.testing import assert_almost_equal
+from kwant.physics import leads
+import kwant
+
+modes_se = leads.selfenergy
+
+def h_slice(t, w, e):
+ h = (4 * t - e) * np.identity(w)
+ h.flat[1 :: w + 1] = -t
+ h.flat[w :: w + 1] = -t
+ return h
+
+
+def test_analytic_numeric():
+ w = 5 # width
+ t = 0.78 # hopping element
+ e = 1.3 # Fermi energy
+
+ assert_almost_equal(leads.square_selfenergy(w, t, e),
+ modes_se(h_slice(t, w, e), -t * np.identity(w)))
+
+
+def test_regular_fully_degenerate():
+ """Selfenergy with an invertible hopping matrix, and degenerate bands."""
+
+ w = 6 # width
+ t = 0.5 # hopping element
+ e = 1.3 # Fermi energy
+
+ h_hop_s = -t * np.identity(w)
+ h_onslice_s = h_slice(t, w, e)
+
+ h_hop = np.zeros((2*w, 2*w))
+ h_hop[:w, :w] = h_hop_s
+ h_hop[w:, w:] = h_hop_s
+
+ h_onslice = np.zeros((2*w, 2*w))
+ h_onslice[:w, :w] = h_onslice_s
+ h_onslice[w:, w:] = h_onslice_s
+
+ g = np.zeros((2*w, 2*w), dtype=complex)
+ g[:w, :w] = leads.square_selfenergy(w, t, e)
+ g[w:, w:] = leads.square_selfenergy(w, t, e)
+
+ assert_almost_equal(g, modes_se(h_onslice, h_hop))
+
+
+def test_regular_degenerate_with_crossing():
+ """This is a testcase with invertible hopping matrices,
+ and degenerate k-values with a crossing such that one
+ mode has a positive velocity, and one a negative velocity
+
+ For this case the fall-back technique must be used.
+ """
+
+ w = 4 # width
+ t = 0.5 # hopping element
+ e = 1.8 # Fermi energy
+
+ global h_hop
+ h_hop_s = -t * np.identity(w)
+ h_onslice_s = h_slice(t, w, e)
+
+ hop = np.zeros((2*w, 2*w))
+ hop[:w, :w] = h_hop_s
+ hop[w:, w:] = -h_hop_s
+
+ h_onslice = np.zeros((2*w, 2*w))
+ h_onslice[:w, :w] = h_onslice_s
+ h_onslice[w:, w:] = -h_onslice_s
+
+ g = np.zeros((2*w, 2*w), dtype=complex)
+ g[:w, :w] = leads.square_selfenergy(w, t, e)
+ g[w:, w:] = -np.conj(leads.square_selfenergy(w, t, e))
+
+ assert_almost_equal(g, modes_se(h_onslice, hop))
+
+
+def test_singular():
+ """This testcase features a rectangular (and hence singular)
+ hopping matrix without degeneracies.
+
+ This case can be treated with the Schur technique."""
+
+ w = 5 # width
+ t = .5 # hopping element
+ e = 0.4 # Fermi energy
+
+ h_hop_s = -t * np.identity(w)
+ h_onslice_s = h_slice(t, w, e)
+
+ h_hop = np.zeros((2*w, w))
+ h_hop[w:, :w] = h_hop_s
+
+ h_onslice = np.zeros((2*w, 2*w))
+ h_onslice[:w, :w] = h_onslice_s
+ h_onslice[:w, w:] = h_hop_s
+ h_onslice[w:, :w] = h_hop_s
+ h_onslice[w:, w:] = h_onslice_s
+ g = leads.square_selfenergy(w, t, e)
+
+ assert_almost_equal(g, modes_se(h_onslice, h_hop))
+
+
+def test_singular_but_square():
+ """This testcase features a singular, square hopping matrices
+ without degeneracies.
+
+ This case can be treated with the Schur technique."""
+
+ w = 5 # width
+ t = 0.9 # hopping element
+ e = 2.38 # Fermi energy
+
+ h_hop_s = -t * np.identity(w)
+ h_onslice_s = h_slice(t, w, e)
+
+ h_hop = np.zeros((2*w, 2*w))
+ h_hop[w:, :w] = h_hop_s
+
+ h_onslice = np.zeros((2*w, 2*w))
+ h_onslice[:w, :w] = h_onslice_s
+ h_onslice[:w, w:] = h_hop_s
+ h_onslice[w:, :w] = h_hop_s
+ h_onslice[w:, w:] = h_onslice_s
+
+ g = np.zeros((2*w, 2*w), dtype=complex)
+ g[:w, :w] = leads.square_selfenergy(w, t, e)
+ assert_almost_equal(g, modes_se(h_onslice, h_hop))
+
+
+def test_singular_fully_degenerate():
+ """This testcase features a rectangular (and hence singular)
+ hopping matrix with complete degeneracy.
+
+ This case can still be treated with the Schur technique."""
+
+ w = 5 # width
+ t = 1.5 # hopping element
+ e = 3.3 # Fermi energy
+
+ h_hop_s = -t * np.identity(w)
+ h_onslice_s = h_slice(t, w, e)
+
+ h_hop = np.zeros((4*w, 2*w))
+ h_hop[2*w:3*w, :w] = h_hop_s
+ h_hop[3*w:4*w, w:2*w] = h_hop_s
+
+ h_onslice = np.zeros((4*w, 4*w))
+ h_onslice[:w, :w] = h_onslice_s
+ h_onslice[:w, 2*w:3*w] = h_hop_s
+ h_onslice[w:2*w, w:2*w] = h_onslice_s
+ h_onslice[w:2*w, 3*w:4*w] = h_hop_s
+ h_onslice[2*w:3*w, :w] = h_hop_s
+ h_onslice[2*w:3*w, 2*w:3*w] = h_onslice_s
+ h_onslice[3*w:4*w, w:2*w] = h_hop_s
+ h_onslice[3*w:4*w, 3*w:4*w] = h_onslice_s
+
+ g = np.zeros((2*w, 2*w), dtype=complex)
+ g[:w, :w] = leads.square_selfenergy(w, t, e)
+ g[w:, w:] = leads.square_selfenergy(w, t, e)
+
+ assert_almost_equal(g, modes_se(h_onslice, h_hop))
+
+
+def test_singular_degenerate_with_crossing():
+ """This testcase features a rectangular (and hence singular)
+ hopping matrix with degeneracy k-values including a crossing
+ with velocities of opposite sign.
+
+ This case must be treated with the fall-back technique."""
+
+ w = 5 # width
+ t = 20.5 # hopping element
+ e = 3.3 # Fermi energy
+
+ h_hop_s = -t * np.identity(w)
+ h_onslice_s = h_slice(t, w, e)
+
+ h_hop = np.zeros((4*w, 2*w))
+ h_hop[2*w:3*w, :w] = h_hop_s
+ h_hop[3*w:4*w, w:2*w] = -h_hop_s
+
+ h_onslice = np.zeros((4*w, 4*w))
+ h_onslice[:w, :w] = h_onslice_s
+ h_onslice[:w, 2*w:3*w] = h_hop_s
+ h_onslice[w:2*w, w:2*w] = -h_onslice_s
+ h_onslice[w:2*w, 3*w:4*w] = -h_hop_s
+ h_onslice[2*w:3*w, :w] = h_hop_s
+ h_onslice[2*w:3*w, 2*w:3*w] = h_onslice_s
+ h_onslice[3*w:4*w, w:2*w] = -h_hop_s
+ h_onslice[3*w:4*w, 3*w:4*w] = -h_onslice_s
+
+ g = np.zeros((2*w, 2*w), dtype=complex)
+ g[:w, :w] = leads.square_selfenergy(w, t, e)
+ g[w:, w:] = -np.conj(leads.square_selfenergy(w, t, e))
+
+ assert_almost_equal(g, modes_se(h_onslice, h_hop))
+
+
+def test_singular_h_and_t():
+ h = 0.1 * np.identity(6)
+ t = np.eye(6, 6, 4)
+ sigma = modes_se(h, t)
+ sigma_should_be = np.zeros((6,6))
+ sigma_should_be[4, 4] = sigma_should_be[5, 5] = -10
+ assert_almost_equal(sigma, sigma_should_be)
+
+
+def test_modes():
+ h, t = .3, .7
+ vecs, vecslinv, nrprop, svd = leads.modes(np.array([[h]]), np.array([[t]]))
+ assert nrprop == 1
+ assert svd is None
+ np.testing.assert_almost_equal((vecs[0] * vecslinv[0].conj()).imag,
+ [0.5, -0.5])
+
+
+def test_modes_bearded_ribbon():
+ # Check if bearded graphene ribbons work.
+ lat = kwant.lattice.honeycomb()
+ sys = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
+ sys[lat.shape((lambda pos: -20 < pos[1] < 20),
+ (0, 0))] = 0.3
+ sys[lat.neighbors()] = -1
+ sys = sys.finalized()
+ h, t = sys.slice_hamiltonian(), sys.inter_slice_hopping()
+ assert leads.modes(h, t).nmodes == 8 # Calculated by plotting dispersion.
+
+
+def test_algorithm_equivalence():
+ np.random.seed(400)
+ n = 12
+ h = np.random.randn(n, n) + 1j * np.random.randn(n, n)
+ h += h.T.conj()
+ t = np.random.randn(n, n) + 1j * np.random.randn(n, n)
+ u, s, vh = np.linalg.svd(t)
+ u, v = u * np.sqrt(s), vh.T.conj() * np.sqrt(s)
+ prop_vecs = []
+ evan_vecs = []
+ algos = [(True,)] + list(product(*([(False,)] + 2 * [(True, False)])))
+ for algo in algos:
+ result = leads.modes(h, t, algorithm=algo)
+
+ vecs, vecslmbdainv = result.vecs, result.vecslmbdainv
+
+ # Bring the calculated vectors to real space
+ if not algo[0]:
+ vecs = np.dot(v, vecs)
+ np.testing.assert_almost_equal(result.sqrt_hop, v)
+ else:
+ vecslmbdainv = np.dot(v.T.conj(), vecslmbdainv)
+ full_vecs = np.r_[vecslmbdainv, vecs]
+
+ prop_vecs.append(full_vecs[:, : 2 * result.nmodes])
+ evan_vecs.append(full_vecs[:, 2 * result.nmodes :])
+
+ msg = 'Algorithm {0} failed.'
+ for vecs, algo in izip(prop_vecs, algos):
+ # Propagating modes should have identical ordering, and only vary
+ # By a phase
+ np.testing.assert_allclose(np.abs(np.sum(vecs/prop_vecs[0],
+ axis=0)), vecs.shape[0],
+ err_msg=msg.format(algo))
+
+ for vecs, algo in izip(evan_vecs, algos):
+ # Evanescent modes must span the same linear space.
+ assert np.linalg.matrix_rank(np.c_[vecs, evan_vecs[0]], tol=1e-12) == \
+ vecs.shape[1], msg.format(algo)
diff --git a/kwant/physics/tests/test_noise.py b/kwant/physics/tests/test_noise.py
index bf5f1d58f5a67ee9d58e1d25fce1b8a68c954d4f..4a9c25bf24c439801f74a8740ab315cfa4f71956 100644
--- a/kwant/physics/tests/test_noise.py
+++ b/kwant/physics/tests/test_noise.py
@@ -43,9 +43,9 @@ class LeadWithOnlySelfEnergy(object):
def __init__(self, lead):
self.lead = lead
- def self_energy(self, energy, args=()):
+ def selfenergy(self, energy, args=()):
assert args == ()
- return self.lead.self_energy(energy)
+ return self.lead.selfenergy(energy)
def test_twoterminal():
diff --git a/kwant/physics/tests/test_selfenergy.py b/kwant/physics/tests/test_selfenergy.py
deleted file mode 100644
index 2975c1415f3abd1c571a48e43a98e74c0d0b7677..0000000000000000000000000000000000000000
--- a/kwant/physics/tests/test_selfenergy.py
+++ /dev/null
@@ -1,250 +0,0 @@
-# Copyright 2011-2013 kwant authors.
-#
-# This file is part of kwant. It is subject to the license terms in the
-# LICENSE file found in the top-level directory of this distribution and at
-# http://kwant-project.org/license. A list of kwant authors can be found in
-# the AUTHORS file at the top-level directory of this distribution and at
-# http://kwant-project.org/authors.
-
-from __future__ import division
-import numpy as np
-from numpy.testing import assert_almost_equal
-import kwant.physics.selfenergy as se
-
-def test_analytic_numeric():
- w = 5 # width
- t = 0.5 # hopping element
- v = 2 # potential
- e = 3.3 # Fermi energy
-
- h_hop = -t * np.identity(w)
- h_onslice = ((v + 4 * t - e)
- * np.identity(w))
- h_onslice.flat[1 :: w + 1] = -t
- h_onslice.flat[w :: w + 1] = -t
-
- assert_almost_equal(se.square_self_energy(w, t, v, e),
- se.self_energy(h_onslice, h_hop))
-
-def test_regular_fully_degenerate():
- """This testcase features an invertible hopping matrix,
- and bands that are fully degenerate.
-
- This case can still be treated with the Schur technique."""
-
- w = 5 # width
- t = 0.5 # hopping element
- v = 2 # potential
- e = 3.3 # Fermi energy
-
- h_hop_s = -t * np.identity(w)
- h_onslice_s = ((v + 4 * t - e)
- * np.identity(w))
- h_onslice_s.flat[1 :: w + 1] = -t
- h_onslice_s.flat[w :: w + 1] = -t
-
- h_hop = np.zeros((2*w, 2*w))
- h_hop[0:w, 0:w] = h_hop_s
- h_hop[w:2*w, w:2*w] = h_hop_s
-
- h_onslice = np.zeros((2*w, 2*w))
- h_onslice[0:w, 0:w] = h_onslice_s
- h_onslice[w:2*w, w:2*w] = h_onslice_s
-
- g = np.zeros((2*w, 2*w), dtype=complex)
- g[0:w, 0:w] = se.square_self_energy(w, t, v, e)
- g[w:2*w, w:2*w] = se.square_self_energy(w, t, v, e)
-
- assert_almost_equal(g,
- se.self_energy(h_onslice, h_hop))
-
-def test_regular_degenerate_with_crossing():
- """This is a testcase with invertible hopping matrices,
- and degenerate k-values with a crossing such that one
- mode has a positive velocity, and one a negative velocity
-
- For this case the fall-back technique must be used.
- """
-
- w = 5 # width
- t = 0.5 # hopping element
- v = 2 # potential
- e = 3.3 # Fermi energy
-
- h_hop_s = -t * np.identity(w)
- h_onslice_s = ((v + 4 * t - e)
- * np.identity(w))
- h_onslice_s.flat[1 :: w + 1] = -t
- h_onslice_s.flat[w :: w + 1] = -t
-
- h_hop = np.zeros((2*w, 2*w))
- h_hop[0:w, 0:w] = h_hop_s
- h_hop[w:2*w, w:2*w] = -h_hop_s
-
- h_onslice = np.zeros((2*w, 2*w))
- h_onslice[0:w, 0:w] = h_onslice_s
- h_onslice[w:2*w, w:2*w] = -h_onslice_s
-
- g = np.zeros((2*w, 2*w), dtype=complex)
- g[0:w, 0:w] = se.square_self_energy(w, t, v, e)
- g[w:2*w, w:2*w] = -np.conj(se.square_self_energy(w, t, v, e))
-
- assert_almost_equal(g,
- se.self_energy(h_onslice, h_hop))
-
-def test_singular():
- """This testcase features a rectangular (and hence singular)
- hopping matrix without degeneracies.
-
- This case can be treated with the Schur technique."""
-
- w = 5 # width
- t = 0.5 # hopping element
- v = 2 # potential
- e = 3.3 # Fermi energy
-
- h_hop_s = -t * np.identity(w)
- h_onslice_s = ((v + 4 * t - e)
- * np.identity(w))
- h_onslice_s.flat[1 :: w + 1] = -t
- h_onslice_s.flat[w :: w + 1] = -t
-
- h_hop = np.zeros((2*w, w))
- h_hop[w:2*w, 0:w] = h_hop_s
-
- h_onslice = np.zeros((2*w, 2*w))
- h_onslice[0:w, 0:w] = h_onslice_s
- h_onslice[0:w, w:2*w] = h_hop_s
- h_onslice[w:2*w, 0:w] = h_hop_s
- h_onslice[w:2*w, w:2*w] = h_onslice_s
-
- assert_almost_equal(se.square_self_energy(w, t, v, e),
- se.self_energy(h_onslice, h_hop))
-
-def test_singular_but_square():
- """This testcase features a singular, square hopping matrices
- without degeneracies.
-
- This case can be treated with the Schur technique."""
-
- w = 5 # width
- t = 0.5 # hopping element
- v = 2 # potential
- e = 3.3 # Fermi energy
-
- h_hop_s = -t * np.identity(w)
- h_onslice_s = ((v + 4 * t - e)
- * np.identity(w))
- h_onslice_s.flat[1 :: w + 1] = -t
- h_onslice_s.flat[w :: w + 1] = -t
-
- h_hop = np.zeros((2*w, 2*w))
- h_hop[w:2*w, 0:w] = h_hop_s
-
- h_onslice = np.zeros((2*w, 2*w))
- h_onslice[0:w, 0:w] = h_onslice_s
- h_onslice[0:w, w:2*w] = h_hop_s
- h_onslice[w:2*w, 0:w] = h_hop_s
- h_onslice[w:2*w, w:2*w] = h_onslice_s
-
- g = np.zeros((2*w, 2*w), dtype=complex)
- g[0:w, 0:w] = se.square_self_energy(w, t, v, e)
-
- assert_almost_equal(g,
- se.self_energy(h_onslice, h_hop))
-
-def test_singular_fully_degenerate():
- """This testcase features a rectangular (and hence singular)
- hopping matrix with complete degeneracy.
-
- This case can still be treated with the Schur technique."""
-
- w = 5 # width
- t = 0.5 # hopping element
- v = 2 # potential
- e = 3.3 # Fermi energy
-
- h_hop_s = -t * np.identity(w)
- h_onslice_s = ((v + 4 * t - e)
- * np.identity(w))
- h_onslice_s.flat[1 :: w + 1] = -t
- h_onslice_s.flat[w :: w + 1] = -t
-
- h_hop = np.zeros((4*w, 2*w))
- h_hop[2*w:3*w, 0:w] = h_hop_s
- h_hop[3*w:4*w, w:2*w] = h_hop_s
-
- h_onslice = np.zeros((4*w, 4*w))
- h_onslice[0:w, 0:w] = h_onslice_s
- h_onslice[0:w, 2*w:3*w] = h_hop_s
- h_onslice[w:2*w, w:2*w] = h_onslice_s
- h_onslice[w:2*w, 3*w:4*w] = h_hop_s
- h_onslice[2*w:3*w, 0:w] = h_hop_s
- h_onslice[2*w:3*w, 2*w:3*w] = h_onslice_s
- h_onslice[3*w:4*w, w:2*w] = h_hop_s
- h_onslice[3*w:4*w, 3*w:4*w] = h_onslice_s
-
- g = np.zeros((2*w, 2*w), dtype=complex)
- g[0:w, 0:w] = se.square_self_energy(w, t, v, e)
- g[w:2*w, w:2*w] = se.square_self_energy(w, t, v, e)
-
- assert_almost_equal(g,
- se.self_energy(h_onslice, h_hop))
-
-def test_singular_degenerate_with_crossing():
- """This testcase features a rectangular (and hence singular)
- hopping matrix with degeneracy k-values including a crossing
- with velocities of opposite sign.
-
- This case must be treated with the fall-back technique."""
-
- w = 5 # width
- t = 0.5 # hopping element
- v = 2 # potential
- e = 3.3 # Fermi energy
-
- h_hop_s = -t * np.identity(w)
- h_onslice_s = ((v + 4 * t - e)
- * np.identity(w))
- h_onslice_s.flat[1 :: w + 1] = -t
- h_onslice_s.flat[w :: w + 1] = -t
-
- h_hop = np.zeros((4*w, 2*w))
- h_hop[2*w:3*w, 0:w] = h_hop_s
- h_hop[3*w:4*w, w:2*w] = -h_hop_s
-
- h_onslice = np.zeros((4*w, 4*w))
- h_onslice[0:w, 0:w] = h_onslice_s
- h_onslice[0:w, 2*w:3*w] = h_hop_s
- h_onslice[w:2*w, w:2*w] = -h_onslice_s
- h_onslice[w:2*w, 3*w:4*w] = -h_hop_s
- h_onslice[2*w:3*w, 0:w] = h_hop_s
- h_onslice[2*w:3*w, 2*w:3*w] = h_onslice_s
- h_onslice[3*w:4*w, w:2*w] = -h_hop_s
- h_onslice[3*w:4*w, 3*w:4*w] = -h_onslice_s
-
- g = np.zeros((2*w, 2*w), dtype=complex)
- g[0:w, 0:w] = se.square_self_energy(w, t, v, e)
- g[w:2*w, w:2*w] = -np.conj(se.square_self_energy(w, t, v, e))
-
- assert_almost_equal(g,
- se.self_energy(h_onslice, h_hop))
-
-def test_singular_h_and_t():
- h = 0.1 * np.identity(6)
- t = np.eye(6, 6, 4)
- sigma = se.self_energy(h, t)
- sigma_should_be = np.zeros((6,6))
- sigma_should_be[4, 4] = sigma_should_be[5, 5] = -10
- assert_almost_equal(se.self_energy(h, t), sigma_should_be)
-
-def test_modes():
- h, t = .3, .7
- vecs, vecslinv, nrpop, svd = se.modes(np.array([[h]]), np.array([[t]]))
- l = (np.sqrt(h**2 - 4 * t**2 + 0j) - h) / (2 * t)
- current = np.sqrt(4 * t**2 - h**2)
- assert nrpop == 1
- assert svd is None
- np.testing.assert_almost_equal(vecs, [2 * [1/np.sqrt(current)]])
- np.testing.assert_almost_equal(vecslinv,
- vecs * np.array([1/l, 1/l.conj()]))
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index b8d25e9dbb86ddac3b75ec19ef2d70892029cfc7..2f3b654ddd028416027f34de2c0867f499a7da1a 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -94,8 +94,7 @@ class SparseSolver(object):
def _make_linear_sys(self, sys, out_leads, in_leads, energy=0,
force_realspace=False, check_hermiticity=True,
args=()):
- """
- Make a sparse linear system of equations defining a scattering
+ """Make a sparse linear system of equations defining a scattering
problem.
Parameters
@@ -129,9 +128,9 @@ class SparseSolver(object):
a small part of the complete solution).
lead_info : list of objects
Contains one entry for each lead. For a lead defined as a
- tight-binding system, this is an instance of `kwant.physics.Modes`
- (as returned by `kwant.physics.modes`), otherwise the lead
- self-energy matrix.
+ tight-binding system, this is an instance of
+ `~kwant.physics.ModesTuple` with a corresponding format,
+ otherwise the lead self-energy matrix.
Notes
-----
@@ -139,6 +138,7 @@ class SparseSolver(object):
or a low-level translationally invariant system.
The system of equations that is created is described in
kwant/doc/other/linear_system.pdf
+
"""
splhsmat = getattr(sp, self.lhsformat + '_matrix')
@@ -168,27 +168,18 @@ class SparseSolver(object):
lead_info = []
for leadnum, interface in enumerate(sys.lead_interfaces):
lead = sys.leads[leadnum]
- if isinstance(lead, system.InfiniteSystem) and not force_realspace:
- h = lead.slice_hamiltonian(args=args)
-
- if check_hermiticity:
- if not np.allclose(h, h.T.conj(), rtol=1e-13):
- msg = "Lead number {0} has a non-Hermitian " \
- "slice Hamiltonian."
- raise ValueError(msg.format(leadnum))
-
- h -= energy * np.identity(h.shape[0])
- v = lead.inter_slice_hopping(args=args)
- modes = physics.modes(h, v)
+ if hasattr(lead, 'modes') and not force_realspace:
+ modes = lead.modes(energy, args=args)
lead_info.append(modes)
+ u, ulinv, nprop, svd_v = modes
- # Note: np.any(v) returns (at least from NumPy 1.6.1 -
- # 1.8-devel) False if v is purely imaginary
- if not (np.any(v.real) or np.any(v.imag)):
+ if len(u) == 0:
# See comment about zero-shaped sparse matrices at the top.
rhs.append(np.zeros((lhs.shape[1], 0)))
continue
- u, ulinv, nprop, svd = modes
+
+ if svd_v is not None:
+ svd_v = svd_v.T.conj()
if leadnum in out_leads:
kept_vars.append(
@@ -206,27 +197,25 @@ class SparseSolver(object):
transf = sp.csc_matrix((np.ones(iface_orbs.size), coords),
shape=(iface_orbs.size, lhs.shape[0]))
- if svd is not None:
- v_sp = sp.csc_matrix(svd[2].T.conj()) * transf
+ if svd_v is not None:
+ v_sp = sp.csc_matrix(svd_v) * transf
vdaguout_sp = transf.T * \
- sp.csc_matrix(np.dot(svd[2] * svd[1], u_out))
+ sp.csc_matrix(np.dot(svd_v.T.conj(), u_out))
lead_mat = - ulinv_out
else:
v_sp = transf
- vdaguout_sp = transf.T * sp.csc_matrix(np.dot(v.T.conj(),
- u_out))
+ vdaguout_sp = transf.T * sp.csc_matrix(u_out)
lead_mat = - ulinv_out
lhs = sp.bmat([[lhs, vdaguout_sp], [v_sp, lead_mat]],
format=self.lhsformat)
if leadnum in in_leads and nprop > 0:
- if svd:
+ if svd_v is not None:
vdaguin_sp = transf.T * sp.csc_matrix(
- -np.dot(svd[2] * svd[1], u_in))
+ -np.dot(svd_v.T.conj(), u_in))
else:
- vdaguin_sp = transf.T * sp.csc_matrix(
- -np.dot(v.T.conj(), u_in))
+ vdaguin_sp = transf.T * sp.csc_matrix(-u_in)
# defer formation of the real matrix until the proper
# system size is known
@@ -235,7 +224,7 @@ class SparseSolver(object):
# See comment about zero-shaped sparse matrices at the top.
rhs.append(np.zeros((lhs.shape[1], 0)))
else:
- sigma = lead.self_energy(energy, args)
+ sigma = lead.selfenergy(energy, args)
lead_info.append(sigma)
indices = np.r_[tuple(slice(offsets[i], offsets[i + 1])
for i in interface)]
@@ -404,7 +393,7 @@ class SparseSolver(object):
self._make_linear_sys(fsys, [], xrange(len(fsys.leads)), energy,
args=args)
- Modes = physics.Modes
+ Modes = physics.ModesTuple
num_extra_vars = sum(li.vecs.shape[1] - li.nmodes
for li in lead_info if isinstance(li, Modes))
num_orb = h.shape[0] - num_extra_vars
@@ -467,7 +456,7 @@ class WaveFunction(object):
(h, self.rhs, kept_vars), lead_info = \
solver._make_linear_sys(sys, [], xrange(len(sys.leads)),
energy, args=args)
- Modes = physics.Modes
+ Modes = physics.ModesTuple
num_extra_vars = sum(li.vecs.shape[1] - li.nmodes
for li in lead_info if isinstance(li, Modes))
self.solver = solver
diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py
index 8945a9323ad839ca539f5cbf128f8991478387b0..d6f61f7f2ac27e38e9bc477057bf6cf589d35d90 100644
--- a/kwant/solvers/tests/_test_sparse.py
+++ b/kwant/solvers/tests/_test_sparse.py
@@ -21,9 +21,9 @@ class LeadWithOnlySelfEnergy(object):
def __init__(self, lead):
self.lead = lead
- def self_energy(self, energy, args=()):
+ def selfenergy(self, energy, args=()):
assert args == ()
- return self.lead.self_energy(energy)
+ return self.lead.selfenergy(energy)
# Test output sanity: that an error is raised if no output is requested,
@@ -55,7 +55,7 @@ def test_output(solve):
s2, modes2 = result2.data, result2.lead_info
assert s2.shape == (modes2[1][2], modes2[0][2])
assert_almost_equal(s1, s2)
- assert_almost_equal(np.dot(s.conjugate().transpose(), s),
+ assert_almost_equal(np.dot(s.T.conj(), s),
np.identity(s.shape[0]))
assert_raises(ValueError, solve, fsys, 0, [])
modes = solve(fsys).lead_info
@@ -272,7 +272,7 @@ def test_tricky_singular_hopping(solve):
# Test equivalence between self-energy and scattering matrix representations.
# Also check that transmission works.
-def test_self_energy(solve):
+def test_selfenergy(solve):
np.random.seed(4)
system = kwant.Builder()
left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
@@ -315,7 +315,7 @@ def test_self_energy(solve):
assert_almost_equal(t_should_be, sol.transmission(1, 0))
-def test_self_energy_reflection(solve):
+def test_selfenergy_reflection(solve):
np.random.seed(4)
system = kwant.Builder()
left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
diff --git a/kwant/solvers/tests/test_mumps.py b/kwant/solvers/tests/test_mumps.py
index 9692420c4c5a722d4e50b0c4cd16c1c62a8a910a..2b1b8b15d1467285a00d8aa443bc3fcf4fd2b827 100644
--- a/kwant/solvers/tests/test_mumps.py
+++ b/kwant/solvers/tests/test_mumps.py
@@ -6,7 +6,7 @@
# the AUTHORS file at the top-level directory of this distribution and at
# http://kwant-project.org/authors.
-from nose.plugins.skip import Skip, SkipTest
+#from nose.plugins.skip import Skip, SkipTest
from numpy.testing.decorators import skipif
try:
from kwant.solvers.mumps import \
@@ -81,19 +81,19 @@ def test_tricky_singular_hopping():
@skipif(_no_mumps)
-def test_self_energy():
+def test_selfenergy():
for opts in opt_list:
reset_options()
options(**opts)
- _test_sparse.test_self_energy(solve)
+ _test_sparse.test_selfenergy(solve)
@skipif(_no_mumps)
-def test_self_energy_reflection():
+def test_selfenergy_reflection():
for opts in opt_list:
reset_options()
options(**opts)
- _test_sparse.test_self_energy_reflection(solve)
+ _test_sparse.test_selfenergy_reflection(solve)
@skipif(_no_mumps)
diff --git a/kwant/solvers/tests/test_sparse.py b/kwant/solvers/tests/test_sparse.py
index 9f3edd694d9d8101f0810250906d4fb6ad86daa9..46b58a16587af765bd5bf84ecb836e81c36b2847 100644
--- a/kwant/solvers/tests/test_sparse.py
+++ b/kwant/solvers/tests/test_sparse.py
@@ -39,12 +39,12 @@ def test_tricky_singular_hopping():
_test_sparse.test_tricky_singular_hopping(solve)
-def test_self_energy():
- _test_sparse.test_self_energy(solve)
+def test_selfenergy():
+ _test_sparse.test_selfenergy(solve)
-def test_self_energy_reflection():
- _test_sparse.test_self_energy_reflection(solve)
+def test_selfenergy_reflection():
+ _test_sparse.test_selfenergy_reflection(solve)
def test_very_singular_leads():
diff --git a/kwant/system.py b/kwant/system.py
index 925117881531eb2cba14654965c57300b526e801..fe931001ba308cc676f225abb16713dc27b199ca 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -13,6 +13,7 @@ __all__ = ['System', 'FiniteSystem', 'InfiniteSystem']
import abc
import types
+from copy import copy
import numpy as np
from . import physics, _system
@@ -67,19 +68,22 @@ class FiniteSystem(System):
Instance Variables
------------------
- leads : sequence of lead objects
- Each lead object has to provide a method ``self_energy(energy, args)``.
+ leads : sequence of leads
+ Each lead has to provide a method
+ ``selfenergy(energy, args)``.
+ It may provide ``modes(energy, args)`` as well.
lead_interfaces : sequence of sequences of integers
- Each sub-sequence contains the indices of the system sites to which the
- lead is connected.
+ Each sub-sequence contains the indices of the system sites
+ to which the lead is connected.
Notes
-----
The length of `leads` must be equal to the length of `lead_interfaces`.
- For lead ``n``, the method leads[n].self_energy must return a square matrix
- whose size is ``sum(self.num_orbitals(neighbor) for neighbor in
- self.lead_interfaces[n])``.
+ For lead ``n``, the method leads[n].selfenergy must return a square matrix
+ whose size is ``sum(self.num_orbitals(neighbor)`` for neighbor in
+ self.lead_interfaces[n])``. The output format for ``leads[n].modes`` has to
+ be as described in `~kwant.physics.ModesTuple`.
Often, the elements of `leads` will be instances of `InfiniteSystem`. If
this is the case for lead ``n``, the sites ``lead_interfaces[n]`` match
@@ -87,6 +91,57 @@ class FiniteSystem(System):
"""
__metaclass__ = abc.ABCMeta
+ def precalculate(self, energy, args=(), leads=None,
+ calculate_selfenergy=True):
+ """
+ Precalculate modes or self-energies in the leads.
+
+ Construct a copy of the system, with the lead modes precalculated,
+ which may significantly speed up calculations where only the system
+ is changing.
+
+ Parameters
+ ----------
+ energy : float
+ Energy at which the modes or self-energies have to be
+ evaluated.
+ args : sequence
+ Additional parameters required for calculating the Hamiltionians
+ leads : list of integers or None
+ Numbers of the leads to be precalculated. If `None`, all are
+ precalculated.
+ calculate_selfenergy : bool
+ Whether to calculate self-energy if modes are available.
+
+ Returns
+ -------
+ sys : FiniteSystem
+ A copy of the original system with some leads precalculated.
+
+ Notes
+ -----
+ If the leads are precalculated at certain `energy` or `args` values,
+ they might give wrong results if used to solve the system with
+ different parameter values. Use this function with caution.
+ """
+ result = copy(self)
+ if leads is None:
+ leads = range(len(self.leads))
+ new_leads = []
+ for nr, lead in enumerate(self.leads):
+ if nr not in leads:
+ new_leads.append(lead)
+ continue
+ modes, selfenergy = None, None
+ try:
+ modes = lead.modes(energy, args)
+ if calculate_selfenergy:
+ selfenergy = physics.selfenergy_from_modes(modes)
+ except AttributeError:
+ selfenergy = lead.selfenergy(energy, args)
+ new_leads.append(PrecalculatedLead(modes, selfenergy))
+ result.leads = new_leads
+ return result
class InfiniteSystem(System):
"""
@@ -145,7 +200,21 @@ class InfiniteSystem(System):
return self.hamiltonian_submatrix(slice_sites, neighbor_sites,
sparse=sparse, args=args)
- def self_energy(self, energy, args=()):
+ def modes(self, energy, args=()):
+ """Return mode decomposition of the lead
+
+ See documentation of `~kwant.physics.ModesTuple` for the return
+ format details.
+ """
+ ham = self.slice_hamiltonian(args=args)
+ shape = ham.shape
+ assert len(shape) == 2
+ assert shape[0] == shape[1]
+ # Subtract energy from the diagonal.
+ ham.flat[::ham.shape[0] + 1] -= energy
+ return physics.modes(ham, self.inter_slice_hopping(args=args))
+
+ def selfenergy(self, energy, args=()):
"""Return self-energy of a lead.
The returned matrix has the shape (s, s), where s is
@@ -158,4 +227,37 @@ class InfiniteSystem(System):
assert shape[0] == shape[1]
# Subtract energy from the diagonal.
ham.flat[::ham.shape[0] + 1] -= energy
- return physics.self_energy(ham, self.inter_slice_hopping(args=args))
+ return physics.selfenergy(ham, self.inter_slice_hopping(args=args))
+
+
+class PrecalculatedLead(object):
+ def __init__(self, modes=None, selfenergy=None):
+ """A general lead defined by its self energy.
+
+ Parameters
+ ----------
+ modes : kwant.physics.ModesTuple
+ Modes of the lead.
+ selfenergy : numpy array
+ Lead self-energy.
+
+ Notes
+ -----
+ At least one of `modes` and `selfenergy` must be provided.
+ """
+ if modes is None and selfenergy is None:
+ raise ValueError("No precalculated values provided.")
+ self._modes = modes
+ self._selfenergy = selfenergy
+
+ def modes(self, energy, args=()):
+ if self._modes is not None:
+ return self._modes
+ else:
+ raise ValueError("No precalculated modes were provided.")
+
+ def selfenergy(self, energy, args=()):
+ if self._selfenergy is not None:
+ return self._selfenergy
+ else:
+ raise ValueError("No precalculated self-energy was provided.")
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index 72d1caf7921307c3b9c657bb655166a4ad5d19c5..93f6b9585014cf6a73e35603cee85b816e6edd7d 100644
--- a/kwant/tests/test_system.py
+++ b/kwant/tests/test_system.py
@@ -56,6 +56,19 @@ def test_hamiltonian_submatrix():
mat_sp = sys2.hamiltonian_submatrix(sparse=True).todense()
np.testing.assert_array_equal(mat_sp, mat_dense)
+ # Test precalculation of modes.
+ np.random.seed(5)
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
+ lead[gr(0)] = np.zeros((2, 2))
+ lead[gr(0), gr(1)] = np.random.randn(2, 2)
+ sys.attach_lead(lead)
+ sys2 = sys.finalized()
+ smatrix = kwant.solve(sys2, .1).data
+ sys3 = sys2.precalculate(.1, calculate_selfenergy=False)
+ smatrix2 = kwant.solve(sys3, .1).data
+ np.testing.assert_almost_equal(smatrix, smatrix2)
+ assert_raises(ValueError, kwant.solve, sys3, 0.2, None, None, True)
+
# Test for shape errors.
sys[gr(0), gr(2)] = np.array([[1, 2]])
sys2 = sys.finalized()
@@ -86,3 +99,4 @@ def test_hamiltonian_submatrix():
mat = mat[perm, :]
mat = mat[:, perm]
np.testing.assert_array_equal(mat, mat_should_be)
+