kwant/continuum/__init__.py

@@ -10,11 +10,12 @@ try:
     from .discretizer import discretize, discretize_symbolic, build_discretized
     from ._common import sympify, lambdify
     from ._common import momentum_operators, position_operators
+    from .landau_levels import to_landau_basis, discretize_landau, LandauLattice
 except ImportError as error:
     msg = ("'kwant.continuum' is not available because one or more of its "
            "dependencies is not installed.")
     raise ImportError(msg) from error
 
-
 __all__ = ['discretize', 'discretize_symbolic', 'build_discretized',
+           'to_landau_basis', 'discretize_landau', 'LandauLattice',
            'sympify', 'lambdify', 'momentum_operators', 'position_operators']

kwant/continuum/discretizer.py

-# Copyright 2011-2017 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -20,7 +20,7 @@ from sympy.printing.lambdarepr import LambdaPrinter
 from sympy.printing.precedence import precedence
 from sympy.core.function import AppliedUndef
 
-from .. import builder, lattice
+from .. import builder, lattice, system
 from .. import KwantDeprecationWarning
 from .._common import reraise_warnings
 from ._common import (sympify, gcd, position_operators, momentum_operators,
@@ -63,7 +63,7 @@ class _DiscretizedBuilder(builder.Builder):
 
         for key, val in itertools.chain(self.site_value_pairs(),
                                         self.hopping_value_pairs()):
-            if isinstance(key, builder.Site):
+            if isinstance(key, system.Site):
                 result.append("# Onsite element:\n")
             else:
                 a, b = key

kwant/continuum/landau_levels.py

+# Copyright 2011-2019 Kwant authors.
+#
+# This file is part of Kwant.  It is subject to the license terms in the file
+# LICENSE.rst 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 file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from keyword import iskeyword
+import functools
+import operator
+import collections
+
+from math import sqrt
+import tinyarray as ta
+import numpy as np
+import sympy
+
+import kwant.lattice
+import kwant.builder
+
+import kwant.continuum
+import kwant.continuum._common
+import kwant.continuum.discretizer
+from kwant.continuum import momentum_operators, position_operators
+
+__all__ = ["to_landau_basis", "discretize_landau"]
+
+coordinate_vectors = dict(zip("xyz", np.eye(3)))
+
+ladder_lower, ladder_raise = sympy.symbols(r"a a^\dagger", commutative=False)
+
+
+def to_landau_basis(hamiltonian, momenta=None):
+    r"""Replace two momenta by Landau level ladder operators.
+
+    Replaces:
+
+        k_0 -> sqrt(B/2) * (a + a^\dagger)
+        k_1 -> 1j * sqrt(B/2) * (a - a^\dagger)
+
+    Parameters
+    ----------
+    hamiltonian : str or sympy expression
+        The Hamiltonian to which to apply the Landau level transformation.
+    momenta : sequence of str (optional)
+        The momenta to replace with Landau level ladder operators. If not
+        provided, 'k_x' and 'k_y' are used
+
+    Returns
+    -------
+    hamiltonian : sympy expression
+    momenta : sequence of sympy atoms
+        The momentum operators that have been replaced by ladder operators.
+    normal_coordinate : sympy atom
+        The remaining position coordinate. May or may not be present
+        in 'hamiltonian'.
+    """
+    hamiltonian = kwant.continuum.sympify(hamiltonian)
+    momenta = _normalize_momenta(momenta)
+    normal_coordinate = _find_normal_coordinate(hamiltonian, momenta)
+
+    # Substitute ladder operators for Landau level momenta
+    B = sympy.symbols("B")
+    hamiltonian = hamiltonian.subs(
+        {
+            momenta[0]: sympy.sqrt(abs(B) / 2) * (ladder_raise + ladder_lower),
+            momenta[1]: sympy.I
+            * sympy.sqrt(abs(B) / 2)
+            * (ladder_lower - ladder_raise),
+        }
+    )
+
+    return hamiltonian, momenta, normal_coordinate
+
+
+def discretize_landau(hamiltonian, N, momenta=None, grid_spacing=1):
+    """Discretize a Hamiltonian in a basis of Landau levels.
+
+    Parameters
+    ----------
+    hamiltonian : str or sympy expression
+    N : positive integer
+        The number of Landau levels in the basis.
+    momenta : sequence of str (optional)
+        The momenta defining the plane perpendicular to the magnetic field.
+        If not provided, "k_x" and "k_y" are used.
+    grid_spacing : float, default: 1
+        The grid spacing to use when discretizing the normal direction
+        (parallel to the magnetic field).
+
+    Returns
+    -------
+    builder : `~kwant.builder.Builder`
+        'hamiltonian' discretized in a basis of Landau levels in the plane
+        defined by 'momenta'. If a third coordinate is present in 'hamiltonian',
+        then this also has a translational symmetry in that coordinate direction.
+        The builder has a parameter 'B' in addition to any other parameters
+        present in the provided 'hamiltonian'.
+
+    Notes
+    -----
+    The units of magnetic field are :math:`ϕ₀ / 2 π a²` with :math:`ϕ₀ = h/e`
+    the magnetic flux quantum and :math:`a` the unit length.
+    """
+
+    if N <= 0:
+        raise ValueError("N must be positive")
+
+    hamiltonian, momenta, normal_coordinate = to_landau_basis(hamiltonian, momenta)
+
+    # Discretize in normal direction and split terms for onsites/hoppings into
+    # monomials in ladder operators.
+    tb_hamiltonian, _ = kwant.continuum.discretize_symbolic(
+        hamiltonian, coords=[normal_coordinate.name]
+    )
+    tb_hamiltonian = {
+        key: kwant.continuum._common.monomials(value, gens=(ladder_lower, ladder_raise))
+        for key, value in tb_hamiltonian.items()
+    }
+
+    # Replace ladder operator monomials by tuple of integers:
+    # e.g. a^\dagger a^2 -> (+1, -2).
+    tb_hamiltonian = {
+        outer_key: {
+            _ladder_term(inner_key): inner_value
+            for inner_key, inner_value in outer_value.items()
+        }
+        for outer_key, outer_value in tb_hamiltonian.items()
+    }
+
+    # Construct map from LandauLattice HoppingKinds to a sequence of pairs
+    # that encode the ladder operators and multiplying expression.
+    tb_hoppings = collections.defaultdict(list)
+    for outer_key, outer_value in tb_hamiltonian.items():
+        for inner_key, inner_value in outer_value.items():
+            tb_hoppings[(*outer_key, sum(inner_key))] += [(inner_key, inner_value)]
+    # Extract the number of orbitals on each site/hopping
+    random_element = next(iter(tb_hoppings.values()))[0][1]
+    norbs = 1 if isinstance(random_element, sympy.Expr) else random_element.shape[0]
+    tb_onsite = tb_hoppings.pop((0, 0), None)
+
+    # Construct Builder
+    if _has_coordinate(normal_coordinate, hamiltonian):
+        sym = kwant.lattice.TranslationalSymmetry([grid_spacing, 0])
+    else:
+        sym = kwant.builder.NoSymmetry()
+    lat = LandauLattice(grid_spacing, norbs=norbs)
+    syst = kwant.Builder(sym)
+
+    # Add onsites
+    landau_sites = (lat(0, j) for j in range(N))
+    if tb_onsite is None:
+        syst[landau_sites] = ta.zeros((norbs, norbs))
+    else:
+        syst[landau_sites] = _builder_value(
+            tb_onsite, normal_coordinate.name, grid_spacing, is_onsite=True
+        )
+
+    # Add zero hoppings between adjacent Landau levels.
+    # Necessary to be able to use the Landau level builder
+    # to populate another builder using builder.fill().
+    syst[kwant.builder.HoppingKind((0, 1), lat)] = ta.zeros((norbs, norbs))
+
+    # Add the hoppings from the Hamiltonian
+    for hopping, parts in tb_hoppings.items():
+        syst[kwant.builder.HoppingKind(hopping, lat)] = _builder_value(
+            parts, normal_coordinate.name, grid_spacing, is_onsite=False
+        )
+
+    return syst
+
+
+# This has to subclass lattice so that it will work with TranslationalSymmetry.
+class LandauLattice(kwant.lattice.Monatomic):
+    """
+    A `~kwant.lattice.Monatomic` lattice with a Landau level index per site.
+
+    Site tags (see `~kwant.system.SiteFamily`) are pairs of integers, where
+    the first integer describes the real space position and the second the
+    Landau level index.
+
+    The real space Bravais lattice is one dimensional, oriented parallel
+    to the magnetic field. The Landau level index represents the harmonic
+    oscillator basis states used for the Landau quantization in the plane
+    normal to the field.
+
+    Parameters
+    ----------
+    grid_spacing : float
+        Real space lattice spacing (parallel to the magnetic field).
+    offset : float (optional)
+        Displacement of the lattice origin from the real space
+        coordinates origin.
+    """
+
+    def __init__(self, grid_spacing, offset=None, name="", norbs=None):
+        if offset is not None:
+            offset = [offset, 0]
+        # The second vector and second coordinate do not matter (they are
+        # not used in pos())
+        super().__init__([[grid_spacing, 0], [0, 1]], offset, name, norbs)
+
+    def pos(self, tag):
+        return ta.array((self.prim_vecs[0, 0] * tag[0] + self.offset[0],))
+
+    def landau_index(self, tag):
+        return tag[-1]
+
+
+def _builder_value(terms, normal_coordinate, grid_spacing, is_onsite):
+    """Construct an onsite/hopping value function from a list of terms
+
+    Parameters
+    ----------
+    terms : list
+        Each element is a pair (ladder_term, sympy expression/matrix).
+        ladder_term is a tuple of integers that encodes a string of
+        Landau raising/lowering operators and the sympy expression
+        is the rest
+    normal_coordinate : str
+    grid_spacing : float
+        The grid spacing in the normal direction
+    is_onsite : bool
+        True if we are constructing an onsite value function
+    """
+    ladder_term_symbols = [sympy.Symbol(_ladder_term_name(lt)) for lt, _ in terms]
+    # Construct a single expression from the terms, multiplying each part
+    # by the placeholder that represents the prefactor from the ladder operator term.
+    (ladder, (_, part)), *rest = zip(ladder_term_symbols, terms)
+    expr = ladder * part
+    for ladder, (_, part) in rest:
+        expr += ladder * part
+    expr = expr.subs(
+        {kwant.continuum.discretizer._displacements[normal_coordinate]: grid_spacing}
+    )
+    # Construct the return string and temporary variable names
+    # for function calls.
+    return_string, map_func_calls, const_symbols, _cache = kwant.continuum.discretizer._return_string(
+        expr, coords=[normal_coordinate]
+    )
+
+    # Remove the ladder term placeholders, as these are not parameters
+    const_symbols = set(const_symbols)
+    for ladder_term in ladder_term_symbols:
+        const_symbols.discard(ladder_term)
+
+    # Construct the argument part of signature. Arguments
+    # consist of any constants and functions in the return string.
+    arg_names = set.union(
+        {s.name for s in const_symbols}, {str(k.func) for k in map_func_calls}
+    )
+    arg_names.discard("Abs")  # Abs function is not a parameter
+    for arg_name in arg_names:
+        if not (arg_name.isidentifier() and not iskeyword(arg_name)):
+            raise ValueError(
+                "Invalid name in used symbols: {}\n"
+                "Names of symbols used in Hamiltonian "
+                "must be valid Python identifiers and "
+                "may not be keywords".format(arg_name)
+            )
+    arg_names = ", ".join(sorted(arg_names))
+    # Construct site part of the function signature
+    site_string = "from_site" if is_onsite else "to_site, from_site"
+    # Construct function signature
+    if arg_names:
+        function_header = "def _({}, {}):".format(site_string, arg_names)
+    else:
+        function_header = "def _({}):".format(site_string)
+    # Construct function body
+    function_body = []
+    if "B" not in arg_names:
+        # B is not a parameter for terms with no ladder operators but we still
+        # need something to pass to _evaluate_ladder_term
+        function_body.append("B = +1")
+    function_body.extend(
+        [
+            "{}, = from_site.pos".format(normal_coordinate),
+            "_ll_index = from_site.family.landau_index(from_site.tag)",
+        ]
+    )
+    # To get the correct hopping if B < 0, we need to set the Hermitian
+    # conjugate of the ladder operator matrix element, which swaps the
+    # from_site and to_site Landau level indices.
+    if not is_onsite:
+        function_body.extend(
+            ["if B < 0:",
+             "    _ll_index = to_site.family.landau_index(to_site.tag)"
+            ]
+        )
+    function_body.extend(
+        "{} = _evaluate_ladder_term({}, _ll_index, B)".format(symb.name, lt)
+        for symb, (lt, _) in zip(ladder_term_symbols, terms)
+    )
+    function_body.extend(
+        "{} = {}".format(v, kwant.continuum.discretizer._print_sympy(k))
+        for k, v in map_func_calls.items()
+    )
+    function_body.append(return_string)
+    func_code = "\n    ".join([function_header] + function_body)
+    # Add "I" to namespace just in case sympy again would miss to replace it
+    # with Python's 1j as it was the case with SymPy 1.2 when "I" was argument
+    # of some function.
+    namespace = {
+        "pi": np.pi,
+        "I": 1j,
+        "_evaluate_ladder_term": _evaluate_ladder_term,
+        "Abs": abs,
+    }
+    namespace.update(_cache)
+    # Construct full source, including cached arrays
+    source = []
+    for k, v in _cache.items():
+        source.append("{} = (\n{})\n".format(k, repr(np.array(v))))
+    source.append(func_code)
+    exec(func_code, namespace)
+    f = namespace["_"]
+    f._source = "".join(source)
+
+    return f
+
+
+def _ladder_term(operator_string):
+    r"""Return a tuple of integers representing a string of ladder operators
+
+    Parameters
+    ----------
+    operator_string : Sympy expression
+        Monomial in ladder operators, e.g. a^\dagger a^2 a^\dagger.
+
+    Returns
+    -------
+    ladder_term : tuple of int
+        e.g. a^\dagger a^2 -> (+1, -2)
+    """
+    ret = []
+    for factor in operator_string.as_ordered_factors():
+        ladder_op, exponent = factor.as_base_exp()
+        if ladder_op == ladder_lower:
+            sign = -1
+        elif ladder_op == ladder_raise:
+            sign = +1
+        else:
+            sign = 0
+        ret.append(sign * int(exponent))
+    return tuple(ret)
+
+
+def _ladder_term_name(ladder_term):
+    """
+    Parameters
+    ----------
+    ladder_term : tuple of int
+
+    Returns
+    -------
+    ladder_term_name : str
+    """
+
+    def ns(i):
+        if i >= 0:
+            return str(i)
+        else:
+            return "_" + str(-i)
+
+    return "_ladder_{}".format("_".join(map(ns, ladder_term)))
+
+
+def _evaluate_ladder_term(ladder_term, n, B):
+    r"""Evaluates the prefactor for a ladder operator on a landau level.
+
+    Example: a^\dagger a^2 -> (n - 1) * sqrt(n)
+
+    Parameters
+    ----------
+    ladder_term : tuple of int
+        Represents a string of ladder operators. Positive
+        integers represent powers of the raising operator,
+        negative integers powers of the lowering operator.
+    n : non-negative int
+        Landau level index on which to act with ladder_term.
+    B : float
+        Magnetic field with sign
+
+    Returns
+    -------
+    ladder_term_prefactor : float
+    """
+    assert n >= 0
+    # For negative B we swap a and a^\dagger.
+    if B < 0:
+        ladder_term = tuple(-i for i in ladder_term)
+    ret = 1
+    for m in reversed(ladder_term):
+        if m > 0:
+            factors = range(n + 1, n + m + 1)
+        elif m < 0:
+            factors = range(n + m + 1, n + 1)
+            if n == 0:
+                return 0  # a|0> = 0
+        else:
+            factors = (1,)
+        ret *= sqrt(functools.reduce(operator.mul, factors))
+        n += m
+    return ret
+
+
+def _normalize_momenta(momenta=None):
+    """Return Landau level momenta as Sympy atoms
+
+    Parameters
+    ----------
+    momenta : None or pair of int or pair of str
+        The momenta to choose. If None then 'k_x' and 'k_y'
+        are chosen. If integers, then these are the indices
+        of the momenta: 0 → k_x, 1 → k_y, 2 → k_z. If strings,
+        then these name the momenta.
+
+    Returns
+    -------
+    The specified momenta as sympy atoms.
+    """
+
+    # Specify which momenta to substitute for the Landau level basis.
+    if momenta is None:
+        # Use k_x and k_y by default
+        momenta = momentum_operators[:2]
+    else:
+        if len(momenta) != 2:
+            raise ValueError("Two momenta must be specified.")
+
+        k_names = [k.name for k in momentum_operators]
+
+        if all([type(i) is int for i in momenta]) and all(
+            [i >= 0 and i < 3 for i in momenta]
+        ):
+            momenta = [momentum_operators[i] for i in momenta]
+        elif all([isinstance(momentum, str) for momentum in momenta]) and all(
+            [momentum in k_names for momentum in momenta]
+        ):
+            momenta = [
+                momentum_operators[k_names.index(momentum)] for momentum in momenta
+            ]
+        else:
+            raise ValueError("Momenta must all be integers or strings.")
+
+    return tuple(momenta)
+
+
+def _find_normal_coordinate(hamiltonian, momenta):
+    discrete_momentum = next(
+        momentum for momentum in momentum_operators if momentum not in momenta
+    )
+    normal_coordinate = position_operators[momentum_operators.index(discrete_momentum)]
+    return normal_coordinate
+
+
+def _has_coordinate(coord, expr):
+    momentum = momentum_operators[position_operators.index(coord)]
+    atoms = set(expr.atoms())
+    return coord in atoms or momentum in atoms

kwant/continuum/tests/test_discretizer.py

@@ -548,14 +548,14 @@ def test_grid(ham, grid_spacing, grid):
 
 
 @pytest.mark.parametrize('ham, grid_offset, offset, norbs', [
-    ('k_x', None, 0, None),
+    ('k_x', None, 0, 1),
     ('k_x', None, 0, 1),
     ('k_x * eye(2)', None, 0, 2),
-    ('k_x', (0,), 0, None),
-    ('k_x', (1,), 1, None),
-    ('k_x + k_y', None, (0, 0), None),
-    ('k_x + k_y', (0, 0), (0, 0), None),
-    ('k_x + k_y', (1, 2), (1, 2), None),
+    ('k_x', (0,), 0, 1),
+    ('k_x', (1,), 1, 1),
+    ('k_x + k_y', None, (0, 0), 1),
+    ('k_x + k_y', (0, 0), (0, 0), 1),
+    ('k_x + k_y', (1, 2), (1, 2), 1),
 ])
 def test_grid_input(ham, grid_offset, offset, norbs):
     # build appriopriate grid
@@ -579,7 +579,7 @@ def test_grid_input(ham, grid_offset, offset, norbs):
 
 def test_grid_offset_passed_to_functions():
     V = lambda x: x
-    grid = Monatomic([[1, ]], offset=[0.5, ])
+    grid = Monatomic([[1, ]], offset=[0.5, ], norbs=1)
     tb = discretize('V(x)', 'x', grid=grid)
     onsite = tb[tb.lattice(0)]
     bools = [np.allclose(onsite(tb.lattice(i), V), V(tb.lattice(i).pos))
@@ -588,11 +588,11 @@ def test_grid_offset_passed_to_functions():
 
 
 @pytest.mark.parametrize("ham, coords, grid", [
-    ("k_x", None, Monatomic([[1, 0]])),
-    ("k_x", 'xy', Monatomic([[1, 0]])),
+    ("k_x", None, Monatomic([[1, 0]], norbs=1)),
+    ("k_x", 'xy', Monatomic([[1, 0]], norbs=1)),
     ("k_x", None, Monatomic([[1, ]], norbs=2)),
     ("k_x * eye(2)", None, Monatomic([[1, ]], norbs=1)),
-    ("k_x+k_y", None, Monatomic([[1, 0], [1, 1]])),
+    ("k_x+k_y", None, Monatomic([[1, 0], [1, 1]], norbs=1)),
 ])
 def test_grid_constraints(ham, coords, grid):
     with pytest.raises(ValueError):
@@ -606,7 +606,7 @@ def test_check_symbol_names(name):
 
 
 def test_rectangular_grid():
-    lat = Monatomic([[1, 0], [0, 2]])
+    lat = Monatomic([[1, 0], [0, 2]], norbs=1)
 
     tb = discretize("V(x, y)", 'xy', grid=lat)
     assert np.allclose(tb[lat(0, 0)](lat(1, 0), lambda x, y: x), 1)

kwant/continuum/tests/test_landau_levels.py

+# Copyright 2011-2019 Kwant authors.
+#
+# This file is part of Kwant.  It is subject to the license terms in the file
+# LICENSE.rst 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 file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from math import sqrt
+
+import numpy as np
+import sympy
+import pytest
+import itertools
+
+import kwant.builder
+import kwant.lattice
+
+from .._common import position_operators, momentum_operators, sympify
+from ..landau_levels import (
+    to_landau_basis,
+    discretize_landau,
+    _ladder_term,
+    _evaluate_ladder_term,
+    ladder_lower,
+    ladder_raise,
+    LandauLattice,
+)
+
+x, y, z = position_operators
+k_x, k_y, k_z = momentum_operators
+B = sympy.symbols("B")
+V = sympy.symbols("V", cls=sympy.Function)
+a, a_dag = ladder_lower, ladder_raise
+
+
+def test_ladder_term():
+    assert _ladder_term(a ** 2 * a_dag) == (-2, 1)
+    assert _ladder_term(a_dag ** 5 * a ** 3 * a_dag) == (5, -3, 1)
+    # non-ladder operators give a shift of 0
+    assert _ladder_term(B * a ** 2) == (0, -2)
+
+
+def test_evaluate_ladder_term():
+    assert np.isclose(_evaluate_ladder_term((+1, -1), 1, +1), 1)
+    assert np.isclose(_evaluate_ladder_term((+1, -1), 2, +1), 2)
+    assert np.isclose(_evaluate_ladder_term((-2, +3, -2), 5, +1), 4 * 5 * 6 * sqrt(5))
+    # annihilating |0> is always 0
+    assert _evaluate_ladder_term((-1,), 0, +1) == 0
+    assert _evaluate_ladder_term((-2,), 1, +1) == 0
+    assert _evaluate_ladder_term((-3,), 1, +1) == 0
+    assert _evaluate_ladder_term((+3, -2), 1, +1) == 0
+    assert _evaluate_ladder_term((-3, -2, +3), 1, +1) == 0
+    # negative B swaps creation and annihilation operators
+    assert _evaluate_ladder_term((+1, -1), 2, +1) == _evaluate_ladder_term(
+        (-1, +1), 2, -1
+    )
+    assert _evaluate_ladder_term((-2, +3, -2), 5, +1) == _evaluate_ladder_term(
+        (+2, -3, +2), 5, -1
+    )
+
+
+def test_to_landau_basis():
+    # test basic usage
+    ham, momenta, normal_coord = to_landau_basis("k_x**2 + k_y**2")
+    assert sympy.expand(ham) == abs(B) * a * a_dag + abs(B) * a_dag * a
+    assert momenta == (k_x, k_y)
+    assert normal_coord == z
+
+    # test that hamiltonian can be specified as a sympy expression
+    ham, momenta, normal_coord = to_landau_basis(sympify("k_x**2 + k_y**2"))
+    assert sympy.expand(ham) == abs(B) * a * a_dag + abs(B) * a_dag * a
+    assert momenta == (k_x, k_y)
+    assert normal_coord == z
+
+    # test that
+    ham, momenta, normal_coord = to_landau_basis("k_x**2 + k_y**2 + k_z**2 + V(z)")
+    assert sympy.expand(ham) == (
+        abs(B) * a * a_dag + abs(B) * a_dag * a + k_z ** 2 + V(z)
+    )
+    assert momenta == (k_x, k_y)
+    assert normal_coord == z
+
+    # test for momenta explicitly specified
+    ham, momenta, normal_coord = to_landau_basis(
+        "k_x**2 + k_y**2 + k_z**2 + k_x*k_y", momenta=("k_z", "k_y")
+    )
+    assert sympy.expand(ham) == (
+        abs(B) * a * a_dag
+        + abs(B) * a_dag * a
+        + k_x ** 2
+        + sympy.I * sympy.sqrt(abs(B) / 2) * k_x * a
+        - sympy.I * sympy.sqrt(abs(B) / 2) * k_x * a_dag
+    )
+    assert momenta == (k_z, k_y)
+    assert normal_coord == x
+
+
+def test_discretize_landau():
+    n_levels = 10
+    magnetic_field = 1 / 3  # a suitably arbitrary value
+    # Ensure that we can handle products of ladder operators by truncating
+    # several levels higher than the number of levels we actually want.
+    a = np.diag(np.sqrt(np.arange(1, n_levels + 5)), k=1)
+    a_dag = a.conjugate().transpose()
+    k_x = sqrt(magnetic_field / 2) * (a + a_dag)
+    k_y = 1j * sqrt(magnetic_field / 2) * (a - a_dag)
+    sigma_0 = np.eye(2)
+    sigma_x = np.array([[0, 1], [1, 0]])
+    sigma_y = np.array([[0, -1j], [1j, 0]])
+
+    # test that errors are raised on invalid input
+    with pytest.raises(ValueError):
+        discretize_landau("k_x**2 + k_y**2", N=0)
+    with pytest.raises(ValueError):
+        discretize_landau("k_x**2 + k_y**2", N=-1)
+
+    # test a basic Hamiltonian with no normal coordinate dependence
+    syst = discretize_landau("k_x**2 + k_y**2", N=n_levels)
+    lat = LandauLattice(1, norbs=1)
+    assert isinstance(syst.symmetry, kwant.builder.NoSymmetry)
+    syst = syst.finalized()
+    assert set(syst.sites) == {lat(0, j) for j in range(n_levels)}
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=0)), np.zeros((n_levels, n_levels))
+    )
+    should_be = k_x @ k_x + k_y @ k_y
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        should_be[:n_levels, :n_levels],
+    )
+
+    # test negative magnetic field
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=-magnetic_field)),
+        should_be[:n_levels, :n_levels],
+    )
+
+    # test hamiltonian with no onsite elements
+    syst = discretize_landau("k_x", N=n_levels)
+    syst = syst.finalized()
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        k_x[:n_levels, :n_levels],
+    )
+
+    # test a basic Hamiltonian with normal coordinate dependence
+    grid = 1 / 7  # a suitably arbitrary value
+    syst = discretize_landau(
+        "k_x**2 + k_y**2 + k_z**2 + V(z)", N=n_levels, grid_spacing=grid
+    )
+    assert isinstance(syst.symmetry, kwant.lattice.TranslationalSymmetry)
+    syst = syst.finalized()
+    zero_potential = syst.cell_hamiltonian(params=dict(B=magnetic_field, V=lambda z: 0))
+    with_potential = syst.cell_hamiltonian(params=dict(B=magnetic_field, V=lambda z: 1))
+    # extra +2/grid**2 from onsite kinetic energy
+    assert np.allclose(
+        zero_potential,
+        should_be[:n_levels, :n_levels] + (2 / grid ** 2) * np.eye(n_levels),
+    )
+    # V(z) just adds the same energy to each onsite
+    assert np.allclose(with_potential - zero_potential, np.eye(n_levels))
+    # hopping matrix does not exchange landau levels
+    assert np.allclose(
+        syst.inter_cell_hopping(params=dict(B=magnetic_field, V=lambda z: 0)),
+        -np.eye(n_levels) / grid ** 2,
+    )
+
+    # test a Hamiltonian with mixing between Landau levels
+    # and spatial degrees of freedom.
+    syst = discretize_landau("k_x**2 + k_y**2 + k_x*k_z", N=n_levels)
+    syst = syst.finalized()
+    assert np.allclose(
+        syst.inter_cell_hopping(params=dict(B=magnetic_field)),
+        -1j * k_x[:n_levels, :n_levels] / 2,
+    )
+
+    # test a Hamiltonian with extra degrees of freedom
+    syst = discretize_landau("sigma_0 * k_x**2 + sigma_x * k_y**2", N=n_levels)
+    syst = syst.finalized()
+    assert syst.sites[0].family.norbs == 2
+    should_be = np.kron(k_x @ k_x, sigma_0) + np.kron(k_y @ k_y, sigma_x)
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        should_be[: 2 * n_levels, : 2 * n_levels],
+    )
+
+    # test a linear Hamiltonian
+    syst = discretize_landau("sigma_y * k_x - sigma_x * k_y", N=n_levels)
+    syst = syst.finalized()
+    should_be = np.kron(k_x, sigma_y) - np.kron(k_y, sigma_x)
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        should_be[: 2 * n_levels, : 2 * n_levels],
+    )
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        syst.hamiltonian_submatrix(params=dict(B=-magnetic_field)),
+    )
+
+
+def test_analytical_spectrum():
+    hamiltonian = """(k_x**2 + k_y**2) * sigma_0 +
+                    alpha * (k_x * sigma_y - k_y * sigma_x) +
+                    g * B * sigma_z"""
+
+    def exact_Es(n, B, alpha, g):
+        # See e.g. R. Winkler (2003), section 8.4.1
+        sign_B = np.sign(B)
+        B = np.abs(B)
+        Ep = 2*B*(n+1) - 0.5*np.sqrt((2*B - sign_B*2*g*B)**2 + 8*B*alpha**2*(n+1))
+        Em = 2*B*n + 0.5*np.sqrt((2*B - sign_B*2*g*B)**2 + 8*B*alpha**2*n)
+        return Ep, Em
+
+    N = 20
+    syst = discretize_landau(hamiltonian, N)
+    syst = syst.finalized()
+    params = dict(alpha=0.07, g=0.04)
+    for _ in range(5):
+        B = 0.01 + np.random.rand()*3
+        params['B'] = B
+        exact = [exact_Es(n, B, params['alpha'], params['g']) for n in range(N)]
+        # We only check the first N levels - the SOI couples adjacent levels,
+        # so the higher ones are not necessarily determined accurately in the
+        # discretization
+        exact = np.sort([energy for energies in exact for energy in energies])[:N]
+        ll_spect = np.linalg.eigvalsh(syst.hamiltonian_submatrix(params=params))[:len(exact)]
+        assert np.allclose(ll_spect, exact)
+
+
+
+def test_fill():
+
+    def shape(site, lower, upper):
+        (z, )= site.pos
+        return lower <= z < upper
+
+    hamiltonian = "k_x**2 + k_y**2 + k_z**2"
+    N = 6
+    template = discretize_landau(hamiltonian, N)
+
+    syst = kwant.Builder()
+    width = 4
+    syst.fill(template, lambda site: shape(site, 0, width), (0, ));
+
+    correct_tags = [(coordinate, ll_index) for coordinate, ll_index
+                    in itertools.product(range(width), range(N))]
+
+    syst_tags = [site.tag for site in syst.sites()]
+
+    assert len(syst_tags) == len(correct_tags)
+    assert all(tag in correct_tags for tag in syst_tags)

