diff --git a/AUTHORS.rst b/AUTHORS.rst
index 43f37b31613b1670f262012847c73ad3acb426fd..3a39755ad20e78b5eaa7daffb2ad1aeb0856bc74 100644
--- a/AUTHORS.rst
+++ b/AUTHORS.rst
@@ -15,6 +15,7 @@ The principal developers of Kwant are
Contributors to Kwant include
* Jörg Behrmann (FU Berlin)
+* Paul Clisson (CEA Grenoble)
* Mathieu Istas (CEA Grenoble)
* Daniel Jaschke (CEA Grenoble)
* Bas Nijholt (TU Delft)
@@ -25,6 +26,7 @@ Contributors to Kwant include
* Sebastian Rubbert (TU Delft)
* Rafał Skolasiński (TU Delft)
* Adrien Sorgniard (CEA Grenoble)
+* Dániel Varjas (TU Delft)
We thank Christoph Gohlke for the creation of installers for Microsoft Windows.
diff --git a/doc/source/code/figure/closed_system.py.diff b/doc/source/code/figure/closed_system.py.diff
index 7f1de4715780c3c04f7ccf668301027617fff545..21f979f983eae19bf11dd35eab3e7d3d02a594bf 100644
--- a/doc/source/code/figure/closed_system.py.diff
+++ b/doc/source/code/figure/closed_system.py.diff
@@ -1,4 +1,4 @@
-@@ -1,140 +1,157 @@
+@@ -1,144 +1,162 @@
# Tutorial 2.4.2. Closed systems
# ==============================
#
@@ -94,14 +94,18 @@
+ fig.savefig("closed_system_result." + extension, dpi=_defs.dpi)
#HIDDEN_END_yvri
+ def sorted_eigs(ev):
+ evals, evecs = ev
+ evals, evecs = map(np.array, zip(*sorted(zip(evals, evecs.transpose()))))
+ return evals, evecs.transpose()
#HIDDEN_BEGIN_wave
- def plot_wave_function(syst):
+ def plot_wave_function(syst, B=0.001):
+ size = (_defs.figwidth_in, _defs.figwidth_in)
+
# Calculate the wave functions in the system.
- ham_mat = syst.hamiltonian_submatrix(sparse=True)
- evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
+ ham_mat = syst.hamiltonian_submatrix(sparse=True, args=[B])
+ evals, evecs = sorted_eigs(sla.eigsh(ham_mat, k=20, which='SM'))
# Plot the probability density of the 10th eigenmode.
- kwant.plotter.map(syst, np.abs(evecs[:, 9])**2,
@@ -115,16 +119,16 @@
#HIDDEN_BEGIN_current
- def plot_current(syst):
+ def plot_current(syst, B=0.001):
+ size = (_defs.figwidth_in, _defs.figwidth_in)
+
# Calculate the wave functions in the system.
- ham_mat = syst.hamiltonian_submatrix(sparse=True)
- evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
+ ham_mat = syst.hamiltonian_submatrix(sparse=True, args=[B])
+ evals, evecs = sorted_eigs(sla.eigsh(ham_mat, k=20, which='SM'))
# Calculate and plot the local current of the 10th eigenmode.
J = kwant.operator.Current(syst)
- current = J(evecs[:, 9])
+ current = J(evecs[:, 9], args=[B])
- kwant.plotter.current(syst, current, colorbar=False)
+ for extension in ('pdf', 'png'):
+ kwant.plotter.current(
diff --git a/doc/source/code/figure/spin_orbit.py.diff b/doc/source/code/figure/spin_orbit.py.diff
index ba32e81221f498d1b7057502b566632fbdf202df..794c1b5b755e22ea1a175c40a9078bd3e5b3e12b 100644
--- a/doc/source/code/figure/spin_orbit.py.diff
+++ b/doc/source/code/figure/spin_orbit.py.diff
@@ -34,10 +34,10 @@
#HIDDEN_END_hwbt
- def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
+ def make_system(t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
# Start with an empty tight-binding system and a single square lattice.
# `a` is the lattice constant (by default set to 1 for simplicity).
- lat = kwant.lattice.square(a)
+ lat = kwant.lattice.square()
syst = kwant.Builder()
@@ -47,23 +47,23 @@
4 * t * sigma_0 + e_z * sigma_z
# hoppings in x-direction
syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
- -t * sigma_0 - 1j * alpha * sigma_y
+ -t * sigma_0 + 1j * alpha * sigma_y / 2
# hoppings in y-directions
syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
- -t * sigma_0 + 1j * alpha * sigma_x
+ -t * sigma_0 - 1j * alpha * sigma_x / 2
#HIDDEN_END_uxrm
#### Define the left lead. ####
- lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
#HIDDEN_BEGIN_yliu
lead[(lat(0, j) for j in range(W))] = 4 * t * sigma_0 + e_z * sigma_z
# hoppings in x-direction
lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
- -t * sigma_0 - 1j * alpha * sigma_y
+ -t * sigma_0 + 1j * alpha * sigma_y / 2
# hoppings in y-directions
lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
- -t * sigma_0 + 1j * alpha * sigma_x
+ -t * sigma_0 - 1j * alpha * sigma_x / 2
#HIDDEN_END_yliu
#### Attach the leads and return the finalized system. ####
diff --git a/doc/source/tutorial/faq.rst b/doc/source/tutorial/faq.rst
index 3833f4743bdd9450f051a426c3ad745ee9005fa1..1022bf03982fa419a95faba812013ee02f2771d4 100644
--- a/doc/source/tutorial/faq.rst
+++ b/doc/source/tutorial/faq.rst
@@ -434,7 +434,12 @@ Scattering states calculated using `~kwant.solvers.default.wave_function` are re
same order as the "incoming" modes of `~kwant.system.InfiniteSystem.modes`.
Kwant considers that the translational symmetry of a lead points "towards
infinity" (*not* towards the system) which means that the incoming modes are
-those that have *negative* velocities:
+those that have *negative* velocities.
+
+This means that for a lead attached on the left of a scattering region (with
+symmetry vector $(-1, 0)$, for example), the
+positive $k$ direction (when inspecting the lead's band structure) actually
+corresponds to the *negative* $x$ direction.
How does Kwant order components of an individual wavefunction?
diff --git a/doc/source/tutorial/first_steps.rst b/doc/source/tutorial/first_steps.rst
index 9a74435f0ef9fb3784e6b638893c1e7dc2545900..1ddce048c4da222856413471aeb452390ed3e7c7 100644
--- a/doc/source/tutorial/first_steps.rst
+++ b/doc/source/tutorial/first_steps.rst
@@ -315,12 +315,12 @@ subbands that increases with energy.
- Instead of plotting to the screen (which is standard)
`~kwant.plotter.plot` can also write to a file specified by the argument
- ``file``. For the plotting to the screen to work the module
- ``matplotlib.pyplot`` has to be imported. (An informative error message
- will remind you if you forget.) The reason for this is pretty technical:
- matplotlib's "backend" can only be chosen before ``matplotlib.pyplot`` has
- been imported. Would Kwant import that module by itself, it would deprive
- you of the possibility to choose a non-default backend later.
+ ``file``.
+
+ - Due to matplotlib's limitations, Kwant's plotting routines have the
+ side effect of fixing matplotlib's "backend". If you would like to choose
+ a different backend than the standard one, you must do so before asking
+ Kwant to plot anything.
.. rubric:: Footnotes
diff --git a/kwant/_common.py b/kwant/_common.py
index 3af56a12b2618f2cb7f7ba550e0adf13223a6ced..76e952cf6c0d2ffc437d42ca5c093a5825aac0e9 100644
--- a/kwant/_common.py
+++ b/kwant/_common.py
@@ -9,6 +9,8 @@
import numpy as np
import numbers
import inspect
+import warnings
+from contextlib import contextmanager
__all__ = ['KwantDeprecationWarning', 'UserCodeError']
@@ -87,6 +89,14 @@ def ensure_rng(rng=None):
"numpy.random.RandomState interface.")
+@contextmanager
+def reraise_warnings(level=3):
+ with warnings.catch_warnings(record=True) as caught_warnings:
+ yield
+ for warning in caught_warnings:
+ warnings.warn(warning.message, stacklevel=level)
+
+
def get_parameters(func):
"""Get the names of the parameters to 'func' and whether it takes kwargs.
diff --git a/kwant/builder.py b/kwant/builder.py
index d2af7e47dcdee9c5b641ef7ec78ff1f742eacb44..1418b40ae86678d91f64034dd58a7669725c24d1 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -20,7 +20,7 @@ from . import system, graph, KwantDeprecationWarning, UserCodeError
from .linalg import lll
from .operator import Density
from .physics import DiscreteSymmetry
-from ._common import ensure_isinstance, get_parameters
+from ._common import ensure_isinstance, get_parameters, reraise_warnings
__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
@@ -393,6 +393,8 @@ class NoSymmetry(Symmetry):
class HoppingKind(tuple):
"""A pattern for matching hoppings.
+ An alias exists for this common name: ``kwant.HoppingKind``.
+
A hopping ``(a, b)`` matches precisely when the site family of ``a`` equals
`family_a` and that of ``b`` equals `family_b` and ``(a.tag - b.tag)`` is
equal to `delta`. In other words, the matching hoppings have the form:
@@ -442,6 +444,17 @@ class HoppingKind(tuple):
else:
ensure_isinstance(family_b, SiteFamily)
family_b = family_b
+
+ try:
+ Site(family_b, family_a.normalize_tag(delta) - delta)
+ except Exception as e:
+ same_fams = family_b is family_a
+ msg = (str(family_a),
+ 'and {} are'.format(family_b) if not same_fams else ' is',
+ 'not compatible with delta={}'.format(delta),
+ )
+ raise ValueError(' '.join(msg)) from e
+
return tuple.__new__(cls, (delta, family_a, family_b))
def __call__(self, builder):
@@ -700,6 +713,8 @@ def _site_ranges(sites):
class Builder:
"""A tight binding system defined on a graph.
+ An alias exists for this common name: ``kwant.Builder``.
+
This is one of the central types in Kwant. It is used to construct tight
binding systems in a flexible way.
@@ -1514,8 +1529,9 @@ class Builder:
if hop_range > 1:
# Automatically increase the period, potentially warn the user.
new_lead = Builder(sym.subgroup((hop_range,)))
- new_lead.fill(lead_builder, lambda site: True,
- lead_builder.sites(), max_sites=float('inf'))
+ with reraise_warnings():
+ new_lead.fill(lead_builder, lambda site: True,
+ lead_builder.sites(), max_sites=float('inf'))
lead_builder = new_lead
sym = lead_builder.symmetry
H = lead_builder.H
@@ -1565,8 +1581,9 @@ class Builder:
# system (this one is guaranteed to contain a complete unit cell of the
# lead). After flood-fill we remove that domain.
start = {sym.act((max_dom + 1,), site) for site in H}
- all_added = self.fill(lead_builder, shape, start,
- max_sites=float('inf'))
+ with reraise_warnings():
+ all_added = self.fill(lead_builder, shape, start,
+ max_sites=float('inf'))
all_added = [site for site in all_added if site not in start]
del self[start]
@@ -1614,6 +1631,18 @@ class Builder:
raise ValueError('Currently, only builders without or with a 1D '
'translational symmetry can be finalized.')
+ # Protect novice users from confusing error messages if they
+ # forget to finalize their Builder.
+
+ @staticmethod
+ def _require_system(*args, **kwargs):
+ """You need a finalized system; Use Builder.finalized() first."""
+ raise TypeError('You need a finalized system; '
+ 'use Builder.finalized() first.')
+
+ hamiltonian = hamiltonian_submatrix = modes = selfenergy = \
+ inter_cell_hopping = cell_hamiltonian = precalculated = \
+ _require_system
################ Finalized systems
diff --git a/kwant/continuum/_common.py b/kwant/continuum/_common.py
index cb0effb8d12f634f1899cde20a8eb3a846e68b8d..ee9285bc99e83a87dddda983b4daedfd9a6c149e 100644
--- a/kwant/continuum/_common.py
+++ b/kwant/continuum/_common.py
@@ -21,6 +21,11 @@ from sympy.physics.matrices import msigma as _msigma
import warnings
+from .._common import reraise_warnings
+
+# TODO: remove when sympy correctly includes MutableDenseMatrix (lol).
+sympy_classes = set(sympy_classes) | {sympy.MutableDenseMatrix}
+
momentum_operators = sympy.symbols('k_x k_y k_z', commutative=False)
position_operators = sympy.symbols('x y z', commutative=False)
@@ -71,7 +76,8 @@ def lambdify(expr, locals=None):
[ 0. , 0. ]])
"""
- expr = sympify(expr, locals)
+ with reraise_warnings(level=4):
+ expr = sympify(expr, locals)
args = [s.name for s in expr.atoms(sympy.Symbol)]
args += [str(f.func) for f in expr.atoms(AppliedUndef, sympy.Function)]
@@ -149,8 +155,10 @@ def sympify(expr, locals=None):
# if ``expr`` is already a ``sympy`` object we may terminate a code path
if isinstance(expr, tuple(sympy_classes)):
if locals:
- warnings.warn('Input expression is already SymPy object: ' +
- '"locals" will not be used.', RuntimeWarning)
+ warnings.warn('Input expression is already SymPy object: '
+ '"locals" will not be used.',
+ RuntimeWarning,
+ stacklevel=2)
return expr
# if ``expr`` is not a "sympy" then we proceed with sympifying process
diff --git a/kwant/continuum/discretizer.py b/kwant/continuum/discretizer.py
index 31a0c29b1318aa02dcd82d45e76513f48c628a26..2560063cf1656a30bccea464966b2711a78a9586 100644
--- a/kwant/continuum/discretizer.py
+++ b/kwant/continuum/discretizer.py
@@ -6,6 +6,7 @@
# the file AUTHORS.rst at the top-level directory of this distribution and at
# http://kwant-project.org/authors.
+from keyword import iskeyword
from collections import defaultdict
import itertools
import warnings
@@ -22,6 +23,7 @@ from sympy.core.function import AppliedUndef
from .. import builder, lattice
from .. import KwantDeprecationWarning
+from .._common import reraise_warnings
from ._common import (sympify, gcd, position_operators, momentum_operators,
monomials)
@@ -67,7 +69,7 @@ class _DiscretizedBuilder(builder.Builder):
else:
a, b = key
assert a is site
- result.extend(["# Hopping in direction ",
+ result.extend(["# Hopping from ",
str(tuple(b.tag)),
":\n"])
result.append(val._source if callable(val) else repr(val))
@@ -190,10 +192,11 @@ def discretize_symbolic(hamiltonian, coords=None, *, locals=None):
The coordinates that have been discretized.
"""
- hamiltonian = sympify(hamiltonian, locals)
+ with reraise_warnings():
+ hamiltonian = sympify(hamiltonian, locals)
atoms_names = [s.name for s in hamiltonian.atoms(sympy.Symbol)]
- if any( s == 'a' for s in atoms_names):
+ if any(s == 'a' for s in atoms_names):
raise TypeError("'a' is a symbol used internally to represent "
"grid spacing; please use a different symbol.")
@@ -311,8 +314,9 @@ def build_discretized(tb_hamiltonian, coords, *, grid=None, locals=None,
raise ValueError("The argument 'coords' must be sorted.")
# run sympifcation on hamiltonian values
- for k, v in tb_hamiltonian.items():
- tb_hamiltonian[k] = sympify(v, locals)
+ with reraise_warnings():
+ for k, v in tb_hamiltonian.items():
+ tb_hamiltonian[k] = sympify(v, locals)
# generate grid if required, check constraints if provided
random_element = next(iter(tb_hamiltonian.values()))
@@ -644,7 +648,16 @@ def _builder_value(expr, coords, grid_spacing, onsite,
# constants and functions in the sympy input will be passed
# as arguments to the value function
arg_names = set.union({s.name for s in const_symbols},
- {str(k.func) for k in map_func_calls})
+ {str(k.func) for k in map_func_calls})
+
+ # check if all argument names are valid python identifiers
+ for name in arg_names:
+ if not (name.isidentifier() and not iskeyword(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(name))
+
arg_names = ', '.join(sorted(arg_names))
if (not arg_names) and (coords is None):
diff --git a/kwant/continuum/tests/test_discretizer.py b/kwant/continuum/tests/test_discretizer.py
index 7511e5a132daea58cd86365714ff4dd31dcae8dd..56ff8242da1d9932573167cc8158d04f5e18ae1f 100644
--- a/kwant/continuum/tests/test_discretizer.py
+++ b/kwant/continuum/tests/test_discretizer.py
@@ -570,3 +570,9 @@ def test_grid_input(ham, grid_offset, offset, norbs):
def test_grid_constraints(ham, coords, grid):
with pytest.raises(ValueError):
discretize(ham, coords, grid=grid)
+
+
+@pytest.mark.parametrize('name', ['1', '1a', '-a', '+a', 'while', 'for'])
+def test_check_symbol_names(name):
+ with pytest.raises(ValueError):
+ discretize(sympy.Symbol(name), 'x')
diff --git a/kwant/lattice.py b/kwant/lattice.py
index 8cf0d6a20065339ee389de2eb7cc2cc23cad074e..e050552c76e1a7503d4058d08c3c0f32cc660c02 100644
--- a/kwant/lattice.py
+++ b/kwant/lattice.py
@@ -504,8 +504,9 @@ class Monatomic(builder.SiteFamily, Polyatomic):
# point precision issues.
class TranslationalSymmetry(builder.Symmetry):
- """
- A translational symmetry defined in real space.
+ """A translational symmetry defined in real space.
+
+ An alias exists for this common name: ``kwant.TranslationalSymmetry``.
Group elements of this symmetry are integer tuples of appropriate length.
diff --git a/kwant/linalg/mumps.py b/kwant/linalg/mumps.py
index 397b4cc963bc9f9d593bce12fcd2c41904bc5912..465ef1aa359ab3c88adbf5a4a4ccf125c73094dd 100644
--- a/kwant/linalg/mumps.py
+++ b/kwant/linalg/mumps.py
@@ -299,7 +299,9 @@ class MUMPSContext:
if reuse_analysis:
if self.mumps_instance is None:
warnings.warn("Missing analysis although reuse_analysis=True. "
- "New analysis is performed.", RuntimeWarning)
+ "New analysis is performed.",
+ RuntimeWarning,
+ stacklevel=2)
self.analyze(a, ordering=ordering, overwrite_a=overwrite_a)
else:
dtype, row, col, data = _make_assembled_from_coo(a,
diff --git a/kwant/operator.pyx b/kwant/operator.pyx
index 8fc66e81c5e856b274e2e422f7466e38fbe61952..2787b28171504499c718df51cc7f0f4b106bbce5 100644
--- a/kwant/operator.pyx
+++ b/kwant/operator.pyx
@@ -428,7 +428,9 @@ cdef class _LocalOperator:
If True, checks that ``onsite``, as well as any relevant parts
of the Hamiltonian are hermitian.
sum : bool, default: False
- If True, then calling this operator will return a single scalar.
+ If True, then calling this operator will return a single scalar,
+ otherwise a vector will be returned (see
+ `~kwant.operator._LocalOperator.__call__` for details).
"""
cdef public int check_hermiticity, sum
@@ -468,14 +470,27 @@ cdef class _LocalOperator:
>>> A(phi, psi)
- to compute the matrix element :math:`\bra{φ} A \ket{ψ}`. Note that
- these quantities may be vectors (e.g. *local* charge or current
- density).
-
- For an operator :math:`Q_{iαβ}`, ``bra`` :math:`φ_α` and
- ``ket`` :math:`ψ_β` this computes :math:`q_i = ∑_{αβ} φ^*_α Q_{iαβ}
- ψ_β` if ``self.sum`` is False, otherwise computes :math:`q = ∑_{iαβ}
- φ^*_α Q_{iαβ} ψ_β`.
+ to compute the matrix element :math:`\bra{φ} A \ket{ψ}`.
+
+ If ``sum=True`` was provided when constructing the operator, then
+ a scalar is returned. If ``sum=False``, then a vector is returned.
+ The vector is defined over the sites of the system if the operator
+ is a `~kwant.operator.Density`, or over the hoppings if it is a
+ `~kwant.operator.Current` or `~kwant.operator.Source`. By default,
+ the returned vector is ordered in the same way as the sites
+ (for `~kwant.operator.Density`) or hoppings in the graph of the
+ system (for `~kwant.operator.Current` or `~kwant.operator.Density`).
+ If the keyword parameter ``where`` was provided when constructing
+ the operator, then the returned vector is instead defined only over
+ the sites or hoppings specified, and is ordered in the same way
+ as ``where``.
+
+ Alternatively stated, for an operator :math:`Q_{iαβ}`, ``bra``
+ :math:`φ_α` and ``ket`` :math:`ψ_β` this computes
+ :math:`q_i = ∑_{αβ} φ^*_α Q_{iαβ} ψ_β` if ``self.sum`` is False,
+ otherwise computes :math:`q = ∑_{iαβ} φ^*_α Q_{iαβ} ψ_β`. where
+ :math:`i` runs over all sites or hoppings, and
+ :math:`α` and :math:`β` run over all the degrees of freedom.
Parameters
----------
@@ -681,8 +696,8 @@ cdef class Density(_LocalOperator):
"""An operator for calculating general densities.
An instance of this class can be called like a function to evaluate the
- expectation value with a wavefunction. See the documentation of the
- ``__call__`` method for more details.
+ expectation value with a wavefunction. See
+ `~kwant.operator.Density.__call__` for details.
Parameters
----------
@@ -703,7 +718,9 @@ cdef class Density(_LocalOperator):
Hermitian, then an error will be raised when the operator is
evaluated.
sum : bool, default: False
- If True, then calling this operator will return a single scalar.
+ If True, then calling this operator will return a single scalar,
+ otherwise a vector will be returned (see
+ `~kwant.operator.Density.__call__` for details).
Notes
-----
@@ -818,8 +835,8 @@ cdef class Current(_LocalOperator):
r"""An operator for calculating general currents.
An instance of this class can be called like a function to evaluate the
- expectation value with a wavefunction. See the documentation of the
- ``__call__`` method for more details.
+ expectation value with a wavefunction. See
+ `~kwant.operator.Current.__call__` for details.
Parameters
----------
@@ -842,7 +859,9 @@ cdef class Current(_LocalOperator):
is not Hermitian, then an error will be raised when the
operator is evaluated.
sum : bool, default: False
- If True, then calling this operator will return a single scalar.
+ If True, then calling this operator will return a single scalar,
+ otherwise a vector will be returned (see
+ `~kwant.operator.Current.__call__` for details).
Notes
-----
@@ -945,8 +964,8 @@ cdef class Source(_LocalOperator):
"""An operator for calculating general sources.
An instance of this class can be called like a function to evaluate the
- expectation value with a wavefunction. See the documentation of the
- ``__call__`` method for more details.
+ expectation value with a wavefunction. See
+ `~kwant.operator.Source.__call__` for details.
Parameters
----------
@@ -968,7 +987,9 @@ cdef class Source(_LocalOperator):
Hermitian, then an error will be raised when the operator is
evaluated.
sum : bool, default: False
- If True, then calling this operator will return a single scalar.
+ If True, then calling this operator will return a single scalar,
+ otherwise a vector will be returned (see
+ `~kwant.operator.Source.__call__` for details).
Notes
-----
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index 4fcdd73ba925a8ff27cf0071cd4fc48ee4d4728c..b156f0292c3feafbc0a718d4ca80b73e97786533 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -121,6 +121,9 @@ class PropagatingModes:
velocity, `k` is momentum and `conserved_quantity` is the conservation
law eigenvalue.
+ In the above, the positive velocity and momentum directions are defined
+ with respect to the translational symmetry direction of the system.
+
The first dimension of `wave_functions` corresponds to the orbitals of all
the sites in a unit cell, the second one to the number of the mode. Each
mode is normalized to carry unit current. If several modes have the same
@@ -490,7 +493,7 @@ def phs_symmetrization(wfs, particle_hole):
----------
wfs : numpy array
A matrix of propagating wave functions at a TRIM that all have the same
- velocity. The wave functions form the columns of this matrix.
+ velocity. The orthonormal wave functions form the columns of this matrix.
particle_hole : numpy array
The matrix representation of the unitary part of the particle-hole
operator, expressed in the tight binding basis.
@@ -508,26 +511,44 @@ def phs_symmetrization(wfs, particle_hole):
"""Apply the particle-hole operator to an array. """
return particle_hole.dot(mat.conj())
- # P always squares to 1 or -1.
- P_squared = np.sign(particle_hole[0,:].dot(particle_hole[:,0].conj()))
- # np.sign returns the same data type as its argument. Make sure
- # that the comparison with integers is okay.
- assert P_squared in (-1, 1)
+ # Take P in the subspace of W = wfs: U = W^+ @ P @ W^*.
+ U = wfs.T.conj().dot(Pdot(wfs))
+ # Check that wfs are orthonormal and the space spanned
+ # by them is closed under ph, meaning U is unitary.
+ if not np.allclose(U.dot(U.T.conj()), np.eye(U.shape[0])):
+ raise ValueError('wfs are not orthonormal or not closed under particle_hole.')
+ P_squared = U.dot(U.conj())
+ if np.allclose(P_squared, np.eye(U.shape[0])):
+ P_squared = 1
+ elif np.allclose(P_squared, -np.eye(U.shape[0])):
+ P_squared = -1
+ else:
+ raise ValueError('particle_hole must square to +1 or -1. P_squared = {}'.format(P_squared))
if P_squared == 1:
- # Make particle hole eigenstates.
- # Phase factor ensures they are not numerically close.
- phases = np.diag([np.exp(1j*np.angle(wf.T.conj().dot(
- Pdot(wf)))*0.5) for wf in wfs.T])
- new_wfs = wfs.dot(phases) + Pdot(wfs.dot(phases))
- # Orthonormalize the modes using QR on the matrix of eigenstates of P.
- # So long as the matrix of coefficients R is purely real, any linear
- # combination of these modes remains an eigenstate of P. From the way
- # we construct eigenstates of P, the coefficients of R are real.
- new_wfs, r = la.qr(new_wfs, mode='economic', pivoting=True)[:2]
- if not np.allclose(r.imag, np.zeros(r.shape)): # skip coverage
- raise RuntimeError("Numerical instability in finding particle-hole \
- symmetric modes.")
+ # Use the matrix square root method from
+ # Applied Mathematics and Computation 234 (2014) 380-384.
+ assert np.allclose(U, U.T)
+ # Schur decomposition guarantees that vecs are orthonormal.
+ vals, vecs = la.schur(U)
+ # U should be normal, so vals is diagonal.
+ assert np.allclose(np.diag(np.diag(vals)), vals)
+ vals = np.diag(vals)
+ # Need to take safe square root of U, the branch cut should not go
+ # through any eigenvalues. Otherwise the square root may not be symmetric.
+ # Find largest gap between eigenvalues
+ phases = np.sort(np.angle(vals))
+ dph = np.append(np.diff(phases), phases[0] + 2*np.pi - phases[-1])
+ i = np.argmax(dph)
+ shift = -np.pi - (phases[i] + dph[i]/2)
+ # Take matrix square root with branch cut in largest gap
+ vals = np.sqrt(vals * np.exp(1j * shift)) * np.exp(-0.5j * shift)
+ sqrtU = vecs.dot(np.diag(vals)).dot(vecs.T.conj())
+ # For symmetric U sqrt(U) is also symmetric.
+ assert np.allclose(sqrtU, sqrtU.T)
+ # We want a new basis W_new such that W_new^+ @ P @ W_new^* = 1.
+ # This is achieved by W_new = W @ sqrt(U).
+ new_wfs = wfs.dot(sqrtU)
# If P^2 = 1, there is no need to sort the modes further.
TRIM_sort = np.zeros((wfs.shape[1],), dtype=int)
else:
@@ -540,7 +561,7 @@ def phs_symmetrization(wfs, particle_hole):
# If there are only two modes in this subspace, they are orthogonal
# so we replace the second one with the P applied to the first one.
if N_modes == 2:
- wf = wfs[:,0]
+ wf = wfs[:, 0]
# Store psi_n and P psi_n.
new_wfs.append(wf)
new_wfs.append(Pdot(wf))
@@ -553,22 +574,22 @@ def phs_symmetrization(wfs, particle_hole):
for i in iterations:
# Take a mode psi_n from the basis - the first column
# of the matrix of remaining modes.
- wf = wfs[:,0]
+ wf = wfs[:, 0]
# Store psi_n and P psi_n.
new_wfs.append(wf)
P_wf = Pdot(wf)
new_wfs.append(P_wf)
# Remove psi_n and P psi_n from the basis matrix of modes.
# First remove psi_n.
- wfs = wfs[:,1:]
+ wfs = wfs[:, 1:]
# Now we project the remaining modes onto the orthogonal
- # complement of P psi_n. Projector:
- Projector = wfs.dot(wfs.T.conj()) - \
+ # complement of P psi_n. projector:
+ projector = wfs.dot(wfs.T.conj()) - \
np.outer(P_wf, P_wf.T.conj())
# After the projection, the mode matrix is rank deficient -
# the span of the column space has dimension one less than
# the number of columns.
- wfs = Projector.dot(wfs)
+ wfs = projector.dot(wfs)
wfs = la.qr(wfs, mode='economic', pivoting=True)[0]
# Remove the redundant column.
wfs = wfs[:, :-1]
@@ -577,13 +598,12 @@ def phs_symmetrization(wfs, particle_hole):
# the projection.
if i == iterations[-1]:
assert wfs.shape[1] == 2
- wf = wfs[:,0]
+ wf = wfs[:, 0]
# Store psi_n and P psi_n.
new_wfs.append(wf)
- P_wf = Pdot(wf)
- new_wfs.append(P_wf)
+ new_wfs.append(Pdot(wf))
assert np.allclose(wfs.T.conj().dot(wfs),
- np.eye(wfs.shape[1]))
+ np.eye(wfs.shape[1]))
new_wfs = np.hstack([col.reshape(len(col), 1)/npl.norm(col) for
col in new_wfs])
assert np.allclose(new_wfs[:, 1::2], Pdot(new_wfs[:, ::2]))
diff --git a/kwant/physics/tests/test_leads.py b/kwant/physics/tests/test_leads.py
index cee31e82436cc45390e524b9e7fe0ade08e695c6..bc41f671811bd6c2cdae99d48c7b82fa8ebeeb3d 100644
--- a/kwant/physics/tests/test_leads.py
+++ b/kwant/physics/tests/test_leads.py
@@ -638,15 +638,11 @@ def test_PHS_TRIM():
for nmodes in (1, 3, n//4, n//2, n):
# Random matrix of 'modes.' Take part of a unitary
# matrix to ensure that the modes form a basis.
- modes = (rng.random_sample((n, n))
- + 1j*rng.random_sample((n, n)))
- modes = la.expm(1j*(modes
- + modes.T.conj()))[:n, :nmodes]
+ modes = kwant.rmt.circular(n, 'A', rng=rng)[:, :nmodes]
# Ensure modes are particle-hole symmetric and
- # normalized
+ # orthonormal
modes = modes + p_mat.dot(modes.conj())
- modes = np.array([col/np.linalg.norm(col) for col
- in modes.T]).T
+ modes = la.qr(modes, mode='economic')[0]
# Mix the modes with a random unitary transformation
U = kwant.rmt.circular(nmodes, 'A', rng=rng)
modes = modes.dot(U)
@@ -666,18 +662,14 @@ def test_PHS_TRIM():
for nmodes in (2, 4, n//2, n):
# Random matrix of 'modes.' Take part of a unitary
# matrix to ensure that the modes form a basis.
- modes = rng.rand(n, n) + 1j * rng.rand(n, n)
- modes = la.expm(1j*(modes +
- modes.T.conj()))[:n, :nmodes]
+ modes = kwant.rmt.circular(n, 'A', rng=rng)[:, :nmodes]
# Ensure modes are particle-hole symmetric and
# orthonormal.
modes[:, nmodes//2:] = \
p_mat.dot(modes[:, :nmodes//2].conj())
modes = la.qr(modes, mode='economic')[0]
# Mix the modes with a random unitary transformation
- U = (rng.random_sample((nmodes, nmodes))
- + 1j*rng.random_sample((nmodes, nmodes)))
- U = la.expm(1j*(U + U.T.conj()))
+ U = kwant.rmt.circular(nmodes, 'A', rng=rng)
modes = modes.dot(U)
# Make the modes PHS symmetric using the method for a
# TRIM.
@@ -690,7 +682,30 @@ def test_PHS_TRIM():
np.eye(phs_modes.shape[1]),
err_msg='Modes are not orthonormal,'
' TRIM PHS in ' + sym)
-
+ # Test the off-diagonal case when p_mat = sigma_x
+ p_mat = np.array([[0, 1], [1, 0]])
+ p_mat = np.kron(p_mat, np.identity(n // len(p_mat)))
+ for nmodes in (1, 3, n//4, n//2):
+ if nmodes > n//2:
+ continue
+ # Random matrix of 'modes.' Take part of a unitary
+ # matrix to ensure that the modes form a basis, all modes
+ # are only in half of the space
+ modes = kwant.rmt.circular(n//2, 'A', rng=rng)[:, :nmodes]
+ modes = np.vstack((modes, np.zeros((n//2, nmodes))))
+ # Add an orthogonal set of vectors that are ph images
+ modes = np.hstack((modes, p_mat.dot(modes.conj())))
+ # Make the modes PHS symmetric using the method for a
+ # TRIM.
+ phs_modes = leads.phs_symmetrization(modes, p_mat)[0]
+ assert_almost_equal(phs_modes,
+ p_mat.dot(phs_modes.conj()),
+ err_msg='PHS broken at a TRIM in '
+ 'off-diagonal test')
+ assert_almost_equal(phs_modes.T.conj().dot(phs_modes),
+ np.eye(phs_modes.shape[1]),
+ err_msg='Modes are not orthonormal,'
+ 'off-diagonal test')
def random_onsite_hop(n, rng=0):
diff --git a/kwant/plotter.py b/kwant/plotter.py
index 92998f772ffee6d1ff9a5fcf9f1a882326157eab..b1a2d537e8325bb0be573c90f7d27aa4f84db600 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -16,6 +16,7 @@ system in two or three dimensions.
"""
from collections import defaultdict
+import sys
import itertools
import functools
import warnings
@@ -35,9 +36,8 @@ try:
import matplotlib.cm
from matplotlib.figure import Figure
from matplotlib import collections
- from matplotlib.backends.backend_agg import FigureCanvasAgg
from . import _colormaps
- mpl_enabled = True
+ mpl_available = True
try:
from mpl_toolkits import mplot3d
has3d = True
@@ -47,7 +47,7 @@ try:
except ImportError:
warnings.warn("matplotlib is not available, only iterator-providing "
"functions will work.", RuntimeWarning)
- mpl_enabled = False
+ mpl_available = False
from . import system, builder, _common
@@ -77,7 +77,7 @@ def _sample_array(array, n_samples, rng=None):
return array[rng.choice(range(la), min(n_samples, la))]
-if mpl_enabled:
+if mpl_available:
class LineCollection(collections.LineCollection):
def __init__(self, segments, reflen=None, **kwargs):
super().__init__(segments, **kwargs)
@@ -285,8 +285,6 @@ if mpl_enabled:
self.set_paths(paths[indx])
if len(self.orig_transforms) > 1:
- self.transforms = np.resize(self.orig_transforms,
- (vs.shape[1],))
self.transforms = self.transforms[indx]
lw_orig = self.linewidths_orig
@@ -613,10 +611,9 @@ def output_fig(fig, output_mode='auto', file=None, savefile_opts=None,
The output mode to be used. Can be one of the following:
'pyplot' : attach the figure to pyplot, with the same behavior as if
pyplot.plot was called to create this figure.
- 'ipython' : attach a `FigureCanvasAgg` to the figure and return it.
- 'return' : return the figure.
- 'file' : same as 'ipython', but also save the figure into a file.
- 'auto' : if fname is given, save to a file, else if pyplot
+ 'return' : attach a `FigureCanvasAgg` to the figure and return it.
+ 'file' : same as 'return', but also save the figure into a file.
+ 'auto' : if fname is given, save to a file, otherwise like pyplot
is imported, attach to pyplot, otherwise just return. See also the
notes below.
file : string or a file object
@@ -634,24 +631,25 @@ def output_fig(fig, output_mode='auto', file=None, savefile_opts=None,
matplotlib in that the `dpi` attribute of the figure is used by defaul
instead of the matplotlib config setting.
"""
- if not mpl_enabled:
+ if not mpl_available:
raise RuntimeError('matplotlib is not installed.')
+
+ # We import backends and pyplot only at the last possible moment (=now)
+ # because this has the side effect of selecting the matplotlib backend for
+ # good. Warn if backend has not been set yet. This check is the same as
+ # the one performed inside matplotlib.use.
+ if 'matplotlib.backends' not in sys.modules:
+ warnings.warn("Kwant's plotting functions have\nthe side effect of "
+ "selecting the matplotlib backend. To avoid this "
+ "warning,\nimport matplotlib.pyplot, "
+ "matplotlib.backends or call matplotlib.use().",
+ RuntimeWarning, stacklevel=3)
+
if output_mode == 'auto':
- if file is not None:
- output_mode = 'file'
- else:
- try:
- matplotlib.pyplot.get_backend()
- output_mode = 'pyplot'
- except AttributeError:
- output_mode = 'pyplot'
+ output_mode = 'pyplot' if file is None else 'file'
if output_mode == 'pyplot':
- try:
- fake_fig = matplotlib.pyplot.figure()
- except AttributeError:
- msg = ('matplotlib.pyplot is unavailable. Execute `import '
- 'matplotlib.pyplot` or use a different output mode.')
- raise RuntimeError(msg)
+ from matplotlib import pyplot
+ fake_fig = pyplot.figure()
fake_fig.canvas.figure = fig
fig.canvas = fake_fig.canvas
for ax in fig.axes:
@@ -660,22 +658,19 @@ def output_fig(fig, output_mode='auto', file=None, savefile_opts=None,
except AttributeError:
pass
if show:
- matplotlib.pyplot.show()
- return fig
- elif output_mode == 'return':
- canvas = FigureCanvasAgg(fig)
- fig.canvas = canvas
- return fig
- elif output_mode == 'file':
- canvas = FigureCanvasAgg(fig)
- if savefile_opts is None:
- savefile_opts = ([], {})
- if 'dpi' not in savefile_opts[1]:
- savefile_opts[1]['dpi'] = fig.dpi
- canvas.print_figure(file, *savefile_opts[0], **savefile_opts[1])
- return fig
+ pyplot.show()
+ elif output_mode in ['return', 'file']:
+ from matplotlib.backends.backend_agg import FigureCanvasAgg
+ fig.canvas = FigureCanvasAgg(fig)
+ if output_mode == 'file':
+ if savefile_opts is None:
+ savefile_opts = ([], {})
+ if 'dpi' not in savefile_opts[1]:
+ savefile_opts[1]['dpi'] = fig.dpi
+ fig.canvas.print_figure(file, *savefile_opts[0], **savefile_opts[1])
else:
- assert False, 'Unknown output_mode'
+ raise ValueError('Unknown output_mode')
+ return fig
# Extracting necessary data from the system.
@@ -979,6 +974,8 @@ def plot(sys, num_lead_cells=2, unit='nn',
show=True, dpi=None, fig_size=None, ax=None):
"""Plot a system in 2 or 3 dimensions.
+ An alias exists for this common name: ``kwant.plot``.
+
Parameters
----------
sys : kwant.builder.Builder or kwant.system.FiniteSystem
@@ -1114,7 +1111,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
its aspect ratio.
"""
- if not mpl_enabled:
+ if not mpl_available:
raise RuntimeError("matplotlib was not found, but is required "
"for plot()")
@@ -1448,7 +1445,7 @@ def mask_interpolate(coords, values, a=None, method='nearest', oversampling=3):
if min_dist < 1e-6 * np.linalg.norm(cmax - cmin):
warnings.warn("Some sites have nearly coinciding positions, "
"interpolation may be confusing.",
- RuntimeWarning)
+ RuntimeWarning, stacklevel=2)
if a is None:
a = min_dist
@@ -1546,7 +1543,7 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None,
correspond to exactly one pixel.
"""
- if not mpl_enabled:
+ if not mpl_available:
raise RuntimeError("matplotlib was not found, but is required "
"for map()")
@@ -1567,7 +1564,8 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None,
raise ValueError('List of values is only allowed as input '
'for finalized systems.')
value = np.array(value)
- img, min, max = mask_interpolate(coords, value, a, method, oversampling)
+ with _common.reraise_warnings():
+ img, min, max = mask_interpolate(coords, value, a, method, oversampling)
border = 0.5 * (max - min) / (np.asarray(img.shape) - 1)
min -= border
max += border
@@ -1642,7 +1640,7 @@ def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
See `~kwant.physics.Bands` for the calculation of dispersion without plotting.
"""
- if not mpl_enabled:
+ if not mpl_available:
raise RuntimeError("matplotlib was not found, but is required "
"for bands()")
@@ -1717,7 +1715,7 @@ def spectrum(syst, x, y=None, params=None, mask=None, file=None,
A figure with the output if `ax` is not set, else None.
"""
- if not mpl_enabled:
+ if not mpl_available:
raise RuntimeError("matplotlib was not found, but is required "
"for plot_spectrum()")
if y is not None and not has3d:
@@ -2071,7 +2069,7 @@ def streamplot(field, box, cmap=None, bgcolor=None, linecolor='k',
fig : matplotlib figure
A figure with the output if `ax` is not set, else None.
"""
- if not mpl_enabled:
+ if not mpl_available:
raise RuntimeError("matplotlib was not found, but is required "
"for current()")
@@ -2173,8 +2171,9 @@ def current(syst, current, relwidth=0.05, **kwargs):
A figure with the output if `ax` is not set, else None.
"""
- return streamplot(*interpolate_current(syst, current, relwidth),
- **kwargs)
+ with _common.reraise_warnings(4):
+ return streamplot(*interpolate_current(syst, current, relwidth),
+ **kwargs)
# TODO (Anton): Fix plotting of parts of the system using color = np.nan.
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index f5f6197e9b791b2fb02c236ed0e26153fdc93f97..57c6fd35ec733165bd1aaf5ba04f0cfe6f598031 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -302,6 +302,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
"""
Compute the scattering matrix of a system.
+ An alias exists for this common name: ``kwant.smatrix``.
+
Parameters
----------
sys : `kwant.system.FiniteSystem`
@@ -394,6 +396,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
"""
Compute the retarded Green's function of the system between its leads.
+ An alias exists for this common name: ``kwant.greens_function``.
+
Parameters
----------
sys : `kwant.system.FiniteSystem`
@@ -489,6 +493,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
"""
Calculate the local density of states of a system at a given energy.
+ An alias exists for this common name: ``kwant.ldos``.
+
Parameters
----------
sys : `kwant.system.FiniteSystem`
@@ -553,6 +559,8 @@ class SparseSolver(metaclass=abc.ABCMeta):
Return a callable object for the computation of the wave function
inside the scattering region.
+ An alias exists for this common name: ``kwant.wave_function``.
+
Parameters
----------
sys : `kwant.system.FiniteSystem`
@@ -570,13 +578,23 @@ class SparseSolver(metaclass=abc.ABCMeta):
Notes
-----
-
The returned object can be itself called like a function. Given a lead
number, it returns a 2d NumPy array that contains the wave function
within the scattering region due to each incoming mode of the given
- lead. Index 0 is the mode number, index 1 is the orbital number. The
- modes appear in the same order as incoming modes in
- `kwant.physics.modes`.
+ lead. Index 0 is the mode number, index 1 is the orbital number.
+
+ The modes appear in the same order as the negative velocity modes in
+ `kwant.physics.PropagatingModes`. In Kwant's convention leads are attached
+ so that their translational symmetry points *away* from the scattering
+ region::
+
+ left lead SR right lead
+ /---------\ /---\ /---------\
+ ...-3-2-1-0-X-X-X-0-1-2-3-...
+
+ This means that incoming modes (coming from infinity towards the
+ scattering region) have *negative* velocity with respect to the
+ lead's symmetry direction.
Examples
--------
diff --git a/kwant/system.py b/kwant/system.py
index 62158d41819ef8f2b62a3e2e8c69f1ceed10f0f6..dfa28558d0efdb08a7ef8af788299e474e6b2759 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -224,6 +224,12 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
See documentation of `~kwant.physics.PropagatingModes` and
`~kwant.physics.StabilizedModes` for the return format details.
+
+ The wave functions of the returned modes are defined over the
+ *unit cell* of the system, which corresponds to the degrees of
+ freedom on the first ``cell_sites`` sites of the system
+ (recall that infinite systems store first the sites in the unit
+ cell, then connected sites in the neighboring unit cell).
"""
from . import physics # Putting this here avoids a circular import.
ham = self.cell_hamiltonian(args, params=params)
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index dee7ced6adb11640066ed49bf1b8b2289475c937..20879c44bed0156df810a16228f9eaf2c02b4573 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -1012,9 +1012,9 @@ def test_HoppingKind():
assert a.tag - b.tag == delta
# test hashability and equality
- hk = builder.HoppingKind((1, 0), g)
- hk2 = builder.HoppingKind((1, 0), g)
- hk3 = builder.HoppingKind((1, 0), g, h)
+ hk = builder.HoppingKind((1, 0, 0), g)
+ hk2 = builder.HoppingKind((1, 0, 0), g)
+ hk3 = builder.HoppingKind((1, 0, 0), g, h)
assert hk == hk2
assert hash(hk) == hash(hk2)
assert hk != hk3
@@ -1022,6 +1022,21 @@ def test_HoppingKind():
assert len({hk: 0, hk2:1, hk3: 2}) == 2
+def test_invalid_HoppingKind():
+ g = kwant.lattice.general(ta.identity(3))
+ h = kwant.lattice.general(np.identity(3)[:-1]) # 2D lattice in 3D
+
+ delta = (1, 0, 0)
+
+ # families have incompatible tags
+ with raises(ValueError):
+ builder.HoppingKind(delta, g, h)
+
+ # delta is incompatible with tags
+ with raises(ValueError):
+ builder.HoppingKind(delta, h)
+
+
def test_ModesLead_and_SelfEnergyLead():
lat = builder.SimpleSiteFamily()
hoppings = [builder.HoppingKind((1, 0), lat),
diff --git a/kwant/tests/test_plotter.py b/kwant/tests/test_plotter.py
index ca3b9b5129b37c65dda8683cf7cea42d4ddb6677..23b0d01a591be30145080734b8202aa6ea95fabe 100644
--- a/kwant/tests/test_plotter.py
+++ b/kwant/tests/test_plotter.py
@@ -15,14 +15,32 @@ from math import cos, sin
import scipy.integrate
import scipy.stats
import pytest
+import sys
import kwant
-from kwant import plotter
from kwant._common import ensure_rng
-if plotter.mpl_enabled:
+try:
from mpl_toolkits import mplot3d
+ import matplotlib
+
+ # This check is the same as the one performed inside matplotlib.use.
+ matplotlib_backend_chosen = 'matplotlib.backends' in sys.modules
+ # If the user did not already choose a backend, then choose
+ # the one with the least dependencies.
+ if not matplotlib_backend_chosen:
+ matplotlib.use('Agg')
+
from matplotlib import pyplot # pragma: no flakes
+except ImportError:
+ pass
+
+from kwant import plotter
+
+
+def test_matplotlib_backend_unset():
+ """Simply importing Kwant should not set the matplotlib backend."""
+ assert matplotlib_backend_chosen is False
def test_importable_without_matplotlib():
@@ -93,7 +111,7 @@ def syst_3d(W=3, r1=2, r2=4, a=1, t=1.0):
return syst
-@pytest.mark.skipif(not plotter.mpl_enabled, reason="No matplotlib available.")
+@pytest.mark.skipif(not plotter.mpl_available, reason="Matplotlib unavailable.")
def test_plot():
plot = plotter.plot
syst2d = syst_2d()
@@ -143,7 +161,7 @@ def bad_transform(pos):
x, y = pos
return x, y, 0
-@pytest.mark.skipif(not plotter.mpl_enabled, reason="No matplotlib available.")
+@pytest.mark.skipif(not plotter.mpl_available, reason="Matplotlib unavailable.")
def test_map():
syst = syst_2d()
with tempfile.TemporaryFile('w+b') as out:
@@ -179,7 +197,7 @@ def test_mask_interpolate():
coords, np.ones(2 * len(coords)))
-@pytest.mark.skipif(not plotter.mpl_enabled, reason="No matplotlib available.")
+@pytest.mark.skipif(not plotter.mpl_available, reason="Matplotlib unavailable.")
def test_bands():
syst = syst_2d().finalized().leads[0]
@@ -194,7 +212,7 @@ def test_bands():
plotter.bands(syst, ax=ax, file=out)
-@pytest.mark.skipif(not plotter.mpl_enabled, reason="No matplotlib available.")
+@pytest.mark.skipif(not plotter.mpl_available, reason="Matplotlib unavailable.")
def test_spectrum():
def ham_1d(a, b, c):
@@ -405,7 +423,7 @@ def test_current_interpolation():
assert scipy.stats.linregress(np.log(data))[2] < -0.8
-@pytest.mark.skipif(not plotter.mpl_enabled, reason="No matplotlib available.")
+@pytest.mark.skipif(not plotter.mpl_available, reason="Matplotlib unavailable.")
def test_current():
syst = syst_2d().finalized()
J = kwant.operator.Current(syst)
diff --git a/kwant/tests/test_wraparound.py b/kwant/tests/test_wraparound.py
index dc8ed3a927ea2248d4991d856b9f00310b812474..eb981d2798fa1b8fd42e36d26a0cd88039274b11 100644
--- a/kwant/tests/test_wraparound.py
+++ b/kwant/tests/test_wraparound.py
@@ -17,7 +17,7 @@ from kwant import plotter
from kwant.wraparound import wraparound, plot_2d_bands
from kwant._common import get_parameters
-if plotter.mpl_enabled:
+if plotter.mpl_available:
from mpl_toolkits import mplot3d # pragma: no flakes
from matplotlib import pyplot # pragma: no flakes
@@ -201,7 +201,7 @@ def test_symmetry():
assert np.all(orig == new)
-@pytest.mark.skipif(not plotter.mpl_enabled, reason="No matplotlib available.")
+@pytest.mark.skipif(not plotter.mpl_available, reason="Matplotlib unavailable.")
def test_plot_2d_bands():
chain = kwant.lattice.chain()
square = kwant.lattice.square()