kwant/graph/__init__.py

@@ -9,8 +9,7 @@
 """Functionality for graphs"""
 
 # Merge the public interface of all submodules.
-__all__ = []
-for module in ['core', 'defs']:
-    exec('from . import {0}'.format(module))
-    exec('from .{0} import *'.format(module))
-    exec('__all__.extend({0}.__all__)'.format(module))
+from .core import *
+from .defs import *
+
+__all__ = [core.__all__ + defs.__all__]

kwant/kpm.py

 # -*- coding: utf-8 -*-
-# Copyright 2011-2016 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -731,23 +731,25 @@ class Correlator:
         g_kernel[0] /= 2
         mu_kernel = np.outer(g_kernel, g_kernel) * self.moments_matrix
 
-        e = (self.energies - self._b) / self._a
+        e_scaled = (self.energies - self._b) / self._a
 
-        # arrays for faster calculation
-        sqrt_e = np.sqrt(1 - e ** 2)
-        arccos_e = np.arccos(e)
+        m_array = np.arange(n_moments)
+        def _integral_factor(e):
+            # arrays for faster calculation
+            sqrt_e = np.sqrt(1 - e ** 2)
+            arccos_e = np.arccos(e)
 
-        exp_n = np.exp(1j * np.outer(arccos_e, np.arange(n_moments)))
-        t_n = np.real(exp_n)
+            exp_n = np.exp(1j * arccos_e * m_array)
+            t_n = np.real(exp_n)
 
-        e_plus = (np.outer(e, np.ones(n_moments)) -
-                  1j * np.outer(sqrt_e, np.arange(n_moments)))
-        e_plus = e_plus * exp_n
+            e_plus = (e - 1j * sqrt_e * m_array)
+            e_plus = e_plus * exp_n
 
-        big_gamma = e_plus[:, None, :] * t_n[:, :, None]
-        big_gamma += big_gamma.conj().swapaxes(1, 2)
-
-        self._integral_factor = np.tensordot(mu_kernel, big_gamma.T)
+            big_gamma = e_plus[None, :] * t_n[:, None]
+            big_gamma += big_gamma.conj().T
+            return np.tensordot(mu_kernel, big_gamma.T)
+        self._integral_factor = np.array([_integral_factor(e)
+                                          for e in e_scaled]).T
 
 
 def conductivity(hamiltonian, alpha='x', beta='x', positions=None, **kwargs):
@@ -926,7 +928,7 @@ def RandomVectors(syst, where=None, rng=None):
     ----------
     syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
         If a system is passed, it should contain no leads.
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
         Spatial range of the vectors produced. If ``syst`` is a
         `~kwant.system.FiniteSystem`, where behaves as in
         `~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
@@ -948,7 +950,7 @@ class LocalVectors:
     ----------
     syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
         If a system is passed, it should contain no leads.
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
         Spatial range of the vectors produced. If ``syst`` is a
         `~kwant.system.FiniteSystem`, where behaves as in
         `~kwant.operator.Density`. If ``syst`` is a matrix, ``where``

kwant/lattice.py

-# Copyright 2011-2013 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -13,7 +13,7 @@ from math import sqrt
 from itertools import product
 import numpy as np
 import tinyarray as ta
-from . import builder
+from . import builder, system
 from .linalg import lll
 from ._common import ensure_isinstance
 
@@ -70,7 +70,7 @@ class Polyatomic:
     A Bravais lattice with an arbitrary number of sites in the basis.
 
     Contains `Monatomic` sublattices.  Note that an instance of ``Polyatomic`` is
-    not itself a `~kwant.builder.SiteFamily`, only its sublattices are.
+    not itself a `~kwant.system.SiteFamily`, only its sublattices are.
 
     Parameters
     ----------
@@ -171,7 +171,7 @@ class Polyatomic:
         >>> syst[lat.neighbors()] = 1
         """
         def shape_sites(symmetry=None):
-            Site = builder.Site
+            Site = system.Site
 
             if symmetry is None:
                 symmetry = builder.NoSymmetry()
@@ -306,7 +306,7 @@ class Polyatomic:
         """
         # This algorithm is not designed to be fast and can be improved,
         # however there is no real need.
-        Site = builder.Site
+        Site = system.Site
         sls = self.sublattices
         shortest_hopping = sls[0].n_closest(
             sls[0].pos(([0] * sls[0].lattice_dim)), 2)[-1]
@@ -396,12 +396,12 @@ def short_array_str(array):
     return full[1 : -1]
 
 
-class Monatomic(builder.SiteFamily, Polyatomic):
+class Monatomic(system.SiteFamily, Polyatomic):
     """
     A Bravais lattice with a single site in the basis.
 
-    Instances of this class provide the `~kwant.builder.SiteFamily` interface.
-    Site tags (see `~kwant.builder.SiteFamily`) are sequences of integers and
+    Instances of this class provide the `~kwant.system.SiteFamily` interface.
+    Site tags (see `~kwant.system.SiteFamily`) are sequences of integers and
     describe the lattice coordinates of a site.
 
     ``Monatomic`` instances are used as site families on their own or as
@@ -470,6 +470,13 @@ class Monatomic(builder.SiteFamily, Polyatomic):
     def __str__(self):
         return self.cached_str
 
+    def normalize_tags(self, tags):
+        tags = np.asarray(tags, int)
+        tags.flags.writeable = False
+        if tags.shape[1] != self.lattice_dim:
+            raise ValueError("Dimensionality mismatch.")
+        return tags
+
     def normalize_tag(self, tag):
         tag = ta.array(tag, int)
         if len(tag) != self.lattice_dim:
@@ -495,6 +502,10 @@ class Monatomic(builder.SiteFamily, Polyatomic):
         """
         return ta.array(self.n_closest(pos)[0])
 
+    def positions(self, tags):
+        """Return the real-space positions of the sites with the given tags."""
+        return tags @ self._prim_vecs + self.offset
+
     def pos(self, tag):
         """Return the real-space position of the site with a given tag."""
         return ta.dot(tag, self._prim_vecs) + self.offset
@@ -687,26 +698,48 @@ class TranslationalSymmetry(builder.Symmetry):
 
     def which(self, site):
         det_x_inv_m_part, det_m = self._get_site_family_data(site.family)[-2:]
-        result = ta.dot(det_x_inv_m_part, site.tag) // det_m
+        if isinstance(site, system.Site):
+            result = ta.dot(det_x_inv_m_part, site.tag) // det_m
+        elif isinstance(site, system.SiteArray):
+            result = np.dot(det_x_inv_m_part, site.tags.transpose()) // det_m
+        else:
+            raise TypeError("'site' must be a Site or a SiteArray")
+
         return -result if self.is_reversed else result
 
     def act(self, element, a, b=None):
-        element = ta.array(element)
-        if element.dtype is not int:
+        is_site = isinstance(a, system.Site)
+        # Tinyarray for small arrays (single site) else numpy
+        array_mod = ta if is_site else np
+        element = array_mod.array(element)
+        if not np.issubdtype(element.dtype, np.integer):
             raise ValueError("group element must be a tuple of integers")
+        if (len(element.shape) == 2 and is_site):
+            raise ValueError("must provide a single group element when "
+                             "acting on single sites.")
+        if (len(element.shape) == 1 and not is_site):
+            raise ValueError("must provide a sequence of group elements "
+                             "when acting on site arrays.")
         m_part = self._get_site_family_data(a.family)[0]
         try:
-            delta = ta.dot(m_part, element)
+            delta = array_mod.dot(m_part, element)
         except ValueError:
             msg = 'Expecting a {0}-tuple group element, but got `{1}` instead.'
             raise ValueError(msg.format(self.num_directions, element))
         if self.is_reversed:
             delta = -delta
         if b is None:
-            return builder.Site(a.family, a.tag + delta, True)
+            if is_site:
+                return system.Site(a.family, a.tag + delta, True)
+            else:
+                return system.SiteArray(a.family, a.tags + delta.transpose())
         elif b.family == a.family:
-            return (builder.Site(a.family, a.tag + delta, True),
-                    builder.Site(b.family, b.tag + delta, True))
+            if is_site:
+                return (system.Site(a.family, a.tag + delta, True),
+                        system.Site(b.family, b.tag + delta, True))
+            else:
+                return (system.SiteArray(a.family, a.tags + delta.transpose()),
+                        system.SiteArray(b.family, b.tags + delta.transpose()))
         else:
             m_part = self._get_site_family_data(b.family)[0]
             try:
@@ -717,8 +750,12 @@ class TranslationalSymmetry(builder.Symmetry):
                 raise ValueError(msg.format(self.num_directions, element))
             if self.is_reversed:
                 delta2 = -delta2
-            return (builder.Site(a.family, a.tag + delta, True),
-                    builder.Site(b.family, b.tag + delta2, True))
+            if is_site:
+                return (system.Site(a.family, a.tag + delta, True),
+                        system.Site(b.family, b.tag + delta2, True))
+            else:
+                return (system.SiteArray(a.family, a.tags + delta.transpose()),
+                        system.SiteArray(b.family, b.tags + delta2.transpose()))
 
     def reversed(self):
         """Return a reversed copy of the symmetry.

kwant/linalg/__init__.py

@@ -10,7 +10,10 @@ __all__ = ['lapack']
 from . import lapack
 
 # Merge the public interface of the other submodules.
-for module in ['decomp_lu', 'decomp_ev', 'decomp_schur']:
-    exec('from . import {0}'.format(module))
-    exec('from .{0} import *'.format(module))
-    exec('__all__.extend({0}.__all__)'.format(module))
+from .decomp_lu import *
+from .decomp_schur import *
+from .decomp_ev import *
+
+__all__.extend([decomp_lu.__all__,
+                decomp_ev.__all__,
+                decomp_schur.__all__])

kwant/linalg/tests/test_mumps.py

@@ -65,7 +65,7 @@ def test_schur_complement_with_dense():
 def test_error_minus_9(r=10):
     """Test if MUMPSError -9 is properly caught by increasing memory"""
 
-    graphene = honeycomb()
+    graphene = honeycomb(norbs=1)
     a, b = graphene.sublattices
 
     def circle(pos):

kwant/operator.pyx

@@ -13,6 +13,7 @@ import cython
 from operator import itemgetter
 import functools as ft
 import collections
+import warnings
 
 import numpy as np
 import tinyarray as ta
@@ -21,10 +22,13 @@ from scipy.sparse import coo_matrix
 from libc cimport math
 
 from .graph.core cimport EdgeIterator
+from .graph.core import DisabledFeatureError, NodeDoesNotExistError
 from .graph.defs cimport gint
 from .graph.defs import gint_dtype
-from .system import InfiniteSystem
-from ._common import UserCodeError, get_parameters, deprecate_args
+from .system import is_infinite, Site
+from ._common import (
+    UserCodeError, KwantDeprecationWarning, get_parameters, deprecate_args
+)
 
 
 ################ Generic Utility functions
@@ -148,7 +152,7 @@ def _get_all_orbs(gint[:, :] where, gint[:, :] site_ranges):
 
 def _get_tot_norbs(syst):
     cdef gint _unused, tot_norbs
-    is_infinite_system = isinstance(syst, InfiniteSystem)
+    is_infinite_system = is_infinite(syst)
     n_sites = syst.cell_size if is_infinite_system else syst.graph.num_nodes
     _get_orbs(np.asarray(syst.site_ranges, dtype=gint_dtype),
               n_sites, &tot_norbs, &_unused)
@@ -159,34 +163,41 @@ def _normalize_site_where(syst, where):
     """Normalize the format of `where` when `where` contains sites.
 
     If `where` is None, then all sites in the system are returned.
-    If it is a general iterator then it is expanded into an array. If `syst`
-    is a finalized Builder then `where` should contain `Site` objects,
+    If it is a general sequence then it is expanded into an array. If `syst`
+    is a finalized Builder then `where` may contain `Site` objects,
     otherwise it should contain integers.
     """
     if where is None:
-        size = (syst.cell_size
-                if isinstance(syst, InfiniteSystem) else syst.graph.num_nodes)
-        _where = list(range(size))
+        if is_infinite(syst):
+            where = list(range(syst.cell_size))
+        else:
+            where = list(range(syst.graph.num_nodes))
     elif callable(where):
         try:
-            _where = [syst.id_by_site[a] for a in filter(where, syst.sites)]
+            where = [syst.id_by_site[s] for s in filter(where, syst.sites)]
         except AttributeError:
-            _where = list(filter(where, range(syst.graph.num_nodes)))
+            if is_infinite(syst):
+                where = [s for s in range(syst.cell_size) if where(s)]
+            else:
+                where = [s for s in range(syst.graph.num_nodes) if where(s)]
     else:
-        try:
-            _where = list(syst.id_by_site[s] for s in where)
-        except AttributeError:
-            _where = list(where)
-            if any(w < 0 or w >= syst.graph.num_nodes for w in _where):
-                raise ValueError('`where` contains sites that are not in the '
-                                 'system.')
+        if isinstance(where[0], Site):
+            try:
+                where = [syst.id_by_site[s] for s in where]
+            except AttributeError:
+                raise TypeError("'where' contains Sites, but the system is not "
+                                "a finalized Builder.")
 
-    if isinstance(syst, InfiniteSystem):
-        if any(w >= syst.cell_size for w in _where):
-            raise ValueError('Only sites in the fundamental domain may be '
-                             'specified using `where`.')
+    where = np.asarray(where, dtype=gint_dtype).reshape(-1, 1)
 
-    return np.asarray(_where, dtype=gint_dtype).reshape(len(_where), 1)
+    if is_infinite(syst) and np.any(where >= syst.cell_size):
+        raise ValueError('Only sites in the fundamental domain may be '
+                         'specified using `where`.')
+    if np.any(np.logical_or(where < 0, where >= syst.graph.num_nodes)):
+        raise ValueError('`where` contains sites that are not in the '
+                         'system.')
+
+    return where
 
 
 def _normalize_hopping_where(syst, where):
@@ -194,42 +205,50 @@ def _normalize_hopping_where(syst, where):
 
     If `where` is None, then all hoppings in the system are returned.
     If it is a general iterator then it is expanded into an array. If `syst` is
-    a finalized Builder then `where` should contain pairs of `Site` objects,
+    a finalized Builder then `where` may contain pairs of `Site` objects,
     otherwise it should contain pairs of integers.
     """
     if where is None:
         # we cannot extract the hoppings in the same order as they are in the
         # graph while simultaneously excluding all inter-cell hoppings
-        if isinstance(syst, InfiniteSystem):
+        if is_infinite(syst):
             raise ValueError('`where` must be provided when calculating '
                              'current in an InfiniteSystem.')
-        _where = list(syst.graph)
+        where = list(syst.graph)
     elif callable(where):
         if hasattr(syst, "sites"):
-            def idx_where(hop):
+            def idxwhere(hop):
                 a, b = hop
                 return where(syst.sites[a], syst.sites[b])
-            _where = list(filter(idx_where, syst.graph))
+            where = list(filter(idxwhere, syst.graph))
         else:
-            _where = list(filter(lambda h: where(*h), syst.graph))
+            where = list(filter(lambda h: where(*h), syst.graph))
     else:
+        if isinstance(where[0][0], Site):
+            try:
+                where = list((syst.id_by_site[a], syst.id_by_site[b])
+                               for a, b in where)
+            except AttributeError:
+                raise TypeError("'where' contains Sites, but the system is not "
+                                "a finalized Builder.")
+        # NOTE: if we ever have operators that contain elements that are
+        #       not in the system graph, then we should modify this check
         try:
-            _where = list((syst.id_by_site[a], syst.id_by_site[b])
-                           for a, b in where)
-        except AttributeError:
-            _where = list(where)
-            # NOTE: if we ever have operators that contain elements that are
-            #       not in the system graph, then we should modify this check
+            error = ValueError('`where` contains hoppings that are not '
+                               'in the system.')
             if any(not syst.graph.has_edge(*w) for w in where):
-                raise ValueError('`where` contains hoppings that are not in the '
-                                 'system.')
+                raise error
+        # If where contains: negative integers, or integers > # of sites
+        except (NodeDoesNotExistError, DisabledFeatureError):
+            raise error
+
+    where = np.asarray(where, dtype=gint_dtype)
 
-    if isinstance(syst, InfiniteSystem):
-        if any(a > syst.cell_size or b > syst.cell_size for a, b in _where):
-            raise ValueError('Only intra-cell hoppings may be specified '
-                             'using `where`.')
+    if is_infinite(syst) and np.any(where > syst.cell_size):
+        raise ValueError('Only intra-cell hoppings may be specified '
+                         'using `where`.')
 
-    return np.asarray(_where, dtype=gint_dtype)
+    return where
 
 
 ## These two classes are here to avoid using closures, as these will
@@ -683,7 +702,11 @@ cdef class _LocalOperator:
             return mat
 
         offsets, norbs = _get_all_orbs(self.where, self._site_ranges)
-        return  BlockSparseMatrix(self.where, offsets, norbs, get_ham)
+        # TODO: update operators to use 'hamiltonian_term' rather than
+        #       'hamiltonian'.
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=KwantDeprecationWarning)
+            return  BlockSparseMatrix(self.where, offsets, norbs, get_ham)
 
     def __getstate__(self):
         return (
@@ -716,10 +739,10 @@ cdef class Density(_LocalOperator):
         maps from site families to square matrices. If a function is given it
         must take the same arguments as the onsite Hamiltonian functions of the
         system.
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
         then this should be a sequence of integers. If a function is provided,
-        it should take a single `int` or `~kwant.builder.Site` (if ``syst`` is
+        it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
         a finalized builder) and return True or False.  If not provided, the
         operator will be calculated over all sites in the system.
     check_hermiticity: bool
@@ -865,11 +888,11 @@ cdef class Current(_LocalOperator):
         matrices (scalars are allowed if the site family has 1 orbital per
         site). If a function is given it must take the same arguments as the
         onsite Hamiltonian functions of the system.
-    where : sequence of pairs of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of pairs of `int` or `~kwant.system.Site`, or callable, optional
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
         then this should be a sequence of pairs of integers. If a function is
         provided, it should take a pair of integers or a pair of
-        `~kwant.builder.Site` (if ``syst`` is a finalized builder) and return
+        `~kwant.system.Site` (if ``syst`` is a finalized builder) and return
         True or False.  If not provided, the operator will be calculated over
         all hoppings in the system.
     check_hermiticity : bool
@@ -997,10 +1020,10 @@ cdef class Source(_LocalOperator):
         matrices (scalars are allowed if the site family has 1 orbital per
         site). If a function is given it must take the same arguments as the
         onsite Hamiltonian functions of the system.
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
         then this should be a sequence of integers. If a function is provided,
-        it should take a single `int` or `~kwant.builder.Site` (if ``syst`` is
+        it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
         a finalized builder) and return True or False.  If not provided, the
         operator will be calculated over all sites in the system.
     check_hermiticity : bool

kwant/physics/__init__.py

@@ -8,8 +8,14 @@
 """Physics-related algorithms"""
 
 # Merge the public interface of all submodules.
-__all__ = []
-for module in ['leads', 'dispersion', 'noise', 'symmetry', 'gauge']:
-    exec('from . import {0}'.format(module))
-    exec('from .{0} import *'.format(module))
-    exec('__all__.extend({0}.__all__)'.format(module))
+from .leads import *
+from .dispersion import *
+from .noise import *
+from .symmetry import *
+from .gauge import *
+
+__all__ = [leads.__all__
+           + dispersion.__all__
+           + noise.__all__
+           + symmetry.__all__
+           + gauge.__all__]

kwant/physics/gauge.py

-# Copyright 2011-2018 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -303,7 +303,7 @@ def _loops_in_infinite(syst):
     extended_sites : callable : int -> Site
         Given a site index in the extended system consisting of
         two unit cells, returns the associated high-level
-        `kwant.builder.Site`.
+        `kwant.system.Site`.
     """
     assert isinstance(syst, system.InfiniteSystem)
     _check_infinite_syst(syst)
@@ -375,7 +375,7 @@ def _loops_in_composite(syst):
         -1 if the site is part of the reduced scattering region (see notes).
     extended_sites : callable : int -> Site
         Given a site index in the extended scattering region (see notes),
-        returns the associated high-level `kwant.builder.Site`.
+        returns the associated high-level `kwant.system.Site`.
 
     Notes
     -----
@@ -401,7 +401,7 @@ def _loops_in_composite(syst):
     # Get distance matrix for the extended scattering region,
     # a function that maps sites to their lead patches (-1 for sites
     # in the reduced scattering region), and a function that maps sites
-    # to high-level 'kwant.builder.Site' objects.
+    # to high-level 'kwant.system.Site' objects.
     distance_matrix, which_patch, extended_sites =\
         _extended_scattering_region(syst)
 
@@ -439,7 +439,7 @@ def _extended_scattering_region(syst):
         -1 if the site is part of the reduced scattering region.
     extended_sites : callable : int -> Site
         Given a site index in the extended scattering region, returns
-        the associated high-level `kwant.builder.Site`.
+        the associated high-level `kwant.system.Site`.
 
     Notes
     -----
@@ -1027,7 +1027,7 @@ class magnetic_gauge:
             The first callable computes the Peierls phase in the scattering
             region and the remaining callables compute the Peierls phases
             in each of the leads. Each callable takes a pair of
-            `~kwant.builder.Site` (a hopping) and returns a unit complex
+            `~kwant.system.Site` (a hopping) and returns a unit complex
             number (Peierls phase) that multiplies that hopping.
         """
         return self._peierls(syst_field, *lead_fields, tol=tol, average=False)

kwant/physics/leads.py

@@ -173,11 +173,11 @@ class StabilizedModes:
         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.
+    sqrt_hop : numpy array
+        Part of the SVD of `h_hop`.
     """
 
-    def __init__(self, vecs, vecslmbdainv, nmodes, sqrt_hop=None):
+    def __init__(self, vecs, vecslmbdainv, nmodes, sqrt_hop):
         kwargs = locals()
         kwargs.pop('self')
         self.__dict__.update(kwargs)

kwant/physics/tests/test_dispersion.py

@@ -17,7 +17,7 @@ from math import pi, cos, sin
 
 def make_lead():
     syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
-    lat = kwant.lattice.square()
+    lat = kwant.lattice.square(norbs=1)
     syst[[lat(0, 0), lat(0, 1)]] = 3
     syst[lat(0, 1), lat(0, 0)] = -1
     syst[((lat(1, y), lat(0, y)) for y in range(2))] = -1
@@ -35,7 +35,7 @@ def test_band_energies(N=5):
 
 def test_same_as_lead():
     syst = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
-    lat = kwant.lattice.chain()
+    lat = kwant.lattice.chain(norbs=1)
     syst[lat(0)] = 0
     syst[lat(0), lat(1)] = complex(cos(0.2), sin(0.2))
 
@@ -49,7 +49,7 @@ def test_same_as_lead():
 
 def test_raise_nonhermitian():
     syst = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
-    lat = kwant.lattice.chain()
+    lat = kwant.lattice.chain(norbs=1)
     syst[lat(0)] = 1j
     syst[lat(0), lat(1)] = complex(cos(0.2), sin(0.2))
     syst = syst.finalized()
@@ -58,7 +58,7 @@ def test_raise_nonhermitian():
 
 def test_band_velocities():
     syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
-    lat = kwant.lattice.square()
+    lat = kwant.lattice.square(norbs=1)
     syst[lat(0, 0)] = 1
     syst[lat(0, 1)] = 3
     syst[lat(1, 0), lat(0, 0)] = -1
@@ -75,7 +75,7 @@ def test_band_velocities():
 
 def test_band_velocity_derivative():
     syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
-    lat = kwant.lattice.square()
+    lat = kwant.lattice.square(norbs=1)
     syst[lat(0, 0)] = 1
     syst[lat(0, 1)] = 3
     syst[lat(1, 0), lat(0, 0)] = -1

kwant/physics/tests/test_leads.py

@@ -267,7 +267,7 @@ def test_modes():
 
 def test_modes_bearded_ribbon():
     # Check if bearded graphene ribbons work.
-    lat = kwant.lattice.honeycomb()
+    lat = kwant.lattice.honeycomb(norbs=1)
     syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
     syst[lat.shape((lambda pos: -20 < pos[1] < 20),
                   (0, 0))] = 0.3
@@ -417,7 +417,7 @@ def test_zero_hopping():
 
 def make_clean_lead(W, E, t):
     syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
-    lat = kwant.lattice.square()
+    lat = kwant.lattice.square(norbs=1)
     syst[(lat(0, j) for j in range(W))] = E
     syst[lat.neighbors()] = -t
     return syst.finalized()

kwant/physics/tests/test_noise.py

@@ -14,7 +14,7 @@ from kwant.physics import two_terminal_shotnoise
 from kwant._common import ensure_rng
 
 n = 5
-chain = kwant.lattice.chain()
+chain = kwant.lattice.chain(norbs=n)
 
 def twoterminal_system():
     rng = ensure_rng(11)

kwant/plotter.py

 # -*- coding: utf-8 -*-
-# Copyright 2011-2018 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -402,7 +402,7 @@ def sys_leads_sites(sys, num_lead_cells=2):
     Returns
     -------
     sites : list of (site, lead_number, copy_number) tuples
-        A site is a `~kwant.builder.Site` instance if the system is not finalized,
+        A site is a `~kwant.system.Site` instance if the system is not finalized,
         and an integer otherwise.  For system sites `lead_number` is `None` and
         `copy_number` is `0`, for leads both are integers.
     lead_cells : list of slices
@@ -427,7 +427,7 @@ def sys_leads_sites(sys, num_lead_cells=2):
                               lead.builder.sites() for i in
                               range(num_lead_cells)))
             lead_cells.append(slice(start, len(sites)))
-    elif isinstance(syst, system.FiniteSystem):
+    elif system.is_finite(syst):
         sites = [(i, None, 0) for i in range(syst.graph.num_nodes)]
         for leadnr, lead in enumerate(syst.leads):
             start = len(sites)
@@ -528,7 +528,7 @@ def sys_leads_hoppings(sys, num_lead_cells=2):
     Returns
     -------
     hoppings : list of (hopping, lead_number, copy_number) tuples
-        A site is a `~kwant.builder.Site` instance if the system is not finalized,
+        A site is a `~kwant.system.Site` instance if the system is not finalized,
         and an integer otherwise.  For system sites `lead_number` is `None` and
         `copy_number` is `0`, for leads both are integers.
     lead_cells : list of slices
@@ -735,17 +735,21 @@ def plot(sys, num_lead_cells=2, unit='nn',
         symbols specifications (only for kwant.system.FiniteSystem).
     site_size : number, function, array, or `None`
         Relative (linear) size of the site symbol.
+        An array may not be used when 'syst' is a kwant.Builder.
     site_color : ``matplotlib`` color description, function, array, or `None`
         A color used for plotting a site in the system. If a colormap is used,
         it should be a function returning single floats or a one-dimensional
         array of floats. By default sites are colored by their site family,
         using the current matplotlib color cycle.
+        An array of colors may not be used when 'syst' is a kwant.Builder.
     site_edgecolor : ``matplotlib`` color description, function, array, or `None`
         Color used for plotting the edges of the site symbols. Only
         valid matplotlib color descriptions are allowed (and no
         combination of floats and colormap as for site_color).
+        An array of colors may not be used when 'syst' is a kwant.Builder.
     site_lw : number, function, array, or `None`
         Linewidth of the site symbol edges.
+        An array may not be used when 'syst' is a kwant.Builder.
     hop_color : ``matplotlib`` color description or a function
         Same as `site_color`, but for hoppings.  A function is passed two sites
         in this case. (arrays are not allowed in this case).
@@ -808,7 +812,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
 
     - The meaning of "site" depends on whether the system to be plotted is a
       builder or a low level system.  For builders, a site is a
-      kwant.builder.Site object.  For low level systems, a site is an integer
+      kwant.system.Site object.  For low level systems, a site is an integer
       -- the site number.
 
     - color and symbol definitions may be tuples, but not lists or arrays.
@@ -911,7 +915,11 @@ def plot(sys, num_lead_cells=2, unit='nn',
             raise ValueError('Invalid value of unit argument.')
 
     # make all specs proper: either constant or lists/np.arrays:
-    def make_proper_site_spec(spec, fancy_indexing=False):
+    def make_proper_site_spec(spec_name,  spec, fancy_indexing=False):
+        if _p.isarray(spec) and isinstance(syst, builder.Builder):
+            raise TypeError('{} cannot be an array when plotting'
+                            ' a Builder; use a function instead.'
+                            .format(spec_name))
         if callable(spec):
             spec = [spec(i[0]) for i in sites if i[1] is None]
         if (fancy_indexing and _p.isarray(spec)
@@ -933,7 +941,8 @@ def plot(sys, num_lead_cells=2, unit='nn',
                 spec = np.asarray(spec, dtype='object')
         return spec
 
-    site_symbol = make_proper_site_spec(site_symbol)
+
+    site_symbol = make_proper_site_spec('site_symbol', site_symbol)
     if site_symbol is None: site_symbol = defaults['site_symbol'][dim]
     # separate different symbols (not done in 3D, the separation
     # would mess up sorting)
@@ -952,7 +961,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
 
     if site_color is None:
         cycle = _color_cycle()
-        if isinstance(syst, (builder.FiniteSystem, builder.InfiniteSystem)):
+        if builder.is_system(syst):
             # Skipping the leads for brevity.
             families = sorted({site.family for site in syst.sites})
             color_mapping = dict(zip(families, cycle))
@@ -967,10 +976,10 @@ def plot(sys, num_lead_cells=2, unit='nn',
             # Unknown finalized system, no sites access.
             site_color = defaults['site_color'][dim]
 
-    site_size = make_proper_site_spec(site_size, fancy_indexing)
-    site_color = make_proper_site_spec(site_color, fancy_indexing)
-    site_edgecolor = make_proper_site_spec(site_edgecolor, fancy_indexing)
-    site_lw = make_proper_site_spec(site_lw, fancy_indexing)
+    site_size = make_proper_site_spec('site_size', site_size, fancy_indexing)
+    site_color = make_proper_site_spec('site_color', site_color, fancy_indexing)
+    site_edgecolor = make_proper_site_spec('site_edgecolor', site_edgecolor, fancy_indexing)
+    site_lw = make_proper_site_spec('site_lw', site_lw, fancy_indexing)
 
     hop_color = make_proper_hop_spec(hop_color)
     hop_lw = make_proper_hop_spec(hop_lw)
@@ -1282,7 +1291,7 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None,
     if callable(value):
         value = [value(site[0]) for site in sites]
     else:
-        if not isinstance(syst, system.FiniteSystem):
+        if not system.is_finite(syst):
             raise ValueError('List of values is only allowed as input '
                              'for finalized systems.')
     value = np.array(value)
@@ -1398,7 +1407,7 @@ def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
                            "for bands()")
 
     syst = sys  # for naming consistency inside function bodies
-    _common.ensure_isinstance(syst, system.InfiniteSystem)
+    _common.ensure_isinstance(syst, (system.InfiniteSystem, system.InfiniteVectorizedSystem))
 
     momenta = np.array(momenta)
     if momenta.ndim != 1:
@@ -1474,7 +1483,7 @@ def spectrum(syst, x, y=None, params=None, mask=None, file=None,
     if y is not None and not _p.has3d:
         raise RuntimeError("Installed matplotlib does not support 3d plotting")
 
-    if isinstance(syst, system.FiniteSystem):
+    if system.is_finite(syst):
         def ham(**kwargs):
             return syst.hamiltonian_submatrix(params=kwargs, sparse=False)
     elif callable(syst):
@@ -1742,7 +1751,7 @@ def interpolate_current(syst, current, relwidth=None, abswidth=None, n=9):
         the extents of `field`: ((x0, x1), (y0, y1), ...)
 
     """
-    if not isinstance(syst, builder.FiniteSystem):
+    if not builder.is_finite_system(syst):
         raise TypeError("The system needs to be finalized.")
 
     if len(current) != syst.graph.num_edges:
@@ -1835,7 +1844,7 @@ def interpolate_density(syst, density, relwidth=None, abswidth=None, n=9,
         the extents of ``field``: ((x0, x1), (y0, y1), ...)
 
     """
-    if not isinstance(syst, builder.FiniteSystem):
+    if not builder.is_finite_system(syst):
         raise TypeError("The system needs to be finalized.")
 
     if len(density) != len(syst.sites):

kwant/solvers/common.py

@@ -210,8 +210,7 @@ class SparseSolver(metaclass=abc.ABCMeta):
                 iface_orbs = np.r_[tuple(slice(offsets[i], offsets[i + 1])
                                         for i in interface)]
 
-                n_lead_orbs = (svd_v.shape[0] if svd_v is not None
-                               else u_out.shape[0])
+                n_lead_orbs = svd_v.shape[0]
                 if n_lead_orbs != len(iface_orbs):
                     msg = ('Lead {0} has hopping with dimensions '
                            'incompatible with its interface dimension.')
@@ -221,25 +220,17 @@ class SparseSolver(metaclass=abc.ABCMeta):
                 transf = sp.csc_matrix((np.ones(len(iface_orbs)), coords),
                                        shape=(iface_orbs.size, lhs.shape[0]))
 
-                if svd_v is not None:
-                    v_sp = sp.csc_matrix(svd_v.T.conj()) * transf
-                    vdaguout_sp = (transf.T *
-                                   sp.csc_matrix(np.dot(svd_v, u_out)))
-                    lead_mat = - ulinv_out
-                else:
-                    v_sp = transf
-                    vdaguout_sp = transf.T * sp.csc_matrix(u_out)
-                    lead_mat = - ulinv_out
+                v_sp = sp.csc_matrix(svd_v.T.conj()) * transf
+                vdaguout_sp = (transf.T *
+                               sp.csc_matrix(np.dot(svd_v, 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_v is not None:
-                        vdaguin_sp = transf.T * sp.csc_matrix(
-                            -np.dot(svd_v, u_in))
-                    else:
-                        vdaguin_sp = transf.T * sp.csc_matrix(-u_in)
+                    vdaguin_sp = transf.T * sp.csc_matrix(
+                        -np.dot(svd_v, u_in))
 
                     # defer formation of the real matrix until the proper
                     # system size is known

kwant/solvers/tests/_test_sparse.py

@@ -15,8 +15,8 @@ import kwant
 from kwant._common import ensure_rng
 
 n = 5
-chain = kwant.lattice.chain()
-sq = square = kwant.lattice.square()
+chain = kwant.lattice.chain(norbs=n)
+sq = square = kwant.lattice.square(norbs=n)
 
 
 class LeadWithOnlySelfEnergy:
@@ -111,6 +111,8 @@ def test_one_lead(smatrix):
 # Test that a system with one lead with no propagating modes has a
 # 0x0 S-matrix.
 def test_smatrix_shape(smatrix):
+    chain = kwant.lattice.chain(norbs=1)
+
     system = kwant.Builder()
     lead0 = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
     lead1 = kwant.Builder(kwant.TranslationalSymmetry((1,)))
@@ -264,6 +266,8 @@ def test_singular_graph_system(smatrix):
 # zero eigenvalues than the lead hopping matrix. Older version of the
 # sparse solver failed here.
 def test_tricky_singular_hopping(smatrix):
+    sq = kwant.lattice.square(norbs=1)
+
     system = kwant.Builder()
     lead = kwant.Builder(kwant.TranslationalSymmetry((4, 0)))
 
@@ -294,6 +298,7 @@ def test_tricky_singular_hopping(smatrix):
 # Test the consistency of transmission and conductance_matrix for a four-lead
 # system without time-reversal symmetry.
 def test_many_leads(*factories):
+    sq = kwant.lattice.square(norbs=1)
     E=2.1
     B=0.01
 
@@ -428,7 +433,7 @@ def test_selfenergy_reflection(greens_function, smatrix):
 
 def test_very_singular_leads(smatrix):
     syst = kwant.Builder()
-    chain = kwant.lattice.chain()
+    chain = kwant.lattice.chain(norbs=2)
     left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
     right_lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
     syst[chain(0)] = left_lead[chain(0)] = right_lead[chain(0)] = np.identity(2)
@@ -443,7 +448,7 @@ def test_very_singular_leads(smatrix):
 
 def test_ldos(ldos):
     syst = kwant.Builder()
-    chain = kwant.lattice.chain()
+    chain = kwant.lattice.chain(norbs=1)
     lead = kwant.Builder(kwant.TranslationalSymmetry(chain.vec((1,))))
     syst[chain(0)] = syst[chain(1)] = lead[chain(0)] = 0
     syst[chain(0), chain(1)] = lead[chain(0), chain(1)] = 1
@@ -512,6 +517,7 @@ def test_arg_passing(wave_function, ldos, smatrix):
     def hopping(site1, site2, a, b):
         return b - a
 
+    square = kwant.lattice.square(norbs=1)
     W = 3
     L = 4