diff --git a/doc/source/pre/whatsnew/0.2.rst b/doc/source/pre/whatsnew/0.2.rst
index 3a0e8e95a4a30a86168855fbd017be1279193cfe..2022810f63e0c00d3df442d91ff1e772f45ba0d8 100644
--- a/doc/source/pre/whatsnew/0.2.rst
+++ b/doc/source/pre/whatsnew/0.2.rst
@@ -42,7 +42,7 @@ description check `~kwant.plotter.plot` documentation.
The behavior of `~kwant.plotter.plot` with low level systems has changed.
Arguments of plot which are functions are given site numbers in place of
-`~kwant.builder.Site` objects when plotting a low level system. This
+`~kwant.system.Site` objects when plotting a low level system. This
provides an easy way to make the appearance of lines and symbols depend on
computation results.
diff --git a/doc/source/pre/whatsnew/1.0.rst b/doc/source/pre/whatsnew/1.0.rst
index 1eb4cf07efa8cda9b030497893b99b263afbabd7..473e414af1052f24ef271e38a2bd2b92ba34d361 100644
--- a/doc/source/pre/whatsnew/1.0.rst
+++ b/doc/source/pre/whatsnew/1.0.rst
@@ -39,7 +39,7 @@ now it suffices to write ::
syst[lat.neighbors()] = t
This is possible because `~kwant.builder.Builder` now accepts *functions* as
-keys in addition to `~kwant.builder.Site` objects and tuples of them
+keys in addition to `~kwant.system.Site` objects and tuples of them
(hoppings). These functions are expected to yield either sites or hoppings,
when given a builder instance as the sole argument. The use of such keys is to
implement sets of sites or hoppings that depend on what is already present in
diff --git a/doc/source/pre/whatsnew/1.5.rst b/doc/source/pre/whatsnew/1.5.rst
index 6204a4e02059822c4ab123ebec43b4e59afd5a5e..4a4d39755a638cfa2a749d19bcd2040b1ceeb5b0 100644
--- a/doc/source/pre/whatsnew/1.5.rst
+++ b/doc/source/pre/whatsnew/1.5.rst
@@ -3,6 +3,40 @@ What's new in Kwant 1.5
This article explains the user-visible changes in Kwant 1.5.0.
+
+Value functions may now be vectorized
+-------------------------------------
+It is now possible to define value functions that act on whole
+arrays of sites at a time.
+::
+
+ def onsite(sites):
+ r = sites.positions()
+ return np.exp(-np.linalg.norm(r, axis=1)**2)
+
+ def hopping(to_sites, from_sites):
+ dr = to_sites.positions() - from_sites.positions()
+ return 1 / np.linalg.norm(dr, axis=1)**2
+
+
+ lat = kwant.lattice.square(norbs=1)
+ syst = kwant.Builder(vectorize=True)
+ syst[(lat(i, j) for i in range(4) for j in range(10)] = onsite
+ syst[lat.neighbors()] = hopping
+ syst = syst.finalized()
+
+ syst.hamiltonian_submatrix()
+
+Notice that when creating the Builder we provide the ``vectorize=True`` flag,
+and that the value functions now take `~kwant.system.SiteArray` objects
+(which have a ``positions`` method, rather than a ``pos`` property as for
+`~kwant.system.Site`). We can now operate on the array of site positions
+using ``numpy``. Using vectorized value functions in this way can produce an
+order of magnitude speedup when evaluating system Hamiltonians (though this
+speedup may be masked by other parts of your computation e.g. calculating
+the scattering matrix).
+
+
Deprecation of leaving 'norbs' unset for site families
------------------------------------------------------
When constructing site families (e.g. lattices) it is now deprecated to
diff --git a/doc/source/reference/kwant.builder.rst b/doc/source/reference/kwant.builder.rst
index 327afca9e9f44c797378b918baa0799ef30af45f..6b470427c75c5ddb2594cdd8edacb705a36187ee 100644
--- a/doc/source/reference/kwant.builder.rst
+++ b/doc/source/reference/kwant.builder.rst
@@ -9,7 +9,6 @@ Types
:toctree: generated/
Builder
- Site
HoppingKind
SimpleSiteFamily
BuilderLead
@@ -17,13 +16,14 @@ Types
ModesLead
FiniteSystem
InfiniteSystem
+ FiniteVectorizedSystem
+ InfiniteVectorizedSystem
Abstract base classes
---------------------
.. autosummary::
:toctree: generated/
- SiteFamily
Symmetry
Lead
@@ -33,3 +33,11 @@ Functions
:toctree: generated/
add_peierls_phase
+
+Mixin Classes
+-------------
+.. autosummary::
+ :toctree: generated/
+
+ _FinalizedBuilderMixin
+ _VectorizedFinalizedBuilderMixin
diff --git a/doc/source/reference/kwant.system.rst b/doc/source/reference/kwant.system.rst
index 3b6f26274da393807c37c6c2826ce5518fe41c7d..2a9bd89c2529d65e5cc3cad63a144dc31e7a03f1 100644
--- a/doc/source/reference/kwant.system.rst
+++ b/doc/source/reference/kwant.system.rst
@@ -19,10 +19,32 @@ In practice, very often Kwant systems are finalized
`kwant.builder.FiniteSystem` or `kwant.builder.InfiniteSystem`) and offer some
additional attributes.
+Sites
+-----
+.. autosummary::
+ :toctree: generated/
+
+ Site
+ SiteArray
+ SiteFamily
+
+Systems
+-------
.. autosummary::
:toctree: generated/
System
+ VectorizedSystem
InfiniteSystem
+ InfiniteVectorizedSystem
FiniteSystem
+ FiniteVectorizedSystem
PrecalculatedLead
+
+Mixin Classes
+-------------
+.. autosummary::
+ :toctree: generated/
+
+ FiniteSystemMixin
+ InfiniteSystemMixin
diff --git a/doc/source/tutorial/first_steps.rst b/doc/source/tutorial/first_steps.rst
index 41ed21bb7fdfef0279b95f000db6d9199fa9ed63..b0760b6432e4d22003aa8bdfb9c18e5352999972 100644
--- a/doc/source/tutorial/first_steps.rst
+++ b/doc/source/tutorial/first_steps.rst
@@ -358,13 +358,13 @@ subbands that increases with energy.
**sites** as indices. Sites themselves have a certain type, and
belong to a **site family**. A site family is also used to convert
something that represents a site (like a tuple) into a
- proper `~kwant.builder.Site` object that can be used with
+ proper `~kwant.system.Site` object that can be used with
`~kwant.builder.Builder`.
In the above example, `lat` is the site family. ``lat(i, j)``
then translates the description of a lattice site in terms of two
integer indices (which is the natural way to do here) into
- a proper `~kwant.builder.Site` object.
+ a proper `~kwant.system.Site` object.
The concept of site families and sites allows `~kwant.builder.Builder`
to mix arbitrary lattices and site families
diff --git a/doc/source/tutorial/operators.rst b/doc/source/tutorial/operators.rst
index ed346f178fc3582a2f8909c8bc4cc76de03cad53..1fdd27871a4f3852df6b32ea92f1b9d633bb6897 100644
--- a/doc/source/tutorial/operators.rst
+++ b/doc/source/tutorial/operators.rst
@@ -438,7 +438,7 @@ the `~kwant.operator.Current`, instead of a constant matrix:
# evaluate the operator
m_current = J_m(psi, params=dict(r0=25, delta=10))
-The function must take a `~kwant.builder.Site` as its first parameter,
+The function must take a `~kwant.system.Site` as its first parameter,
and may optionally take other parameters (i.e. it must have the same
signature as a Hamiltonian onsite function), and must return the square
matrix that defines the operator we wish to calculate.
@@ -529,8 +529,8 @@ We see that the probability current is conserved across the scattering region,
but the z-projected spin current is not due to the fact that the Hamiltonian
does not commute with :math:`σ_z` everywhere in the scattering region.
-.. note:: ``where`` can also be provided as a sequence of `~kwant.builder.Site`
- or a sequence of hoppings (i.e. pairs of `~kwant.builder.Site`),
+.. note:: ``where`` can also be provided as a sequence of `~kwant.system.Site`
+ or a sequence of hoppings (i.e. pairs of `~kwant.system.Site`),
rather than a function.
diff --git a/doc/source/tutorial/plotting.rst b/doc/source/tutorial/plotting.rst
index bf6b51ae7e148fc6fdad8c37c5f24ed57421a079..48b6ed4796251b64aa9de10922ce911cd3aade86 100644
--- a/doc/source/tutorial/plotting.rst
+++ b/doc/source/tutorial/plotting.rst
@@ -105,7 +105,7 @@ the start end end site of hopping as arguments:
kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw);
Note that since we are using an unfinalized Builder, a ``site`` is really an
-instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
+instance of `~kwant.system.Site`. With these adjustments we arrive at a plot
that carries the same information, but is much easier to interpret.
Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
@@ -175,7 +175,7 @@ this is interpreted as the radius of the inner circle).
Finally, note that since we are dealing with a finalized system now, a site `i`
is represented by an integer. In order to obtain the original
-`~kwant.builder.Site`, ``syst.sites[i]`` can be used.
+`~kwant.system.Site`, ``syst.sites[i]`` can be used.
The way how data is presented of course influences what features of the data
are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
diff --git a/doc/source/tutorial/spin_potential_shape.rst b/doc/source/tutorial/spin_potential_shape.rst
index b87eaa7f6afced106110ca10823e48a6cb667069..23d90bf8bbf384af9d3644a34867c635c5eaf115 100644
--- a/doc/source/tutorial/spin_potential_shape.rst
+++ b/doc/source/tutorial/spin_potential_shape.rst
@@ -285,7 +285,7 @@ define the potential profile of a quantum well as:
else:
return 0
-This function takes two arguments: the first of type `~kwant.builder.Site`,
+This function takes two arguments: the first of type `~kwant.system.Site`,
from which you can get the real-space coordinates using ``site.pos``, and the
value of the potential as the second. Note that in `potential` we can access
variables `L` and `L_well` that are defined globally.
@@ -375,10 +375,10 @@ oscillatory transmission behavior through resonances in the quantum well.
.. specialnote:: Technical details
- Functions can also be used for hoppings. In this case, they take
- two `~kwant.builder.Site`'s as arguments and then an arbitrary number
+ two `~kwant.system.Site`'s as arguments and then an arbitrary number
of additional arguments.
- - Apart from the real-space position `pos`, `~kwant.builder.Site` has also an
+ - Apart from the real-space position `pos`, `~kwant.system.Site` has also an
attribute `tag` containing the lattice indices of the site.
.. _tutorial-abring:
@@ -449,7 +449,7 @@ that returns ``True`` whenever a point is inside the shape, and
return (r1 ** 2 < rsq < r2 ** 2)
Note that this function takes a real-space position as argument (not a
-`~kwant.builder.Site`).
+`~kwant.system.Site`).
We can now simply add all of the lattice points inside this shape at
once, using the function `~kwant.lattice.Square.shape`
diff --git a/kwant/_system.pyx b/kwant/_system.pyx
index 71860ed674dbabfade12a69c9eb2df9f31b6f65e..2f58dd4ff3bdb725acdf00238636006af337cde4 100644
--- a/kwant/_system.pyx
+++ b/kwant/_system.pyx
@@ -1,4 +1,4 @@
-# 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
@@ -12,12 +12,16 @@ import numpy as np
from scipy import sparse as sp
from itertools import chain
import types
+import bisect
from .graph.core cimport CGraph, gintArraySlice
from .graph.defs cimport gint
from .graph.defs import gint_dtype
from ._common import deprecate_args
+
+### Non-vectorized methods
+
msg = ('Hopping from site {0} to site {1} does not match the '
'dimensions of onsite Hamiltonians of these sites.')
@@ -352,3 +356,383 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
mat = func(ham, args, params, self.graph, diag, from_sites,
n_by_to_site, to_norb, to_off, from_norb, from_off)
return (mat, to_norb, from_norb) if return_norb else mat
+
+
+### Vectorized methods
+
+
+_shape_error_msg = (
+ "The following hoppings have matrix elements of incompatible shape "
+ "with other matrix elements: {}"
+)
+
+
+@cython.boundscheck(False)
+def _vectorized_make_sparse(subgraphs, hams, long [:] norbs, long [:] orb_offsets,
+ long [:] site_offsets, shape=None):
+ ndata = sum(h.shape[0] * h.shape[1] * h.shape[2] for h in hams)
+
+ cdef long [:] rows_view, cols_view
+ cdef complex [:] data_view
+
+ rows = np.empty((ndata,), dtype=long)
+ cols = np.empty((ndata,), dtype=long)
+ data = np.empty((ndata,), dtype=complex)
+ rows_view = rows
+ cols_view = cols
+ data_view = data
+
+ cdef long i, j, k, m, N, M, P, to_off, from_off,\
+ ta, fa, to_norbs, from_norbs
+ cdef long [:] ts_offs, fs_offs
+ cdef complex [:, :, :] h
+
+ m = 0
+ # This outer loop zip() is pure Python, but that's ok, as it
+ # has very few entries and the inner loops are fully vectorized
+ for ((ta, fa), (ts_offs, fs_offs)), h in zip(subgraphs, hams):
+ N = h.shape[0]
+ M = h.shape[1]
+ P = h.shape[2]
+
+ if norbs[ta] != M or norbs[fa] != P:
+ to_sites = site_offsets[ta] + np.array(ts_offs)
+ from_sites = site_offsets[fa] + np.array(fs_offs)
+ hops = np.array([to_sites, from_sites]).transpose()
+ raise ValueError(_shape_error_msg.format(hops))
+
+ for i in range(N):
+ to_off = orb_offsets[ta] + norbs[ta] * ts_offs[i]
+ from_off = orb_offsets[fa] + norbs[fa] * fs_offs[i]
+ for j in range(M):
+ for k in range(P):
+ rows_view[m] = to_off + j
+ cols_view[m] = from_off + k
+ data_view[m] = h[i, j, k]
+ m += 1
+
+ if shape is None:
+ shape = (orb_offsets[-1], orb_offsets[-1])
+
+ return sp.coo_matrix((data, (rows, cols)), shape=shape)
+
+
+@cython.boundscheck(False)
+def _vectorized_make_dense(subgraphs, hams, long [:] norbs, long [:] orb_offsets,
+ long [:] site_offsets, shape=None):
+ if shape is None:
+ shape = (orb_offsets[-1], orb_offsets[-1])
+ mat = np.zeros(shape, dtype=complex)
+ cdef complex [:, :] mat_view
+ mat_view = mat
+
+ cdef long i, j, k, N, M, P, to_off, from_off,\
+ ta, fa, to_norbs, from_norbs
+ cdef long [:] ts_offs, fs_offs
+ cdef complex [:, :, :] h
+
+ # This outer loop zip() is pure Python, but that's ok, as it
+ # has very few entries and the inner loops are fully vectorized
+ for ((ta, fa), (ts_offs, fs_offs)), h in zip(subgraphs, hams):
+ N = h.shape[0]
+ M = h.shape[1]
+ P = h.shape[2]
+
+ if norbs[ta] != M or norbs[fa] != P:
+ to_sites = site_offsets[ta] + np.array(ts_offs)
+ from_sites = site_offsets[fa] + np.array(fs_offs)
+ hops = np.array([to_sites, from_sites]).transpose()
+ raise ValueError(_shape_error_msg.format(hops))
+
+ for i in range(N):
+ to_off = orb_offsets[ta] + norbs[ta] * ts_offs[i]
+ from_off = orb_offsets[fa] + norbs[fa] * fs_offs[i]
+ for j in range(M):
+ for k in range(P):
+ mat_view[to_off + j, from_off + k] = h[i, j, k]
+ return mat
+
+
+def _count_norbs(syst, site_offsets, hams, args=(), params=None):
+ """Return the norbs and orbital offset of each site family in 'syst'
+
+ Parameters
+ ----------
+ syst : `kwant.system.System`
+ site_offsets : sequence of int
+ Contains the index of the first site for each site array
+ in 'syst.site_arrays'.
+ hams : sequence of arrays or 'None'
+ The Hamiltonian for each term in 'syst.terms'. 'None'
+ if the corresponding term has not been evaluated.
+ args, params : positional and keyword arguments to the system Hamiltonian
+ """
+ site_ranges = syst.site_ranges
+ if site_ranges is None:
+ # NOTE: Remove this branch once we can rely on
+ # site families storing the norb information.
+
+ site_arrays = syst.site_arrays
+ family_norbs = {s.family: None for s in site_arrays}
+
+ # Compute the norbs per site family using already evaluated
+ # Hamiltonian terms where possible
+ for ham, t in zip(hams, syst.terms):
+ (to_w, from_w), _ = syst.subgraphs[t.subgraph]
+ if ham is not None:
+ family_norbs[site_arrays[to_w].family] = ham.shape[1]
+ family_norbs[site_arrays[from_w].family] = ham.shape[2]
+
+ # Evaluate Hamiltonian terms where we do not already have them
+ for n, t in enumerate(syst.terms):
+ (to_w, from_w), _ = syst.subgraphs[t.subgraph]
+ to_fam = site_arrays[to_w].family
+ from_fam = site_arrays[from_w].family
+ if family_norbs[to_fam] is None or family_norbs[from_fam] is None:
+ ham = syst.hamiltonian_term(n, args=args, params=params)
+ family_norbs[to_fam] = ham.shape[1]
+ family_norbs[from_fam] = ham.shape[2]
+
+ # This case is technically allowed by the format (some sites present
+ # but no hamiltonian terms that touch them) but is very unlikely.
+ if any(norbs is None for norbs in family_norbs.values()):
+ raise ValueError("Cannot determine the number of orbitals for "
+ "some site families.")
+
+ orb_offset = 0
+ site_ranges = []
+ for start, site_array in zip(site_offsets, syst.site_arrays):
+ norbs = family_norbs[site_array.family]
+ site_ranges.append((start, norbs, orb_offset))
+ orb_offset += len(site_array) * norbs
+ site_ranges.append((site_offsets[-1], 0, orb_offset))
+ site_ranges = np.array(site_ranges)
+
+ _, norbs, orb_offsets = site_ranges.transpose()
+ # The final (extra) element in orb_offsets is the total number of orbitals
+ assert len(orb_offsets) == len(syst.site_arrays) + 1
+
+ return norbs, orb_offsets
+
+
+def _expand_norbs(compressed_norbs, site_offsets):
+ "Return norbs per site, given norbs per site array."
+ norbs = np.empty(site_offsets[-1], int)
+ for start, stop, norb in zip(site_offsets, site_offsets[1:],
+ compressed_norbs):
+ norbs[start:stop] = norb
+ return norbs
+
+
+def _reverse_subgraph(subgraph):
+ (to_sa, from_sa), (to_off, from_off) = subgraph
+ return ((from_sa, to_sa), (from_off, to_off))
+
+
+@deprecate_args
+@cython.binding(True)
+def vectorized_hamiltonian_submatrix(self, args=(), sparse=False,
+ return_norb=False, *, params=None):
+ """Return The system Hamiltonian.
+
+ Parameters
+ ----------
+ args : tuple, defaults to empty
+ Positional arguments to pass to ``hamiltonian_term``. Mutually
+ exclusive with 'params'.
+ sparse : bool
+ Whether to return a sparse or a dense matrix. Defaults to ``False``.
+ return_norb : bool
+ Whether to return arrays of numbers of orbitals. Defaults to ``False``.
+ params : dict, optional
+ Dictionary of parameter names and their values. Mutually exclusive
+ with 'args'.
+
+ Returns
+ -------
+ hamiltonian_part : numpy.ndarray or scipy.sparse.coo_matrix
+ The Hamiltonian of the system.
+ norb : array of integers
+ Numbers of orbitals on each site. Only returned when ``return_norb``
+ is true.
+
+ Notes
+ -----
+ Providing positional arguments via 'args' is deprecated,
+ instead, provide named parameters as a dictionary via 'params'.
+ """
+
+ site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
+
+ subgraphs = [
+ self.subgraphs[t.subgraph]
+ for t in self.terms
+ ]
+ # Add Hermitian conjugates
+ subgraphs += [
+ _reverse_subgraph(self.subgraphs[t.subgraph])
+ for t in self.terms
+ if t.hermitian
+ ]
+
+ hams = [
+ self.hamiltonian_term(n, args=args, params=params)
+ for n, _ in enumerate(self.terms)
+ ]
+ # Add Hermitian conjugates
+ hams += [
+ ham.conjugate().transpose(0, 2, 1)
+ for ham, t in zip(hams, self.terms)
+ if t.hermitian
+ ]
+
+ norbs, orb_offsets = _count_norbs(self, site_offsets, hams,
+ args=args, params=params)
+
+ func = _vectorized_make_sparse if sparse else _vectorized_make_dense
+ mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets)
+
+ if return_norb:
+ return (mat, _expand_norbs(norbs, site_offsets))
+ else:
+ return mat
+
+
+@deprecate_args
+@cython.binding(True)
+def vectorized_cell_hamiltonian(self, args=(), sparse=False, *, params=None):
+ """Hamiltonian of a single cell of the infinite system.
+
+ Providing positional arguments via 'args' is deprecated,
+ instead, provide named parameters as a dictionary via 'params'.
+ """
+
+ site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
+ # Site array where next cell starts
+ next_cell = bisect.bisect(site_offsets, self.cell_size) - 1
+
+ def inside_cell(term):
+ (to_which, from_which), _= self.subgraphs[term.subgraph]
+ return to_which < next_cell and from_which < next_cell
+
+ cell_terms = [
+ n for n, t in enumerate(self.terms)
+ if inside_cell(t)
+ ]
+
+ subgraphs = [
+ self.subgraphs[self.terms[n].subgraph]
+ for n in cell_terms
+ ]
+ # Add Hermitian conjugates
+ subgraphs += [
+ _reverse_subgraph(self.subgraphs[self.terms[n].subgraph])
+ for n in cell_terms
+ if self.terms[n].hermitian
+ ]
+
+ hams = [
+ self.hamiltonian_term(n, args=args, params=params)
+ for n in cell_terms
+ ]
+ # Add Hermitian conjugates
+ hams += [
+ ham.conjugate().transpose(0, 2, 1)
+ for ham, n in zip(hams, cell_terms)
+ if self.terms[n].hermitian
+ ]
+
+ # _count_norbs requires passing hamiltonians for all terms, or 'None'
+ # if it has not been evaluated
+ all_hams = [None] * len(self.terms)
+ for n, ham in zip(cell_terms, hams):
+ all_hams[n] = ham
+ norbs, orb_offsets = _count_norbs(self, site_offsets, all_hams,
+ args=args, params=params)
+
+ shape = (orb_offsets[next_cell], orb_offsets[next_cell])
+ func = _vectorized_make_sparse if sparse else _vectorized_make_dense
+ mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape)
+
+ return mat
+
+
+@deprecate_args
+@cython.binding(True)
+def vectorized_inter_cell_hopping(self, args=(), sparse=False, *, params=None):
+ """Hopping Hamiltonian between two cells of the infinite system.
+
+ Providing positional arguments via 'args' is deprecated,
+ instead, provide named parameters as a dictionary via 'params'.
+ """
+
+ # Take advantage of the fact that there are distinct entries in
+ # onsite_terms for sites inside the unit cell, and distinct entries
+ # in onsite_terms for hoppings between the unit cell sites and
+ # interface sites. This way we only need to check the first entry
+ # in onsite/hopping_terms
+
+ site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
+ # Site array where next cell starts
+ next_cell = bisect.bisect(site_offsets, self.cell_size) - 1
+
+ inter_cell_hopping_terms = [
+ n for n, t in enumerate(self.terms)
+ # *from* site is in interface
+ if (self.subgraphs[t.subgraph][0][1] >= next_cell
+ and self.subgraphs[t.subgraph][0][0] < next_cell)
+ ]
+ reversed_inter_cell_hopping_terms = [
+ n for n, t in enumerate(self.terms)
+ # *to* site is in interface
+ if (self.subgraphs[t.subgraph][0][0] >= next_cell
+ and self.subgraphs[t.subgraph][0][1] < next_cell)
+ ]
+
+ # Evaluate inter-cell hoppings only
+ inter_cell_hams = [
+ self.hamiltonian_term(n, args=args, params=params)
+ for n in inter_cell_hopping_terms
+ ]
+ reversed_inter_cell_hams = [
+ self.hamiltonian_term(n, args=args, params=params)
+ .conjugate().transpose(0, 2, 1)
+ for n in reversed_inter_cell_hopping_terms
+ ]
+
+ hams = inter_cell_hams + reversed_inter_cell_hams
+
+ subgraphs = [
+ self.subgraphs[self.terms[n].subgraph]
+ for n in inter_cell_hopping_terms
+ ]
+ subgraphs += [
+ _reverse_subgraph(self.subgraphs[self.terms[n].subgraph])
+ for n in reversed_inter_cell_hopping_terms
+ ]
+
+ # All the 'from' sites are in the previous domain, but to build a
+ # matrix we need to get the equivalent sites in the fundamental domain.
+ # We set the 'from' site array to the one from the fundamental domain.
+ subgraphs = [
+ ((to_sa, from_sa - next_cell), (to_off, from_off))
+ for (to_sa, from_sa), (to_off, from_off) in subgraphs
+ ]
+
+ # _count_norbs requires passing hamiltonians for all terms, or 'None'
+ # if it has not been evaluated
+ all_hams = [None] * len(self.terms)
+ for n, ham in zip(inter_cell_hopping_terms, inter_cell_hams):
+ all_hams[n] = ham
+ for n, ham in zip(reversed_inter_cell_hopping_terms,
+ reversed_inter_cell_hams):
+ # Transpose to get back correct shape wrt. original term
+ all_hams[n] = ham.transpose(0, 2, 1)
+ norbs, orb_offsets = _count_norbs(self, site_offsets, all_hams,
+ args=args, params=params)
+
+ shape = (orb_offsets[next_cell],
+ orb_offsets[len(self.site_arrays) - next_cell])
+ func = _vectorized_make_sparse if sparse else _vectorized_make_dense
+ mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape)
+ return mat
diff --git a/kwant/builder.py b/kwant/builder.py
index 390ca0ba3c54d30d080d67c5002b43a49b74f4bf..31019c34533e879d383192f377e857b4e74f06d1 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -1,4 +1,4 @@
-# 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
@@ -14,11 +14,14 @@ import copy
from functools import total_ordering, wraps, update_wrapper
from itertools import islice, chain
import textwrap
+import bisect
+import numbers
import inspect
import tinyarray as ta
import numpy as np
from scipy import sparse
from . import system, graph, KwantDeprecationWarning, UserCodeError
+from .system import Site, SiteArray, SiteFamily
from .linalg import lll
from .operator import Density
from .physics import DiscreteSymmetry, magnetic_gauge
@@ -26,182 +29,13 @@ from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
interleave, deprecate_args, memoize)
-__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
- 'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead',
- 'add_peierls_phase']
+__all__ = ['Builder', 'SimpleSiteFamily', 'Symmetry', 'HoppingKind', 'Lead',
+ 'BuilderLead', 'SelfEnergyLead', 'ModesLead', 'add_peierls_phase']
-################ Sites and site families
-
-class Site(tuple):
- """A site, member of a `SiteFamily`.
-
- Sites are the vertices of the graph which describes the tight binding
- system in a `Builder`.
-
- A site is uniquely identified by its family and its tag.
-
- Parameters
- ----------
- family : an instance of `SiteFamily`
- The 'type' of the site.
- tag : a hashable python object
- The unique identifier of the site within the site family, typically a
- vector of integers.
-
- Raises
- ------
- ValueError
- If `tag` is not a proper tag for `family`.
-
- Notes
- -----
- For convenience, ``family(*tag)`` can be used instead of ``Site(family,
- tag)`` to create a site.
-
- The parameters of the constructor (see above) are stored as instance
- variables under the same names. Given a site ``site``, common things to
- query are thus ``site.family``, ``site.tag``, and ``site.pos``.
- """
- __slots__ = ()
-
- family = property(operator.itemgetter(0),
- doc="The site family to which the site belongs.")
- tag = property(operator.itemgetter(1), doc="The tag of the site.")
-
-
- def __new__(cls, family, tag, _i_know_what_i_do=False):
- if _i_know_what_i_do:
- return tuple.__new__(cls, (family, tag))
- try:
- tag = family.normalize_tag(tag)
- except (TypeError, ValueError) as e:
- msg = 'Tag {0} is not allowed for site family {1}: {2}'
- raise type(e)(msg.format(repr(tag), repr(family), e.args[0]))
- return tuple.__new__(cls, (family, tag))
-
- def __repr__(self):
- return 'Site({0}, {1})'.format(repr(self.family), repr(self.tag))
-
- def __str__(self):
- sf = self.family
- return '<Site {0} of {1}>'.format(self.tag, sf.name if sf.name else sf)
-
- def __getnewargs__(self):
- return (self.family, self.tag, True)
-
- @property
- def pos(self):
- """Real space position of the site.
-
- This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
- """
- return self.family.pos(self.tag)
-
+################ Site families
@total_ordering
-class SiteFamily(metaclass=abc.ABCMeta):
- """Abstract base class for site families.
-
- Site families are the 'type' of `Site` objects. Within a family, individual
- sites are uniquely identified by tags. Valid tags must be hashable Python
- objects, further details are up to the family.
-
- Site families must be immutable and fully defined by their initial
- arguments. They must inherit from this abstract base class and call its
- __init__ function providing it with two arguments: a canonical
- representation and a name. The canonical representation will be returned as
- the objects representation and must uniquely identify the site family
- instance. The name is a string used to distinguish otherwise identical site
- families. It may be empty. ``norbs`` defines the number of orbitals
- on sites associated with this site family; it may be `None`, in which case
- the number of orbitals is not specified.
-
-
- All site families must define the method `normalize_tag` which brings a tag
- to the standard format for this site family.
-
- Site families that are intended for use with plotting should also provide a
- method `pos(tag)`, which returns a vector with real-space coordinates of the
- site belonging to this family with a given tag.
-
- If the ``norbs`` of a site family are provided, and sites of this family
- are used to populate a `~kwant.builder.Builder`, then the associated
- Hamiltonian values must have the correct shape. That is, if a site family
- has ``norbs = 2``, then any on-site terms for sites belonging to this
- family should be 2x2 matrices. Similarly, any hoppings to/from sites
- belonging to this family must have a matrix structure where there are two
- rows/columns. This condition applies equally to Hamiltonian values that
- are given by functions. If this condition is not satisfied, an error will
- be raised.
- """
-
- def __init__(self, canonical_repr, name, norbs):
- self.canonical_repr = canonical_repr
- self.hash = hash(canonical_repr)
- self.name = name
- if norbs is None:
- warnings.warn("Not specfying norbs is deprecated. Always specify "
- "norbs when creating site families.",
- KwantDeprecationWarning, stacklevel=3)
- if norbs is not None:
- if int(norbs) != norbs or norbs <= 0:
- raise ValueError('The norbs parameter must be an integer > 0.')
- norbs = int(norbs)
- self.norbs = norbs
-
- def __repr__(self):
- return self.canonical_repr
-
- def __str__(self):
- if self.name:
- msg = '<{0} site family {1}{2}>'
- else:
- msg = '<unnamed {0} site family{2}>'
- orbs = ' with {0} orbitals'.format(self.norbs) if self.norbs else ''
- return msg.format(self.__class__.__name__, self.name, orbs)
-
- def __hash__(self):
- return self.hash
-
- def __eq__(self, other):
- try:
- return self.canonical_repr == other.canonical_repr
- except AttributeError:
- return False
-
- def __ne__(self, other):
- try:
- return self.canonical_repr != other.canonical_repr
- except AttributeError:
- return True
-
- def __lt__(self, other):
- # If this raises an AttributeError, we were trying
- # to compare it to something non-comparable anyway.
- return self.canonical_repr < other.canonical_repr
-
- @abc.abstractmethod
- def normalize_tag(self, tag):
- """Return a normalized version of the tag.
-
- Raises TypeError or ValueError if the tag is not acceptable.
- """
- pass
-
- def __call__(self, *tag):
- """
- A convenience function.
-
- This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).
- """
- # Catch a likely and difficult to find mistake.
- if tag and isinstance(tag[0], tuple):
- raise ValueError('Use site_family(1, 2) instead of '
- 'site_family((1, 2))!')
- return Site(self, tag)
-
-
class SimpleSiteFamily(SiteFamily):
"""A site family used as an example and for testing.
@@ -307,19 +141,44 @@ class Symmetry(metaclass=abc.ABCMeta):
def which(self, site):
"""Calculate the domain of the site.
- Return the group element whose action on a certain site from the
- fundamental domain will result in the given ``site``.
+ Parameters
+ ----------
+ site : `~kwant.system.Site` or `~kwant.system.SiteArray`
+
+ Returns
+ -------
+ group_element : tuple or sequence of tuples
+ A single tuple if ``site`` is a Site, or a sequence of tuples if
+ ``site`` is a SiteArray. The group element(s) whose action
+ on a certain site(s) from the fundamental domain will result
+ in the given ``site``.
"""
pass
@abc.abstractmethod
def act(self, element, a, b=None):
- """Act with a symmetry group element on a site or hopping."""
+ """Act with symmetry group element(s) on site(s) or hopping(s).
+
+ Parameters
+ ----------
+ element : tuple or sequence of tuples
+ Group element(s) with which to act on the provided site(s)
+ or hopping(s)
+ a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
+ If Site then ``element`` is a single tuple, if SiteArray then
+ ``element`` is a sequence of tuples. If only ``a`` is provided then
+ ``element`` acts on the site(s) of ``a``. If ``b`` is also provided
+ then ``element`` acts on the hopping(s) ``(a, b)``.
+ """
pass
def to_fd(self, a, b=None):
"""Map a site or hopping to the fundamental domain.
+ Parameters
+ ----------
+ a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
+
If ``b`` is None, return a site equivalent to ``a`` within the
fundamental domain. Otherwise, return a hopping equivalent to ``(a,
b)`` but where the first element belongs to the fundamental domain.
@@ -329,11 +188,30 @@ class Symmetry(metaclass=abc.ABCMeta):
return self.act(-self.which(a), a, b)
def in_fd(self, site):
- """Tell whether ``site`` lies within the fundamental domain."""
- for d in self.which(site):
- if d != 0:
- return False
- return True
+ """Tell whether ``site`` lies within the fundamental domain.
+
+ Parameters
+ ----------
+ site : `~kwant.system.Site` or `~kwant.system.SiteArray`
+
+ Returns
+ -------
+ in_fd : bool or sequence of bool
+ single bool if ``site`` is a Site, or a sequence of
+ bool if ``site`` is a SiteArray. In the latter case
+ we return whether each site in the SiteArray is in
+ the fundamental domain.
+ """
+ if isinstance(site, Site):
+ for d in self.which(site):
+ if d != 0:
+ return False
+ return True
+ elif isinstance(site, SiteArray):
+ which = self.which(site)
+ return np.logical_and.reduce(which != 0, axis=1)
+ else:
+ raise TypeError("'site' must be a Site or SiteArray")
@abc.abstractmethod
def subgroup(self, *generators):
@@ -412,8 +290,8 @@ class HoppingKind(tuple):
----------
delta : Sequence of integers
The sequence is interpreted as a vector with integer elements.
- family_a : `~kwant.builder.SiteFamily`
- family_b : `~kwant.builder.SiteFamily` or ``None`` (default)
+ family_a : `~kwant.system.SiteFamily`
+ family_b : `~kwant.system.SiteFamily` or ``None`` (default)
The default value means: use the same family as `family_a`.
Notes
@@ -570,7 +448,7 @@ class BuilderLead(Lead):
The tight-binding system of a lead. It has to possess appropriate
symmetry, and it may not contain hoppings between further than
neighboring images of the fundamental domain.
- interface : sequence of `Site` instances
+ interface : sequence of `~kwant.system.Site` instances
Sequence of sites in the scattering region to which the lead is
attached.
@@ -578,9 +456,9 @@ class BuilderLead(Lead):
----------
builder : `Builder`
The tight-binding system of a lead.
- interface : list of `Site` instances
+ interface : list of `~kwant.system.Site` instances
A sorted list of interface sites.
- padding : list of `Site` instances
+ padding : list of `~kwant.system.Site` instances
A sorted list of sites that originate from the lead, have the same
onsite Hamiltonian, and are connected by the same hoppings as the lead
sites.
@@ -607,7 +485,10 @@ class BuilderLead(Lead):
The order of interface sites is kept during finalization.
"""
- return InfiniteSystem(self.builder, self.interface)
+ if self.builder.vectorize:
+ return InfiniteVectorizedSystem(self.builder, self.interface)
+ else:
+ return InfiniteSystem(self.builder, self.interface)
def _ensure_signature(func):
@@ -638,7 +519,7 @@ class SelfEnergyLead(Lead):
selfenergy_func : function
Has the same signature as `selfenergy` (without the ``self``
parameter) and returns the self energy matrix for the interface sites.
- interface : sequence of `Site` instances
+ interface : sequence of `~kwant.system.Site` instances
parameters : sequence of strings
The parameters on which the lead depends.
"""
@@ -669,7 +550,7 @@ class ModesLead(Lead):
and returns the modes of the lead as a tuple of
`~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes`.
interface :
- sequence of `Site` instances
+ sequence of `~kwant.system.Site` instances
parameters : sequence of strings
The parameters on which the lead depends.
"""
@@ -797,10 +678,11 @@ class Builder:
This is one of the central types in Kwant. It is used to construct tight
binding systems in a flexible way.
- The nodes of the graph are `Site` instances. The edges, i.e. the hoppings,
- are pairs (2-tuples) of sites. Each node and each edge has a value
- associated with it. The values associated with nodes are interpreted as
- on-site Hamiltonians, the ones associated with edges as hopping integrals.
+ The nodes of the graph are `~kwant.system.Site` instances. The edges,
+ i.e. the hoppings, are pairs (2-tuples) of sites. Each node and each
+ edge has a value associated with it. The values associated with nodes
+ are interpreted as on-site Hamiltonians, the ones associated with edges
+ as hopping integrals.
To make the graph accessible in a way that is natural within the Python
language it is exposed as a *mapping* (much like a built-in Python
@@ -827,6 +709,14 @@ class Builder:
chiral : 2D array, dictionary, function or `None`
The unitary part of the onsite chiral symmetry operator.
Same format as that of `conservation_law`.
+ vectorize : bool, default: False
+ If True then the finalized Builder will evaluate its Hamiltonian in
+ a vectorized manner. This requires that all value functions take
+ `~kwant.system.SiteArray` as first/second parameter rather than
+ `~kwant.system.Site`, and return an array of values. The returned
+ array must have the same length as the provided SiteArray(s), and
+ may contain either scalars (i.e. a vector of values) or matrices
+ (i.e. a 3d array of values).
Notes
-----
@@ -905,7 +795,7 @@ class Builder:
"""
def __init__(self, symmetry=None, *, conservation_law=None, time_reversal=None,
- particle_hole=None, chiral=None):
+ particle_hole=None, chiral=None, vectorize=False):
if symmetry is None:
symmetry = NoSymmetry()
else:
@@ -915,6 +805,7 @@ class Builder:
self.time_reversal = time_reversal
self.particle_hole = particle_hole
self.chiral = chiral
+ self.vectorize = vectorize
self.leads = []
self.H = {}
@@ -998,6 +889,7 @@ class Builder:
result.time_reversal = self.time_reversal
result.particle_hole = self.particle_hole
result.chiral = self.chiral
+ result.vectorize = self.vectorize
result.leads = self.leads
result.H = self.H
return result
@@ -1452,7 +1344,7 @@ class Builder:
A boolean function of site returning whether the site should be
added to the target builder or not. The shape must be compatible
with the symmetry of the target builder.
- start : `Site` instance or iterable thereof or iterable of numbers
+ start : `~kwant.system.Site` or iterable thereof or iterable of numbers
The site(s) at which the the flood-fill starts. If start is an
iterable of numbers, the starting site will be
``template.closest(start)``.
@@ -1462,7 +1354,7 @@ class Builder:
Returns
-------
- added_sites : list of `Site` objects that were added to the system.
+ added_sites : list of `~kwant.system.Site` that were added to the system.
Raises
------
@@ -1614,7 +1506,7 @@ class Builder:
----------
lead_builder : `Builder` with 1D translational symmetry
Builder of the lead which has to be attached.
- origin : `Site`
+ origin : `~kwant.system.Site`
The site which should belong to a domain where the lead should
begin. It is used to attach a lead inside the system, e.g. to an
inner radius of a ring.
@@ -1624,7 +1516,7 @@ class Builder:
Returns
-------
- added_sites : list of `Site` objects that were added to the system.
+ added_sites : list of `~kwant.system.Site`
Raises
------
@@ -1677,6 +1569,13 @@ class Builder:
if add_cells < 0 or int(add_cells) != add_cells:
raise ValueError('add_cells must be an integer >= 0.')
+ if self.vectorize != lead_builder.vectorize:
+ raise ValueError(
+ "Sites of the lead were added to the scattering "
+ "region, but only one of these systems is vectorized. "
+ "Vectorize both systems to allow attaching leads."
+ )
+
sym = lead_builder.symmetry
H = lead_builder.H
@@ -1697,6 +1596,7 @@ class Builder:
time_reversal=lead_builder.time_reversal,
particle_hole=lead_builder.particle_hole,
chiral=lead_builder.chiral,
+ vectorize=lead_builder.vectorize,
)
with reraise_warnings():
new_lead.fill(lead_builder, lambda site: True,
@@ -1793,9 +1693,15 @@ class Builder:
`Symmetry` can be finalized.
"""
if self.symmetry.num_directions == 0:
- return FiniteSystem(self)
+ if self.vectorize:
+ return FiniteVectorizedSystem(self)
+ else:
+ return FiniteSystem(self)
elif self.symmetry.num_directions == 1:
- return InfiniteSystem(self)
+ if self.vectorize:
+ return InfiniteVectorizedSystem(self)
+ else:
+ return InfiniteSystem(self)
else:
raise ValueError('Currently, only builders without or with a 1D '
'translational symmetry can be finalized.')
@@ -1969,6 +1875,10 @@ def add_peierls_phase(syst, peierls_parameter='phi', fix_gauge=True):
return f
+ if syst.vectorize:
+ raise TypeError("'add_peierls_phase' does not work with "
+ "vectorized Builders")
+
ret = _add_peierls_phase(syst, peierls_parameter).finalized()
if fix_gauge:
@@ -2115,6 +2025,135 @@ class _FinalizedBuilderMixin:
return DiscreteSymmetry(projectors, *(evaluate(symm) for symm in
self._symmetries))
+ def pos(self, i):
+ return self.sites[i].pos
+
+
+class _VectorizedFinalizedBuilderMixin(_FinalizedBuilderMixin):
+ """Common functionality for all vectorized finalized builders
+
+ Attributes
+ ----------
+ _term_values : sequence
+ Each item is either an array of values (if 'param_names is None')
+ or a vectorized value function (in which case 'param_names' is a
+ sequence of strings: the extra parameters to the value function).
+ _onsite_term_by_site_id : sequence of int
+ Maps site (integer) to the term that evaluates the associated
+ onsite matrix element.
+ _hopping_term_by_edge_id : sequence of int
+ Maps edge id in the graph to the term that evaluates the associated
+ hopping matrix element. If negative, then the associated hopping is
+ the Hermitian conjugate of another hopping; if the number stored is
+ 't' (< 0) then the associated hopping is stored in term '-t - 1'.
+ """
+ def hamiltonian(self, i, j, *args, params=None):
+ warnings.warn(
+ "Querying individual matrix elements with 'hamiltonian' is "
+ "deprecated, and now takes O(log N) time where N is the number "
+ "of matrix elements per hamiltonian term. Query many matrix "
+ "elements at once using 'hamiltonian_term'.",
+ KwantDeprecationWarning
+ )
+ site_offsets = np.cumsum([0] + [len(s) for s in self.site_arrays])
+ if i == j:
+ which_term = self._onsite_term_by_site_id[i]
+ (w, _), (off, _) = self.subgraphs[self.terms[which_term].subgraph]
+ i_off = i - site_offsets[w]
+ selector = np.searchsorted(off, i_off)
+ onsite = self.hamiltonian_term(
+ which_term, [selector], args, params=params)
+ return onsite[0]
+ else:
+ edge_id = self.graph.first_edge_id(i, j)
+ which_term = self._hopping_term_by_edge_id[edge_id]
+ herm_conj = which_term < 0
+ if herm_conj:
+ which_term = -which_term - 1
+ term = self.terms[which_term]
+ # To search for hopping (i, j) in hopping_terms, we need
+ # to treat hopping_terms as a record array of integer pairs
+ dtype = np.dtype([('f0', int), ('f1', int)])
+ # For hermitian conjugate terms search through the
+ # *other* term's hoppings because those are sorted.
+
+ (to_w, from_w), hoppings = self.subgraphs[term.subgraph]
+ hoppings = hoppings.transpose() # to get array-of-pairs
+ if herm_conj:
+ hop = (j - site_offsets[to_w], i - site_offsets[from_w])
+ else:
+ hop = (i - site_offsets[to_w], j - site_offsets[from_w])
+ hop = np.array(hop, dtype=dtype)
+ hoppings = hoppings.view(dtype).reshape(-1)
+ selector = np.recarray.searchsorted(hoppings, hop)
+ h = self.hamiltonian_term(
+ which_term, [selector], args, params=params)
+ h = h[0]
+ if herm_conj:
+ h = h.conjugate().transpose()
+ return h
+
+ def hamiltonian_term(self, term_number, selector=slice(None),
+ args=(), params=None):
+ if args and params:
+ raise TypeError("'args' and 'params' are mutually exclusive.")
+
+ term = self.terms[term_number]
+ val = self._term_values[term_number]
+
+ if not callable(val):
+ return val[selector]
+
+ # Construct site arrays to pass to the vectorized value function.
+ subgraph = self.subgraphs[self.terms[term_number].subgraph]
+ (to_which, from_which), (to_off, from_off) = subgraph
+ is_onsite = to_off is from_off
+ to_off = to_off[selector]
+ from_off = from_off[selector]
+ assert len(to_off) == len(from_off)
+
+ to_family = self.site_arrays[to_which].family
+ to_tags = self.site_arrays[to_which].tags
+ to_site_array = SiteArray(to_family, to_tags[to_off])
+ if not is_onsite:
+ from_family = self.site_arrays[from_which].family
+ from_tags = self.site_arrays[from_which].tags
+ from_site_array = SiteArray(from_family, from_tags[from_off])
+
+ # Construct args from params
+ if params:
+ # There was a problem extracting parameter names from the value
+ # function (probably an illegal signature) and we are using
+ # keyword parameters.
+ if self._term_errors[term_number] is not None:
+ raise self._term_errors[term_number]
+ try:
+ args = [params[p] for p in term.parameters]
+ except KeyError:
+ missing = [p for p in term.parameters if p not in params]
+ msg = ('System is missing required arguments: ',
+ ', '.join(map('"{}"'.format, missing)))
+ raise TypeError(''.join(msg))
+
+ if is_onsite:
+ try:
+ ham = val(to_site_array, *args)
+ except Exception as exc:
+ _raise_user_error(exc, val)
+ else:
+ try:
+ ham = val(
+ *self.symmetry.to_fd(
+ to_site_array,
+ from_site_array),
+ *args)
+ except Exception as exc:
+ _raise_user_error(exc, val)
+
+ ham = _normalize_term_value(ham, len(to_site_array))
+
+ return ham
+
# The same (value, parameters) pair will be used for many sites/hoppings,
# so we cache it to avoid wasting extra memory.
@@ -2147,6 +2186,64 @@ def _value_params_pair_cache(nstrip):
return get
+def _make_graph(H, id_by_site):
+ g = graph.Graph()
+ g.num_nodes = len(id_by_site) # Some sites could not appear in any edge.
+ for tail, hvhv in H.items():
+ for head in islice(hvhv, 2, None, 2):
+ if tail == head:
+ continue
+ g.add_edge(id_by_site[tail], id_by_site[head])
+ return g.compressed()
+
+
+def _finalize_leads(leads, id_by_site):
+ #### Connect leads.
+ finalized_leads = []
+ lead_interfaces = []
+ lead_paddings = []
+ for lead_nr, lead in enumerate(leads):
+ try:
+ with warnings.catch_warnings(record=True) as ws:
+ warnings.simplefilter("always")
+ # The following line is the whole "payload" of the entire
+ # try-block.
+ finalized_leads.append(lead.finalized())
+ except ValueError as e:
+ # Re-raise the exception with an additional message.
+ msg = 'Problem finalizing lead {0}:'.format(lead_nr)
+ e.args = (' '.join((msg,) + e.args),)
+ raise
+ else:
+ for w in ws:
+ # Re-raise any warnings with an additional message and the
+ # proper stacklevel.
+ w = w.message
+ msg = 'When finalizing lead {0}:'.format(lead_nr)
+ warnings.warn(w.__class__(' '.join((msg,) + w.args)),
+ stacklevel=3)
+ try:
+ interface = [id_by_site[isite] for isite in lead.interface]
+ except KeyError as e:
+ msg = ("Lead {0} is attached to a site that does not "
+ "belong to the scattering region:\n {1}")
+ raise ValueError(msg.format(lead_nr, e.args[0]))
+
+ lead_interfaces.append(np.array(interface))
+
+ padding = getattr(lead, 'padding', [])
+ # Some padding sites might have been removed after the lead was
+ # attached. Unlike in the case of the interface, this is not a
+ # problem.
+ finalized_padding = [
+ id_by_site[isite] for isite in padding if isite in id_by_site
+ ]
+
+ lead_paddings.append(np.array(finalized_padding))
+
+ return finalized_leads, lead_interfaces, lead_paddings
+
+
class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
"""Finalized `Builder` with leads.
@@ -2155,11 +2252,11 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
Attributes
----------
sites : sequence
- ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
+ ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
to the integer-labeled site ``i`` of the low-level system. The sites
are ordered first by their family and then by their tag.
- id_by_site : dict
- The inverse of ``sites``; maps high-level `~kwant.builder.Site`
+ id_by_site : mapping
+ The inverse of ``sites``; maps high-level `~kwant.system.Site`
instances to their integer label.
Satisfies ``id_by_site[sites[i]] == i``.
"""
@@ -2173,57 +2270,10 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
for site_id, site in enumerate(sites):
id_by_site[site] = site_id
- #### Make graph.
- g = graph.Graph()
- g.num_nodes = len(sites) # Some sites could not appear in any edge.
- for tail, hvhv in builder.H.items():
- for head in islice(hvhv, 2, None, 2):
- if tail == head:
- continue
- g.add_edge(id_by_site[tail], id_by_site[head])
- g = g.compressed()
-
- #### Connect leads.
- finalized_leads = []
- lead_interfaces = []
- lead_paddings = []
- for lead_nr, lead in enumerate(builder.leads):
- try:
- with warnings.catch_warnings(record=True) as ws:
- warnings.simplefilter("always")
- # The following line is the whole "payload" of the entire
- # try-block.
- finalized_leads.append(lead.finalized())
- for w in ws:
- # Re-raise any warnings with an additional message and the
- # proper stacklevel.
- w = w.message
- msg = 'When finalizing lead {0}:'.format(lead_nr)
- warnings.warn(w.__class__(' '.join((msg,) + w.args)),
- stacklevel=3)
- except ValueError as e:
- # Re-raise the exception with an additional message.
- msg = 'Problem finalizing lead {0}:'.format(lead_nr)
- e.args = (' '.join((msg,) + e.args),)
- raise
- try:
- interface = [id_by_site[isite] for isite in lead.interface]
- except KeyError as e:
- msg = ("Lead {0} is attached to a site that does not "
- "belong to the scattering region:\n {1}")
- raise ValueError(msg.format(lead_nr, e.args[0]))
-
- lead_interfaces.append(np.array(interface))
-
- padding = getattr(lead, 'padding', [])
- # Some padding sites might have been removed after the lead was
- # attached. Unlike in the case of the interface, this is not a
- # problem.
- finalized_padding = [
- id_by_site[isite] for isite in padding if isite in id_by_site
- ]
+ graph = _make_graph(builder.H, id_by_site)
- lead_paddings.append(np.array(finalized_padding))
+ finalized_leads, lead_interfaces, lead_paddings =\
+ _finalize_leads(builder.leads, id_by_site)
# Because many onsites/hoppings share the same (value, parameter)
# pairs, we keep them in a cache so that we only store a given pair
@@ -2232,7 +2282,7 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
onsites = [cache(builder.H[site][1]) for site in sites]
cache = _value_params_pair_cache(2)
hoppings = [cache(builder._get_edge(sites[tail], sites[head]))
- for tail, head in g]
+ for tail, head in graph]
# Compute the union of the parameters of onsites and hoppings. Here,
# 'onsites' and 'hoppings' are pairs whose second element is one of
@@ -2252,7 +2302,7 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
else:
parameters = frozenset(parameters)
- self.graph = g
+ self.graph = graph
self.sites = sites
self.site_ranges = _site_ranges(sites)
self.id_by_site = id_by_site
@@ -2265,8 +2315,502 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
self.lead_paddings = lead_paddings
self._init_discrete_symmetries(builder)
- def pos(self, i):
- return self.sites[i].pos
+
+def _make_site_arrays(sites):
+ tags_by_family = {}
+ for family, tag in sites: # Sites are tuples
+ tags_by_family.setdefault(family, []).append(tag)
+ site_arrays = []
+ for family, tags in sorted(tags_by_family.items()):
+ site_arrays.append(SiteArray(family, sorted(tags)))
+ return site_arrays
+
+
+# Wrapper for accessing sites by their sequential number from a
+# list of SiteArrays.
+class _Sites:
+
+ def __init__(self, site_arrays):
+ self._site_arrays = site_arrays
+ lengths = [0] + [len(arr.tags) for arr in site_arrays]
+ self._start_idxs = np.cumsum(lengths)[:-1]
+
+ def __len__(self):
+ return sum(len(arr.tags) for arr in self._site_arrays)
+
+ def __getitem__(self, idx):
+ if not isinstance(idx, numbers.Integral):
+ raise TypeError("Only individual sites may be indexed")
+ if idx < 0 or idx >= len(self):
+ raise IndexError(idx)
+
+ which = bisect.bisect(self._start_idxs, idx) - 1
+ arr = self._site_arrays[which]
+ start = self._start_idxs[which]
+ fam = arr.family
+ tag = arr.tags[idx - start]
+ return Site(fam, tag)
+
+ def __iter__(self):
+ for arr in self._site_arrays:
+ for tag in arr.tags:
+ yield Site(arr.family, tag)
+
+
+# Wrapper for accessing the sequential number of a site, given
+# Site object, from a list of SiteArrays.
+class _IdBySite:
+
+ def __init__(self, site_arrays):
+ self._site_arrays = site_arrays
+ lengths = [0] + [len(arr.tags) for arr in site_arrays]
+ self._start_idxs = np.cumsum(lengths)[:-1]
+
+ def __len__(self):
+ return sum(len(arr.tags) for arr in self._site_arrays)
+
+ def __getitem__(self, site):
+ # treat tags as 1D sequence by defining custom dtype
+ tag_dtype = np.asarray(site.tag).dtype
+ dtype = np.dtype([('f{}'.format(i), tag_dtype)
+ for i in range(len(site.tag))])
+ site_tag = np.asarray(site.tag).view(dtype)[0]
+ # This slightly convoluted logic is necessary because there
+ # may be >1 site_array with the same family for InfiniteSystems.
+ for start, arr in zip(self._start_idxs, self._site_arrays):
+ if arr.family != site.family:
+ continue
+ tags = arr.tags.view(dtype).reshape(-1)
+ idx_in_fam = np.recarray.searchsorted(tags, site_tag)
+ if idx_in_fam >= len(tags) or tags[idx_in_fam] != site_tag:
+ continue
+ return start + idx_in_fam
+
+ raise KeyError(site)
+
+
+def _normalize_term_value(value, expected_n_values):
+ try:
+ value = np.asarray(value, dtype=complex)
+ except TypeError:
+ raise ValueError(
+ "Matrix elements declared with incompatible shapes."
+ ) from None
+ # Upgrade to vector of matrices
+ if len(value.shape) == 1:
+ value = value[:, np.newaxis, np.newaxis]
+ if len(value.shape) != 3:
+ msg = (
+ "Vectorized value functions must return an array of"
+ "scalars or an array of matrices."
+ )
+ raise ValueError(msg)
+ if value.shape[0] != expected_n_values:
+ raise ValueError("Value functions must return a single value per "
+ "onsite/hopping.")
+ return value
+
+
+def _sort_term(term, value):
+ term = np.asarray(term)
+
+ if not callable(value):
+ value = _normalize_term_value(value, len(term))
+ # Ensure that values still correspond to the correct
+ # sites in 'term' once the latter has been sorted.
+ value = value[term.argsort()]
+
+ term.sort()
+
+ return term, value
+
+
+def _sort_hopping_term(term, value):
+ term = np.asarray(term)
+ # We want to sort the hopping terms in the same way
+ # as tuples (i, j), so we need to use a record array.
+ dtype = np.dtype([('f0', term.dtype), ('f1', term.dtype)])
+ term_prime = term.view(dtype=dtype).reshape(-1)
+ # _normalize_term will sort 'term' in-place via 'term_prime'
+ _, value = _sort_term(term_prime, value)
+
+ return term, value
+
+
+def _make_onsite_terms(builder, sites, site_offsets, term_offset):
+ # Construct onsite terms.
+ #
+ # onsite_subgraphs
+ # Will contain tuples of the form described in
+ # kwant.system.VectorizedSystem.subgraphs, but during construction
+ # contains lists of the sites associated with each onsite term.
+ #
+ # onsite_term_values
+ # Contains a callable or array of constant values for each term.
+ #
+ # onsite_terms
+ # Constructed after the main loop. Contains a kwant.system.Term
+ # tuple for each onsite term.
+ #
+ # _onsite_term_by_site_id
+ # Maps the site ID to the number of the term that the site is
+ # a part of.
+ # lists the number of the
+ # Hamiltonian term associated with each site/hopping. For
+ # Hermitian conjugate hoppings "-term - 1" is stored instead.
+
+ onsite_subgraphs = []
+ onsite_term_values = []
+ onsite_term_parameters = []
+ # We just use the cache to handle non/callables and exceptions
+ cache = _value_params_pair_cache(1)
+ # maps onsite key -> onsite term number
+ onsite_to_term_nr = collections.OrderedDict()
+ _onsite_term_by_site_id = []
+ for site_id, site in enumerate(sites):
+ val = builder.H[builder.symmetry.to_fd(site)][1]
+ const_val = not callable(val)
+ which_site_array = bisect.bisect(site_offsets, site_id) - 1
+ # "key" uniquely identifies an onsite term.
+ if const_val:
+ key = (None, which_site_array)
+ else:
+ key = (id(val), which_site_array)
+ if key not in onsite_to_term_nr:
+ # Start a new term
+ onsite_to_term_nr[key] = len(onsite_subgraphs)
+ onsite_subgraphs.append([])
+ if const_val:
+ onsite_term_values.append([])
+ onsite_term_parameters.append(None)
+ else:
+ val, params = cache(val)
+ onsite_term_values.append(val)
+ onsite_term_parameters.append(params)
+ onsite_subgraphs[onsite_to_term_nr[key]].append(site_id)
+ _onsite_term_by_site_id.append(onsite_to_term_nr[key])
+ if const_val:
+ vals = onsite_term_values[onsite_to_term_nr[key]]
+ vals.append(val)
+ # Sort the sites in each term, and normalize any constant
+ # values to arrays of the correct dtype and shape.
+ onsite_subgraphs, onsite_term_values = zip(*(
+ _sort_term(term, value)
+ for term, value in
+ zip(onsite_subgraphs, onsite_term_values)))
+ # Store site array index and site offsets rather than sites themselves
+ tmp = []
+ for (_, which), s in zip(onsite_to_term_nr, onsite_subgraphs):
+ s = s - site_offsets[which]
+ tmp.append(((which, which), (s, s)))
+ onsite_subgraphs = tmp
+ # onsite_term_errors[i] contains an exception if the corresponding
+ # term has a value function with an illegal signature. We only raise
+ # the exception if we actually try to evaluate the offending term
+ # (to maintain backwards compatibility).
+ onsite_term_errors = [
+ err if isinstance(err, Exception) else None
+ for err in onsite_term_parameters
+ ]
+ onsite_terms = [
+ system.Term(
+ subgraph=term_offset + i,
+ hermitian=False,
+ parameters=(
+ params if not isinstance(params, Exception) else None
+ ),
+ )
+ for i, params in enumerate(onsite_term_parameters)
+ ]
+ _onsite_term_by_site_id = np.array(_onsite_term_by_site_id)
+
+ return (onsite_subgraphs, onsite_terms, onsite_term_values,
+ onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id)
+
+
+def _make_hopping_terms(builder, graph, sites, site_offsets, cell_size, term_offset):
+ # Construct the hopping terms
+ #
+ # The logic is the same as for the onsite terms, with the following
+ # differences.
+ #
+ # _hopping_term_by_edge_id
+ # Maps hopping edge IDs to the number of the term that the hopping
+ # is a part of. For Hermitian conjugate hoppings "-term_number -1"
+ # is stored instead.
+
+ hopping_subgraphs = []
+ hopping_term_values = []
+ hopping_term_parameters = []
+ # We just use the cache to handle non/callables and exceptions
+ cache = _value_params_pair_cache(2)
+ # maps hopping key -> hopping term number
+ hopping_to_term_nr = collections.OrderedDict()
+ _hopping_term_by_edge_id = []
+ for tail, head in graph:
+ tail_site, head_site = sites[tail], sites[head]
+ if tail >= cell_size:
+ # The tail belongs to the previous domain. Find the
+ # corresponding hopping with the tail in the fund. domain.
+ tail_site, head_site = builder.symmetry.to_fd(tail_site, head_site)
+ val = builder._get_edge(tail_site, head_site)
+ herm_conj = val is Other
+ if herm_conj:
+ val = builder._get_edge(*builder.symmetry.to_fd(head_site, tail_site))
+ const_val = not callable(val)
+
+ tail_site_array = bisect.bisect(site_offsets, tail) - 1
+ head_site_array = bisect.bisect(site_offsets, head) - 1
+ # "key" uniquely identifies a hopping term.
+ if const_val:
+ key = (None, tail_site_array, head_site_array)
+ else:
+ key = (id(val), tail_site_array, head_site_array)
+ if herm_conj:
+ key = (key[0], head_site_array, tail_site_array)
+
+ if key not in hopping_to_term_nr:
+ # start a new term
+ hopping_to_term_nr[key] = len(hopping_subgraphs)
+ hopping_subgraphs.append([])
+ if const_val:
+ hopping_term_values.append([])
+ hopping_term_parameters.append(None)
+ else:
+ val, params = cache(val)
+ hopping_term_values.append(val)
+ hopping_term_parameters.append(params)
+
+ if herm_conj:
+ # Hopping terms come after all onsite terms, so we need
+ # to offset the term ID by that many
+ term_id = -term_offset - hopping_to_term_nr[key] - 1
+ _hopping_term_by_edge_id.append(term_id)
+ else:
+ # Hopping terms come after all onsite terms, so we need
+ # to offset the term ID by that many
+ term_id = term_offset + hopping_to_term_nr[key]
+ _hopping_term_by_edge_id.append(term_id)
+ hopping_subgraphs[hopping_to_term_nr[key]].append((tail, head))
+ if const_val:
+ vals = hopping_term_values[hopping_to_term_nr[key]]
+ vals.append(val)
+ # Sort the hoppings in each term, and normalize any constant
+ # values to arrays of the correct dtype and shape.
+ if hopping_subgraphs:
+ hopping_subgraphs, hopping_term_values = zip(*(
+ _sort_hopping_term(term, value)
+ for term, value in
+ zip(hopping_subgraphs, hopping_term_values)))
+ # Store site array index and site offsets rather than sites themselves
+ tmp = []
+ for (_, tail_which, head_which), h in zip(hopping_to_term_nr,
+ hopping_subgraphs):
+ start = (site_offsets[tail_which], site_offsets[head_which])
+ # Transpose to get a pair of arrays rather than array of pairs
+ # We use the fact that the underlying array is stored in
+ # array-of-pairs order to search through it in 'hamiltonian'.
+ tmp.append(((tail_which, head_which), (h - start).transpose()))
+ hopping_subgraphs = tmp
+ # hopping_term_errors[i] contains an exception if the corresponding
+ # term has a value function with an illegal signature. We only raise
+ # the exception if this becomes a problem (to maintain backwards
+ # compatibility)
+ hopping_term_errors = [
+ err if isinstance(err, Exception) else None
+ for err in hopping_term_parameters
+ ]
+ hopping_terms = [
+ system.Term(
+ subgraph=term_offset + i,
+ hermitian=True, # Builders are always Hermitian
+ parameters=(
+ params if not isinstance(params, Exception) else None
+ ),
+ )
+ for i, params in enumerate(hopping_term_parameters)
+ ]
+ _hopping_term_by_edge_id = np.array(_hopping_term_by_edge_id)
+
+ return (hopping_subgraphs, hopping_terms, hopping_term_values,
+ hopping_term_parameters, hopping_term_errors,
+ _hopping_term_by_edge_id)
+
+
+class FiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.FiniteVectorizedSystem):
+ """Finalized `Builder` with leads.
+
+ Usable as input for the solvers in `kwant.solvers`.
+
+ Attributes
+ ----------
+ site_arrays : sequence of `~kwant.system.SiteArray`
+ The sites of the system.
+ sites : sequence
+ ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
+ to the integer-labeled site ``i`` of the low-level system. The sites
+ are ordered first by their family and then by their tag.
+ id_by_site : mapping
+ The inverse of ``sites``; maps high-level `~kwant.system.Site`
+ instances to their integer label.
+ Satisfies ``id_by_site[sites[i]] == i``.
+ """
+
+ def __init__(self, builder):
+ assert builder.symmetry.num_directions == 0
+ assert builder.vectorize
+
+ sites = sorted(builder.H)
+ id_by_site = {site: site_id for site_id, site in enumerate(sites)}
+
+ if not all(s.family.norbs for s in sites):
+ raise ValueError("All sites must belong to families with norbs "
+ "specified for vectorized Builders. Specify "
+ "norbs when creating site families.")
+
+ graph = _make_graph(builder.H, id_by_site)
+
+ finalized_leads, lead_interfaces, lead_paddings =\
+ _finalize_leads(builder.leads, id_by_site)
+
+ del id_by_site # cleanup due to large size
+
+ site_arrays = _make_site_arrays(builder.H)
+ # We need this to efficiently find which array a given
+ # site belongs to
+ site_offsets = np.cumsum([0] + [len(arr) for arr in site_arrays])
+
+ (onsite_subgraphs, onsite_terms, onsite_term_values,
+ onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id) =\
+ _make_onsite_terms(builder, sites, site_offsets, term_offset=0)
+
+ (hopping_subgraphs, hopping_terms, hopping_term_values,
+ hopping_term_parameters, hopping_term_errors,
+ _hopping_term_by_edge_id) =\
+ _make_hopping_terms(builder, graph, sites, site_offsets,
+ len(sites), term_offset=len(onsite_terms))
+
+ # Construct the combined onsite/hopping term datastructures
+ subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
+ terms = tuple(onsite_terms) + tuple(hopping_terms)
+ term_values = tuple(onsite_term_values) + tuple(hopping_term_values)
+ term_errors = tuple(onsite_term_errors) + tuple(hopping_term_errors)
+
+ # Construct system parameters
+ parameters = set()
+ for params in chain(onsite_term_parameters, hopping_term_parameters):
+ if params is not None:
+ parameters.update(params)
+ parameters = frozenset(parameters)
+
+ self.site_arrays = site_arrays
+ self.sites = _Sites(self.site_arrays)
+ self.id_by_site = _IdBySite(self.site_arrays)
+ self.graph = graph
+ self.subgraphs = subgraphs
+ self.terms = terms
+ self._term_values = term_values
+ self._term_errors = term_errors
+ self._hopping_term_by_edge_id = _hopping_term_by_edge_id
+ self._onsite_term_by_site_id = _onsite_term_by_site_id
+ self.parameters = parameters
+ self.symmetry = builder.symmetry
+ self.leads = finalized_leads
+ self.lead_interfaces = lead_interfaces
+ self.lead_paddings = lead_paddings
+ self._init_discrete_symmetries(builder)
+
+
+def _make_lead_sites(builder, interface_order):
+ #### For each site of the fundamental domain, determine whether it has
+ #### neighbors in the previous domain or not.
+ sym = builder.symmetry
+ lsites_with = [] # Fund. domain sites with neighbors in prev. dom
+ lsites_without = [] # Remaining sites of the fundamental domain
+ for tail in builder.H: # Loop over all sites of the fund. domain.
+ for head in builder._out_neighbors(tail):
+ fd = sym.which(head)[0]
+ if fd == 1:
+ # Tail belongs to fund. domain, head to the next domain.
+ lsites_with.append(tail)
+ break
+ else:
+ # Tail is a fund. domain site not connected to prev. domain.
+ lsites_without.append(tail)
+
+ if not lsites_with:
+ warnings.warn('Infinite system with disconnected cells.',
+ RuntimeWarning, stacklevel=3)
+
+ ### Create list of sites and a lookup table
+ minus_one = ta.array((-1,))
+ plus_one = ta.array((1,))
+ if interface_order is None:
+ # interface must be sorted
+ interface = [sym.act(minus_one, s) for s in lsites_with]
+ interface.sort()
+ else:
+ lsites_with_set = set(lsites_with)
+ lsites_with = []
+ interface = []
+ if interface_order:
+ shift = ta.array((-sym.which(interface_order[0])[0] - 1,))
+ for shifted_iface_site in interface_order:
+ # Shift the interface domain before the fundamental domain.
+ # That's the right place for the interface of a lead to be, but
+ # the sites of interface_order might live in a different
+ # domain.
+ iface_site = sym.act(shift, shifted_iface_site)
+ lsite = sym.act(plus_one, iface_site)
+
+ try:
+ lsites_with_set.remove(lsite)
+ except KeyError:
+ if (-sym.which(shifted_iface_site)[0] - 1,) != shift:
+ raise ValueError(
+ 'The sites in interface_order do not all '
+ 'belong to the same lead cell.')
+ else:
+ raise ValueError('A site in interface_order is not an '
+ 'interface site:\n' + str(iface_site))
+ interface.append(iface_site)
+ lsites_with.append(lsite)
+ if lsites_with_set:
+ raise ValueError(
+ 'interface_order did not contain all interface sites.')
+ # `interface_order` *must* be sorted, hence `interface` should also
+ if interface != sorted(interface):
+ raise ValueError('Interface sites must be sorted.')
+ del lsites_with_set
+
+ return sorted(lsites_with), sorted(lsites_without), interface
+
+
+def _make_lead_graph(builder, sites, id_by_site, cell_size):
+ sym = builder.symmetry
+ g = graph.Graph()
+ g.num_nodes = len(sites) # Some sites could not appear in any edge.
+ for tail_id, tail in enumerate(sites[:cell_size]):
+ for head in builder._out_neighbors(tail):
+ head_id = id_by_site.get(head)
+ if head_id is None:
+ # Head belongs neither to the fundamental domain nor to the
+ # previous domain. Check that it belongs to the next
+ # domain and ignore it otherwise as an edge corresponding
+ # to this one has been added already or will be added.
+ fd = sym.which(head)[0]
+ if fd != 1:
+ msg = ('Further-than-nearest-neighbor cells '
+ 'are connected by hopping\n{0}.')
+ raise ValueError(msg.format((tail, head)))
+ continue
+ if head_id >= cell_size:
+ # Head belongs to previous domain. The edge added here
+ # correspond to one left out just above.
+ g.add_edge(head_id, tail_id)
+ g.add_edge(tail_id, head_id)
+
+ return g.compressed()
class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
@@ -2275,10 +2819,10 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
Attributes
----------
sites : sequence
- ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
+ ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
to the integer-labeled site ``i`` of the low-level system.
- id_by_site : dict
- The inverse of ``sites``; maps high-level `~kwant.builder.Site`
+ id_by_site : mapping
+ The inverse of ``sites``; maps high-level `~kwant.system.Site`
instances to their integer label.
Satisfies ``id_by_site[sites[i]] == i``.
@@ -2302,69 +2846,12 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
sym = builder.symmetry
assert sym.num_directions == 1
- #### For each site of the fundamental domain, determine whether it has
- #### neighbors in the previous domain or not.
- lsites_with = [] # Fund. domain sites with neighbors in prev. dom
- lsites_without = [] # Remaining sites of the fundamental domain
- for tail in builder.H: # Loop over all sites of the fund. domain.
- for head in builder._out_neighbors(tail):
- fd = sym.which(head)[0]
- if fd == 1:
- # Tail belongs to fund. domain, head to the next domain.
- lsites_with.append(tail)
- break
- else:
- # Tail is a fund. domain site not connected to prev. domain.
- lsites_without.append(tail)
+ lsites_with, lsites_without, interface =\
+ _make_lead_sites(builder, interface_order)
cell_size = len(lsites_with) + len(lsites_without)
- if not lsites_with:
- warnings.warn('Infinite system with disconnected cells.',
- RuntimeWarning, stacklevel=3)
-
- ### Create list of sites and a lookup table
- minus_one = ta.array((-1,))
- plus_one = ta.array((1,))
- if interface_order is None:
- # interface must be sorted
- interface = [sym.act(minus_one, s) for s in lsites_with]
- interface.sort()
- else:
- lsites_with_set = set(lsites_with)
- lsites_with = []
- interface = []
- if interface_order:
- shift = ta.array((-sym.which(interface_order[0])[0] - 1,))
- for shifted_iface_site in interface_order:
- # Shift the interface domain before the fundamental domain.
- # That's the right place for the interface of a lead to be, but
- # the sites of interface_order might live in a different
- # domain.
- iface_site = sym.act(shift, shifted_iface_site)
- lsite = sym.act(plus_one, iface_site)
-
- try:
- lsites_with_set.remove(lsite)
- except KeyError:
- if (-sym.which(shifted_iface_site)[0] - 1,) != shift:
- raise ValueError(
- 'The sites in interface_order do not all '
- 'belong to the same lead cell.')
- else:
- raise ValueError('A site in interface_order is not an '
- 'interface site:\n' + str(iface_site))
- interface.append(iface_site)
- lsites_with.append(lsite)
- if lsites_with_set:
- raise ValueError(
- 'interface_order did not contain all interface sites.')
- # `interface_order` *must* be sorted, hence `interface` should also
- if interface != sorted(interface):
- raise ValueError('Interface sites must be sorted.')
- del lsites_with_set
-
# we previously sorted the interface, so don't sort it again
- sites = sorted(lsites_with) + sorted(lsites_without) + interface
+ sites = lsites_with + lsites_without + interface
del lsites_with
del lsites_without
del interface
@@ -2372,41 +2859,20 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
for site_id, site in enumerate(sites):
id_by_site[site] = site_id
+ graph = _make_lead_graph(builder, sites, id_by_site, cell_size)
+
# In the following, because many onsites/hoppings share the same
# (value, parameter) pairs, we keep them in 'cache' so that we only
# store a given pair in memory *once*. This is like interning strings.
- #### Make graph and extract onsite Hamiltonians.
+ #### Extract onsites
cache = _value_params_pair_cache(1)
- g = graph.Graph()
- g.num_nodes = len(sites) # Some sites could not appear in any edge.
- onsites = []
- for tail_id, tail in enumerate(sites[:cell_size]):
- onsites.append(cache(builder.H[tail][1]))
- for head in builder._out_neighbors(tail):
- head_id = id_by_site.get(head)
- if head_id is None:
- # Head belongs neither to the fundamental domain nor to the
- # previous domain. Check that it belongs to the next
- # domain and ignore it otherwise as an edge corresponding
- # to this one has been added already or will be added.
- fd = sym.which(head)[0]
- if fd != 1:
- msg = ('Further-than-nearest-neighbor cells '
- 'are connected by hopping\n{0}.')
- raise ValueError(msg.format((tail, head)))
- continue
- if head_id >= cell_size:
- # Head belongs to previous domain. The edge added here
- # correspond to one left out just above.
- g.add_edge(head_id, tail_id)
- g.add_edge(tail_id, head_id)
- g = g.compressed()
+ onsites = [cache(builder.H[tail][1]) for tail in sites[:cell_size]]
#### Extract hoppings.
cache = _value_params_pair_cache(2)
hoppings = []
- for tail_id, head_id in g:
+ for tail_id, head_id in graph:
tail = sites[tail_id]
head = sites[head_id]
if tail_id >= cell_size:
@@ -2433,7 +2899,7 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
else:
parameters = frozenset(parameters)
- self.graph = g
+ self.graph = graph
self.sites = sites
self.site_ranges = _site_ranges(sites)
self.id_by_site = id_by_site
@@ -2444,6 +2910,125 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
self.cell_size = cell_size
self._init_discrete_symmetries(builder)
+ def hamiltonian(self, i, j, *args, params=None):
+ cs = self.cell_size
+ if i == j >= cs:
+ i -= cs
+ j -= cs
+ return super().hamiltonian(i, j, *args, params=params)
+
+
+class InfiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.InfiniteVectorizedSystem):
+ """Finalized infinite system, extracted from a `Builder`.
+
+ Attributes
+ ----------
+ site_arrays : sequence of `~kwant.system.SiteArray`
+ The sites of the system.
+ sites : sequence
+ ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
+ to the integer-labeled site ``i`` of the low-level system.
+ id_by_site : mapping
+ The inverse of ``sites``; maps high-level `~kwant.system.Site`
+ instances to their integer label.
+ Satisfies ``id_by_site[sites[i]] == i``.
+
+ Notes
+ -----
+ In infinite systems ``sites`` and ``site_arrays`` consists of 3 parts:
+ sites in the fundamental domain (FD) with hoppings to neighboring cells,
+ sites in the FD with no hoppings to neighboring cells, and sites in FD+1
+ attached to the FD by hoppings. Each of these three subsequences is
+ individually sorted.
+ """
+
+ def __init__(self, builder, interface_order=None):
+ """
+ Finalize a builder instance which has to have exactly a single
+ symmetry direction.
+
+ If interface_order is not set, the order of the interface sites in the
+ finalized system will be arbitrary. If interface_order is set to a
+ sequence of interface sites, this order will be kept.
+ """
+ sym = builder.symmetry
+ assert sym.num_directions == 1
+ assert builder.vectorize
+
+ lsites_with, lsites_without, interface =\
+ _make_lead_sites(builder, interface_order)
+ cell_size = len(lsites_with) + len(lsites_without)
+
+
+ sites = lsites_with + lsites_without + interface
+ id_by_site = {}
+ for site_id, site in enumerate(sites):
+ id_by_site[site] = site_id
+ # these are needed for constructing hoppings
+ lsites_with = frozenset(lsites_with)
+ lsites_without = frozenset(lsites_without)
+ interface = frozenset(interface)
+
+ if not all(s.family.norbs for s in sites):
+ raise ValueError("All sites must belong to families with norbs "
+ "specified for vectorized Builders. Specify "
+ "norbs when creating site families.")
+
+ graph = _make_lead_graph(builder, sites, id_by_site, cell_size)
+
+ del id_by_site # cleanup due to large size
+
+ # In order to conform to the kwant.system.InfiniteVectorizedSystem
+ # interface we need to put the sites that connect to the previous
+ # cell *first*, then the sites that do not couple to the previous
+ # cell, then the sites in the previous cell. Because sites in
+ # a SiteArray are sorted by tag this means that the sites in these
+ # 3 different sets need to be in different SiteArrays.
+ site_arrays = (
+ _make_site_arrays(lsites_with)
+ + _make_site_arrays(lsites_without)
+ + _make_site_arrays(interface)
+ )
+
+ site_offsets = np.cumsum([0] + [len(arr) for arr in site_arrays])
+
+ (onsite_subgraphs, onsite_terms, onsite_term_values,
+ onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id) =\
+ _make_onsite_terms(builder, sites, site_offsets, term_offset=0)
+
+ (hopping_subgraphs, hopping_terms, hopping_term_values,
+ hopping_term_parameters, hopping_term_errors,
+ _hopping_term_by_edge_id) =\
+ _make_hopping_terms(builder, graph, sites, site_offsets,
+ cell_size, term_offset=len(onsite_terms))
+
+ # Construct the combined onsite/hopping term datastructures
+ subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
+ terms = tuple(onsite_terms) + tuple(hopping_terms)
+ term_values = tuple(onsite_term_values) + tuple(hopping_term_values)
+ term_errors = tuple(onsite_term_errors) + tuple(hopping_term_errors)
+
+ # Construct system parameters
+ parameters = set()
+ for params in chain(onsite_term_parameters, hopping_term_parameters):
+ if params is not None:
+ parameters.update(params)
+ parameters = frozenset(parameters)
+
+ self.site_arrays = site_arrays
+ self.sites = _Sites(self.site_arrays)
+ self.id_by_site = _IdBySite(self.site_arrays)
+ self.graph = graph
+ self.subgraphs = subgraphs
+ self.terms = terms
+ self._term_values = term_values
+ self._term_errors = term_errors
+ self._hopping_term_by_edge_id = _hopping_term_by_edge_id
+ self._onsite_term_by_site_id = _onsite_term_by_site_id
+ self.parameters = parameters
+ self.symmetry = builder.symmetry
+ self.cell_size = cell_size
+ self._init_discrete_symmetries(builder)
def hamiltonian(self, i, j, *args, params=None):
cs = self.cell_size
@@ -2452,5 +3037,15 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
j -= cs
return super().hamiltonian(i, j, *args, params=params)
- def pos(self, i):
- return self.sites[i].pos
+
+def is_finite_system(syst):
+ return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
+
+
+def is_infinite_system(syst):
+ return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
+
+
+def is_system(syst):
+ return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem,
+ InfiniteSystem, InfiniteVectorizedSystem))
diff --git a/kwant/continuum/discretizer.py b/kwant/continuum/discretizer.py
index 53f2379cd53f17588dc9e59cb0300a11e90c54c3..ce9e6fe2f65bffeed3da826b02f1618e474a1757 100644
--- a/kwant/continuum/discretizer.py
+++ b/kwant/continuum/discretizer.py
@@ -1,4 +1,4 @@
-# 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
diff --git a/kwant/continuum/landau_levels.py b/kwant/continuum/landau_levels.py
index 13e70c4e6b7224d3dab456d5a4497587514ea8af..c653c55ef2634cd13e932f0c3f9c5fd026c10ff9 100644
--- a/kwant/continuum/landau_levels.py
+++ b/kwant/continuum/landau_levels.py
@@ -176,7 +176,7 @@ class LandauLattice(kwant.lattice.Monatomic):
"""
A `~kwant.lattice.Monatomic` lattice with a Landau level index per site.
- Site tags (see `~kwant.builder.SiteFamily`) are pairs of integers, where
+ 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.
diff --git a/kwant/kpm.py b/kwant/kpm.py
index 5da09b097e394bc57c07c7de97fb958a3665fc3b..44cf9758f4ddfff3f098be074ca929689c562f11 100644
--- a/kwant/kpm.py
+++ b/kwant/kpm.py
@@ -1,5 +1,5 @@
# -*- 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
@@ -928,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``
@@ -950,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``
diff --git a/kwant/lattice.py b/kwant/lattice.py
index 899ce5474216421b97f35b922a51a62cd85b9602..e2557796f127c14ffde5a14ac7e2681a98ba7b60 100644
--- a/kwant/lattice.py
+++ b/kwant/lattice.py
@@ -1,4 +1,4 @@
-# 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.
diff --git a/kwant/operator.pyx b/kwant/operator.pyx
index a10624d1ccf3d7d6a2a175e16d6e00e3126c3bad..2a3e71edd01d492814d97dcf702216833c96fb42 100644
--- a/kwant/operator.pyx
+++ b/kwant/operator.pyx
@@ -13,7 +13,7 @@ import cython
from operator import itemgetter
import functools as ft
import collections
-import numbers
+import warnings
import numpy as np
import tinyarray as ta
@@ -25,8 +25,10 @@ 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
@@ -150,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)
@@ -166,7 +168,7 @@ def _normalize_site_where(syst, where):
otherwise it should contain integers.
"""
if where is None:
- if isinstance(syst, InfiniteSystem):
+ if is_infinite(syst):
where = list(range(syst.cell_size))
else:
where = list(range(syst.graph.num_nodes))
@@ -174,13 +176,12 @@ def _normalize_site_where(syst, where):
try:
where = [syst.id_by_site[s] for s in filter(where, syst.sites)]
except AttributeError:
- if isinstance(syst, InfiniteSystem):
+ 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:
- # Cannot check for builder.Site due to circular imports
- if not isinstance(where[0], numbers.Integral):
+ if isinstance(where[0], Site):
try:
where = [syst.id_by_site[s] for s in where]
except AttributeError:
@@ -189,7 +190,7 @@ def _normalize_site_where(syst, where):
where = np.asarray(where, dtype=gint_dtype).reshape(-1, 1)
- if isinstance(syst, InfiniteSystem) and np.any(where >= syst.cell_size):
+ 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)):
@@ -210,7 +211,7 @@ def _normalize_hopping_where(syst, where):
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)
@@ -223,8 +224,7 @@ def _normalize_hopping_where(syst, where):
else:
where = list(filter(lambda h: where(*h), syst.graph))
else:
- # Cannot check for builder.Site due to circular imports
- if not isinstance(where[0][0], numbers.Integral):
+ if isinstance(where[0][0], Site):
try:
where = list((syst.id_by_site[a], syst.id_by_site[b])
for a, b in where)
@@ -244,7 +244,7 @@ def _normalize_hopping_where(syst, where):
where = np.asarray(where, dtype=gint_dtype)
- if isinstance(syst, InfiniteSystem) and np.any(where > syst.cell_size):
+ if is_infinite(syst) and np.any(where > syst.cell_size):
raise ValueError('Only intra-cell hoppings may be specified '
'using `where`.')
@@ -702,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 (
@@ -735,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
@@ -884,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
@@ -1016,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
diff --git a/kwant/physics/gauge.py b/kwant/physics/gauge.py
index eb4e696164ac9c7c3f9e5ce1c986df601d4d0bb2..cfb927a6fb1e3397ad4d739160f0e4031b736060 100644
--- a/kwant/physics/gauge.py
+++ b/kwant/physics/gauge.py
@@ -1,4 +1,4 @@
-# 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)
diff --git a/kwant/plotter.py b/kwant/plotter.py
index fabe5f61c9234b2323e0c11c860962b3b5b2bd5d..31a0aba6658b18eb3c79fccc3c0ffcdb2804bedf 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -1,5 +1,5 @@
# -*- 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
@@ -812,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.
@@ -961,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))
@@ -1291,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)
@@ -1407,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:
@@ -1483,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):
@@ -1751,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:
@@ -1844,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):
diff --git a/kwant/system.py b/kwant/system.py
index d270691508189866a0d75809917048303443c883..a1c287bb97a982353d79410b8c2cd4fc2ca4bbad 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -1,4 +1,4 @@
-# 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
@@ -8,13 +8,269 @@
"""Low-level interface of systems"""
-__all__ = ['System', 'FiniteSystem', 'InfiniteSystem']
+__all__ = [
+ 'Site', 'SiteArray', 'SiteFamily',
+ 'System', 'VectorizedSystem', 'FiniteSystem', 'FiniteVectorizedSystem',
+ 'InfiniteSystem', 'InfiniteVectorizedSystem',
+]
import abc
import warnings
+import operator
from copy import copy
+from collections import namedtuple
+from functools import total_ordering, lru_cache
+import numpy as np
from . import _system
-from ._common import deprecate_args
+from ._common import deprecate_args, KwantDeprecationWarning
+
+
+
+################ Sites and Site families
+
+class Site(tuple):
+ """A site, member of a `SiteFamily`.
+
+ Sites are the vertices of the graph which describes the tight binding
+ system in a `Builder`.
+
+ A site is uniquely identified by its family and its tag.
+
+ Parameters
+ ----------
+ family : an instance of `SiteFamily`
+ The 'type' of the site.
+ tag : a hashable python object
+ The unique identifier of the site within the site family, typically a
+ vector of integers.
+
+ Raises
+ ------
+ ValueError
+ If `tag` is not a proper tag for `family`.
+
+ Notes
+ -----
+ For convenience, ``family(*tag)`` can be used instead of ``Site(family,
+ tag)`` to create a site.
+
+ The parameters of the constructor (see above) are stored as instance
+ variables under the same names. Given a site ``site``, common things to
+ query are thus ``site.family``, ``site.tag``, and ``site.pos``.
+ """
+ __slots__ = ()
+
+ family = property(operator.itemgetter(0),
+ doc="The site family to which the site belongs.")
+ tag = property(operator.itemgetter(1), doc="The tag of the site.")
+
+
+ def __new__(cls, family, tag, _i_know_what_i_do=False):
+ if _i_know_what_i_do:
+ return tuple.__new__(cls, (family, tag))
+ try:
+ tag = family.normalize_tag(tag)
+ except (TypeError, ValueError) as e:
+ msg = 'Tag {0} is not allowed for site family {1}: {2}'
+ raise type(e)(msg.format(repr(tag), repr(family), e.args[0]))
+ return tuple.__new__(cls, (family, tag))
+
+ def __repr__(self):
+ return 'Site({0}, {1})'.format(repr(self.family), repr(self.tag))
+
+ def __str__(self):
+ sf = self.family
+ return '<Site {0} of {1}>'.format(self.tag, sf.name if sf.name else sf)
+
+ def __getnewargs__(self):
+ return (self.family, self.tag, True)
+
+ @property
+ def pos(self):
+ """Real space position of the site.
+
+ This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
+ """
+ return self.family.pos(self.tag)
+
+
+class SiteArray:
+ """An array of sites, members of a `SiteFamily`.
+
+ Parameters
+ ----------
+ family : an instance of `SiteFamily`
+ The 'type' of the sites.
+ tags : a sequence of python objects
+ Sequence of unique identifiers of the sites within the
+ site array family, typically vectors of integers.
+
+ Raises
+ ------
+ ValueError
+ If `tags` are not proper tags for `family`.
+
+ See Also
+ --------
+ kwant.system.Site
+ """
+
+ def __init__(self, family, tags):
+ self.family = family
+ try:
+ tags = family.normalize_tags(tags)
+ except (TypeError, ValueError) as e:
+ msg = 'Tags {0} are not allowed for site family {1}: {2}'
+ raise type(e)(msg.format(repr(tags), repr(family), e.args[0]))
+ self.tags = tags
+
+ def __repr__(self):
+ return 'SiteArray({0}, {1})'.format(repr(self.family), repr(self.tags))
+
+ def __str__(self):
+ sf = self.family
+ return ('<SiteArray {0} of {1}>'
+ .format(self.tags, sf.name if sf.name else sf))
+
+ def __len__(self):
+ return len(self.tags)
+
+ def __eq__(self, other):
+ if not isinstance(other, SiteArray):
+ raise NotImplementedError()
+ return self.family == other.family and np.all(self.tags == other.tags)
+
+ def positions(self):
+ """Real space position of the site.
+
+ This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
+ """
+ return self.family.positions(self.tags)
+
+
+@total_ordering
+class SiteFamily:
+ """Abstract base class for site families.
+
+ Site families are the 'type' of `Site` objects. Within a family, individual
+ sites are uniquely identified by tags. Valid tags must be hashable Python
+ objects, further details are up to the family.
+
+ Site families must be immutable and fully defined by their initial
+ arguments. They must inherit from this abstract base class and call its
+ __init__ function providing it with two arguments: a canonical
+ representation and a name. The canonical representation will be returned as
+ the objects representation and must uniquely identify the site family
+ instance. The name is a string used to distinguish otherwise identical site
+ families. It may be empty. ``norbs`` defines the number of orbitals
+ on sites associated with this site family; it may be `None`, in which case
+ the number of orbitals is not specified.
+
+
+ All site families must define either 'normalize_tag' or 'normalize_tags',
+ which brings a tag (or, in the latter case, a sequence of tags) to the
+ standard format for this site family.
+
+ Site families may also implement methods ``pos(tag)`` and
+ ``positions(tags)``, which return a vector of realspace coordinates or an
+ array of vectors of realspace coordinates of the site(s) belonging to this
+ family with the given tag(s). These methods are used in plotting routines.
+ ``positions(tags)`` should return an array with shape ``(N, M)`` where
+ ``N`` is the length of ``tags``, and ``M`` is the realspace dimension.
+
+ If the ``norbs`` of a site family are provided, and sites of this family
+ are used to populate a `~kwant.builder.Builder`, then the associated
+ Hamiltonian values must have the correct shape. That is, if a site family
+ has ``norbs = 2``, then any on-site terms for sites belonging to this
+ family should be 2x2 matrices. Similarly, any hoppings to/from sites
+ belonging to this family must have a matrix structure where there are two
+ rows/columns. This condition applies equally to Hamiltonian values that
+ are given by functions. If this condition is not satisfied, an error will
+ be raised.
+ """
+
+ def __init__(self, canonical_repr, name, norbs):
+ self.canonical_repr = canonical_repr
+ self.hash = hash(canonical_repr)
+ self.name = name
+ if norbs is None:
+ warnings.warn("Not specfying norbs is deprecated. Always specify "
+ "norbs when creating site families.",
+ KwantDeprecationWarning, stacklevel=3)
+ if norbs is not None:
+ if int(norbs) != norbs or norbs <= 0:
+ raise ValueError('The norbs parameter must be an integer > 0.')
+ norbs = int(norbs)
+ self.norbs = norbs
+
+ def __init_subclass__(cls, **kwargs):
+ super().__init_subclass__(**kwargs)
+ if (cls.normalize_tag is SiteFamily.normalize_tag
+ and cls.normalize_tags is SiteFamily.normalize_tags):
+ raise TypeError("Must redefine either 'normalize_tag' or "
+ "'normalize_tags'")
+
+ def __repr__(self):
+ return self.canonical_repr
+
+ def __str__(self):
+ if self.name:
+ msg = '<{0} site family {1}{2}>'
+ else:
+ msg = '<unnamed {0} site family{2}>'
+ orbs = ' with {0} orbitals'.format(self.norbs) if self.norbs else ''
+ return msg.format(self.__class__.__name__, self.name, orbs)
+
+ def __hash__(self):
+ return self.hash
+
+ def __eq__(self, other):
+ try:
+ return self.canonical_repr == other.canonical_repr
+ except AttributeError:
+ return False
+
+ def __ne__(self, other):
+ try:
+ return self.canonical_repr != other.canonical_repr
+ except AttributeError:
+ return True
+
+ def __lt__(self, other):
+ # If this raises an AttributeError, we were trying
+ # to compare it to something non-comparable anyway.
+ return self.canonical_repr < other.canonical_repr
+
+ def normalize_tag(self, tag):
+ """Return a normalized version of the tag.
+
+ Raises TypeError or ValueError if the tag is not acceptable.
+ """
+ tag, = self.normalize_tags([tag])
+ return tag
+
+ def normalize_tags(self, tags):
+ """Return a normalized version of the tags.
+
+ Raises TypeError or ValueError if the tags are not acceptable.
+ """
+ return np.array([self.normalize_tag(tag) for tag in tags])
+
+ def __call__(self, *tag):
+ """
+ A convenience function.
+
+ This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).
+ """
+ # Catch a likely and difficult to find mistake.
+ if tag and isinstance(tag[0], tuple):
+ raise ValueError('Use site_family(1, 2) instead of '
+ 'site_family((1, 2))!')
+ return Site(self, tag)
+
+
+
+################ Systems
class System(metaclass=abc.ABCMeta):
@@ -92,12 +348,100 @@ class System(metaclass=abc.ABCMeta):
details = ', and '.join((', '.join(details[:-1]), details[-1]))
return '<{} with {}>'.format(self.__class__.__name__, details)
+ hamiltonian_submatrix = _system.hamiltonian_submatrix
+
+
+Term = namedtuple(
+ "Term",
+ ["subgraph", "hermitian", "parameters"],
+)
+
+
+class VectorizedSystem(System, metaclass=abc.ABCMeta):
+ """Abstract general low-level system with support for vectorization.
+
+ Attributes
+ ----------
+ graph : kwant.graph.CGraph
+ The system graph.
+ subgraphs : sequence of tuples
+ Each subgraph has the form '((idx1, idx2), (offsets1, offsets2))'
+ where 'offsets1' and 'offsets2' index sites within the site arrays
+ indexed by 'idx1' and 'idx2'.
+ terms : sequence of tuples
+ Each tuple has the following structure:
+ (subgraph: int, hermitian: bool, parameters: List(str))
+ 'subgraph' indexes 'subgraphs' and supplies the to/from sites of this
+ term. 'hermitian' is 'True' if the term needs its Hermitian
+ conjugate to be added when evaluating the Hamiltonian, and 'parameters'
+ contains a list of parameter names used when evaluating this term.
+ site_arrays : sequence of SiteArray
+ The sites of the system.
+ site_ranges : None or Nx3 integer array
+ Has 1 row per site array, plus one extra row. Each row consists
+ of ``(first_site, norbs, orb_offset)``: the index of the first
+ site in the site array, the number of orbitals on each site in
+ the site array, and the offset of the first orbital of the first
+ site in the site array. In addition, the final row has the form
+ ``(len(graph.num_nodes), 0, tot_norbs)`` where ``tot_norbs`` is the
+ total number of orbitals in the system. ``None`` if any site array
+ in 'site_arrays' does not have 'norbs' specified. Note 'site_ranges'
+ is directly computable from 'site_arrays'.
+ parameters : frozenset of strings
+ The names of the parameters on which the system depends. This attribute
+ is provisional and may be changed in a future version of Kwant
+
+ Notes
+ -----
+ The sites of the system are indexed by integers ranging from 0 to
+ ``self.graph.num_nodes - 1``.
+
+ Optionally, a class derived from ``System`` can provide a method ``pos`` which
+ is assumed to return the real-space position of a site given its index.
+ """
+ @abc.abstractmethod
+ def hamiltonian_term(self, term_number, selector=slice(None),
+ args=(), params=None):
+ """Return the Hamiltonians for hamiltonian term number k.
-# Add a C-implemented function as an unbound method to class System.
-System.hamiltonian_submatrix = _system.hamiltonian_submatrix
+ Parameters
+ ----------
+ term_number : int
+ The number of the term to evaluate.
+ selector : slice or sequence of int, default: slice(None)
+ The elements of the term to evaluate.
+ args : tuple
+ Positional arguments to the term. (Deprecated)
+ params : dict
+ Keyword parameters to the term
+ Returns
+ -------
+ hamiltonian : 3d complex array
+ Has shape ``(N, P, Q)`` where ``N`` is the number of matrix
+ elements in this term (or the number selected by 'selector'
+ if provided), ``P`` and ``Q`` are the number of orbitals in the
+ 'to' and 'from' site arrays associated with this term.
-class FiniteSystem(System, metaclass=abc.ABCMeta):
+ Providing positional arguments via 'args' is deprecated,
+ instead, provide named parameters as a dictionary via 'params'.
+ """
+ @property
+ @lru_cache(1)
+ def site_ranges(self):
+ site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
+ norbs = [arr.family.norbs for arr in self.site_arrays] + [0]
+ if any(norb is None for norb in norbs):
+ return None
+ orb_offsets = np.cumsum(
+ [0] + [len(arr) * arr.family.norbs for arr in self.site_arrays]
+ )
+ return np.array([site_offsets, norbs, orb_offsets]).transpose()
+
+ hamiltonian_submatrix = _system.vectorized_hamiltonian_submatrix
+
+
+class FiniteSystemMixin(metaclass=abc.ABCMeta):
"""Abstract finite low-level system, possibly with leads.
Attributes
@@ -220,7 +564,19 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
return symmetries.validate(ham)
-class InfiniteSystem(System, metaclass=abc.ABCMeta):
+class FiniteSystem(System, FiniteSystemMixin, metaclass=abc.ABCMeta):
+ pass
+
+
+class FiniteVectorizedSystem(VectorizedSystem, FiniteSystemMixin, metaclass=abc.ABCMeta):
+ pass
+
+
+def is_finite(syst):
+ return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
+
+
+class InfiniteSystemMixin(metaclass=abc.ABCMeta):
"""Abstract infinite low-level system.
An infinite system consists of an infinite series of identical cells.
@@ -261,30 +617,10 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
infinite system. The other scheme has the numbers of site 0 and 1
exchanged, as well as of site 3 and 4.
+ Sites in the fundamental domain cell must belong to a different site array
+ than the sites in the previous cell. In the above example this means that
+ sites '(0, 1, 2)' and '(3, 4)' must belong to different site arrays.
"""
- @deprecate_args
- def cell_hamiltonian(self, args=(), sparse=False, *, params=None):
- """Hamiltonian of a single cell of the infinite system.
-
- Providing positional arguments via 'args' is deprecated,
- instead, provide named parameters as a dictionary via 'params'.
- """
- cell_sites = range(self.cell_size)
- return self.hamiltonian_submatrix(args, cell_sites, cell_sites,
- sparse=sparse, params=params)
-
- @deprecate_args
- def inter_cell_hopping(self, args=(), sparse=False, *, params=None):
- """Hopping Hamiltonian between two cells of the infinite system.
-
- Providing positional arguments via 'args' is deprecated,
- instead, provide named parameters as a dictionary via 'params'.
- """
- cell_sites = range(self.cell_size)
- interface_sites = range(self.cell_size, self.graph.num_nodes)
- return self.hamiltonian_submatrix(args, cell_sites, interface_sites,
- sparse=sparse, params=params)
-
@deprecate_args
def modes(self, energy=0, args=(), *, params=None):
"""Return mode decomposition of the lead
@@ -368,6 +704,41 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
return list(broken)
+class InfiniteSystem(System, InfiniteSystemMixin, metaclass=abc.ABCMeta):
+
+ @deprecate_args
+ def cell_hamiltonian(self, args=(), sparse=False, *, params=None):
+ """Hamiltonian of a single cell of the infinite system.
+
+ Providing positional arguments via 'args' is deprecated,
+ instead, provide named parameters as a dictionary via 'params'.
+ """
+ cell_sites = range(self.cell_size)
+ return self.hamiltonian_submatrix(args, cell_sites, cell_sites,
+ sparse=sparse, params=params)
+
+ @deprecate_args
+ def inter_cell_hopping(self, args=(), sparse=False, *, params=None):
+ """Hopping Hamiltonian between two cells of the infinite system.
+
+ Providing positional arguments via 'args' is deprecated,
+ instead, provide named parameters as a dictionary via 'params'.
+ """
+ cell_sites = range(self.cell_size)
+ interface_sites = range(self.cell_size, self.graph.num_nodes)
+ return self.hamiltonian_submatrix(args, cell_sites, interface_sites,
+ sparse=sparse, params=params)
+
+
+class InfiniteVectorizedSystem(VectorizedSystem, InfiniteSystemMixin, metaclass=abc.ABCMeta):
+ cell_hamiltonian = _system.vectorized_cell_hamiltonian
+ inter_cell_hopping = _system.vectorized_inter_cell_hopping
+
+
+def is_infinite(syst):
+ return isinstance(syst, (InfiniteSystem, InfiniteVectorizedSystem))
+
+
class PrecalculatedLead:
def __init__(self, modes=None, selfenergy=None):
"""A general lead defined by its self energy.
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index 7ba00b403513ffe3100acd184820adf1c788ed15..253c9d7cca43251dce2c608a50150b9fb2da7f4f 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -1,4 +1,4 @@
-# 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
@@ -19,8 +19,8 @@ from pytest import raises, warns
from numpy.testing import assert_almost_equal
import kwant
-from kwant import builder
-from kwant._common import ensure_rng
+from kwant import builder, system
+from kwant._common import ensure_rng, KwantDeprecationWarning
def test_bad_keys():
@@ -290,12 +290,28 @@ def random_hopping_integral(rng):
def check_onsite(fsyst, sites, subset=False, check_values=True):
+ vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))
+
+ if vectorized:
+ site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])
+
freq = {}
for node in range(fsyst.graph.num_nodes):
site = fsyst.sites[node].tag
freq[site] = freq.get(site, 0) + 1
if check_values and site in sites:
- assert fsyst.onsites[node][0] is sites[site]
+ if vectorized:
+ term = fsyst._onsite_term_by_site_id[node]
+ value = fsyst._term_values[term]
+ if callable(value):
+ assert value is sites[site]
+ else:
+ (w, _), (off, _) = fsyst.subgraphs[fsyst.terms[term].subgraph]
+ node_off = node - site_offsets[w]
+ selector = np.searchsorted(off, node_off)
+ assert np.allclose(value[selector], sites[site])
+ else:
+ assert fsyst.onsites[node][0] is sites[site]
if not subset:
# Check that all sites of `fsyst` are in `sites`.
for site in freq.keys():
@@ -306,24 +322,50 @@ def check_onsite(fsyst, sites, subset=False, check_values=True):
def check_hoppings(fsyst, hops):
+ vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))
+
+ if vectorized:
+ site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])
+
assert fsyst.graph.num_edges == 2 * len(hops)
for edge_id, edge in enumerate(fsyst.graph):
- tail, head = edge
- tail = fsyst.sites[tail].tag
- head = fsyst.sites[head].tag
- value = fsyst.hoppings[edge_id][0]
- if value is builder.Other:
- assert (head, tail) in hops
+ i, j = edge
+ tail = fsyst.sites[i].tag
+ head = fsyst.sites[j].tag
+
+ if vectorized:
+ term = fsyst._hopping_term_by_edge_id[edge_id]
+ if term < 0: # Hermitian conjugate
+ assert (head, tail) in hops
+ else:
+ value = fsyst._term_values[term]
+ assert (tail, head) in hops
+ if callable(value):
+ assert value is hops[tail, head]
+ else:
+ dtype = np.dtype([('f0', int), ('f1', int)])
+ subgraph = fsyst.terms[term].subgraph
+ (to_w, from_w), hoppings = fsyst.subgraphs[subgraph]
+ hop = (i - site_offsets[to_w], j - site_offsets[from_w])
+ hop = np.array(hop, dtype=dtype)
+ hoppings = hoppings.transpose().view(dtype).reshape(-1)
+ selector = np.recarray.searchsorted(hoppings, hop)
+ assert np.allclose(value[selector], hops[tail, head])
else:
- assert (tail, head) in hops
- assert value is hops[tail, head]
+ value = fsyst.hoppings[edge_id][0]
+ if value is builder.Other:
+ assert (head, tail) in hops
+ else:
+ assert (tail, head) in hops
+ assert value is hops[tail, head]
def check_id_by_site(fsyst):
for i, site in enumerate(fsyst.sites):
assert fsyst.id_by_site[site] == i
-def test_finalization():
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_finalization(vectorize):
"""Test the finalization of finite and infinite systems.
In order to exactly verify the finalization, low-level features of the
@@ -377,7 +419,7 @@ def test_finalization():
neighbors = sorted(neighbors)
# Build scattering region from blueprint and test it.
- syst = builder.Builder()
+ syst = builder.Builder(vectorize=vectorize)
fam = kwant.lattice.general(ta.identity(2), norbs=1)
for site, value in sr_sites.items():
syst[fam(*site)] = value
@@ -388,7 +430,7 @@ def test_finalization():
check_onsite(fsyst, sr_sites)
check_hoppings(fsyst, sr_hops)
# check that sites are sorted
- assert fsyst.sites == tuple(sorted(fam(*site) for site in sr_sites))
+ assert tuple(fsyst.sites) == tuple(sorted(fam(*site) for site in sr_sites))
# Build lead from blueprint and test it.
lead = builder.Builder(kwant.TranslationalSymmetry((size, 0)))
@@ -421,12 +463,12 @@ def test_finalization():
# Attach lead with improper interface.
syst.leads[-1] = builder.BuilderLead(
- lead, 2 * tuple(builder.Site(fam, n) for n in neighbors))
+ lead, 2 * tuple(system.Site(fam, n) for n in neighbors))
raises(ValueError, syst.finalized)
# Attach lead properly.
syst.leads[-1] = builder.BuilderLead(
- lead, (builder.Site(fam, n) for n in neighbors))
+ lead, (system.Site(fam, n) for n in neighbors))
fsyst = syst.finalized()
assert len(fsyst.lead_interfaces) == 1
assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] ==
@@ -434,15 +476,15 @@ def test_finalization():
# test that we cannot finalize a system with a badly sorted interface order
raises(ValueError, builder.InfiniteSystem, lead,
- [builder.Site(fam, n) for n in reversed(neighbors)])
+ [system.Site(fam, n) for n in reversed(neighbors)])
# site ordering independent of whether interface was specified
- flead_order = builder.InfiniteSystem(lead, [builder.Site(fam, n)
+ flead_order = builder.InfiniteSystem(lead, [system.Site(fam, n)
for n in neighbors])
- assert flead.sites == flead_order.sites
+ assert tuple(flead.sites) == tuple(flead_order.sites)
syst.leads[-1] = builder.BuilderLead(
- lead, (builder.Site(fam, n) for n in neighbors))
+ lead, (system.Site(fam, n) for n in neighbors))
fsyst = syst.finalized()
assert len(fsyst.lead_interfaces) == 1
assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] ==
@@ -457,47 +499,62 @@ def test_finalization():
raises(ValueError, lead.finalized)
-def test_site_ranges():
+def _make_system_from_sites(sites, vectorize):
+ syst = builder.Builder(vectorize=vectorize)
+ for s in sites:
+ norbs = s.family.norbs
+ if norbs:
+ syst[s] = np.eye(s.family.norbs)
+ else:
+ syst[s] = None
+ return syst.finalized()
+
+
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_site_ranges(vectorize):
lat1a = kwant.lattice.chain(norbs=1, name='a')
lat1b = kwant.lattice.chain(norbs=1, name='b')
lat2 = kwant.lattice.chain(norbs=2)
- site_ranges = builder._site_ranges
# simple case -- single site family
for lat in (lat1a, lat2):
sites = list(map(lat, range(10)))
- ranges = site_ranges(sites)
+ syst = _make_system_from_sites(sites, vectorize)
+ ranges = syst.site_ranges
expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)]
- assert ranges == expected
+ assert np.array_equal(ranges, expected)
# pair of site families
- sites = it.chain(map(lat1a, range(4)), map(lat1b, range(6)),
- map(lat1a, range(4)))
- expected = [(0, 1, 0), (4, 1, 4), (10, 1, 10), (14, 0, 14)]
- assert expected == site_ranges(tuple(sites))
+ sites = it.chain(map(lat1a, range(4)), map(lat1b, range(6)))
+ syst = _make_system_from_sites(sites, vectorize)
+ expected = [(0, 1, 0), (4, 1, 4), (10, 0, 10)]
+ assert np.array_equal(expected, syst.site_ranges)
sites = it.chain(map(lat2, range(4)), map(lat1a, range(6)),
map(lat1b, range(4)))
+ syst = _make_system_from_sites(sites, vectorize)
expected = [(0, 2, 0), (4, 1, 4*2), (10, 1, 4*2+6), (14, 0, 4*2+10)]
- assert expected == site_ranges(tuple(sites))
+ assert np.array_equal(expected, syst.site_ranges)
# test with an actual builder
for lat in (lat1a, lat2):
sites = list(map(lat, range(10)))
- syst = kwant.Builder()
+ syst = kwant.Builder(vectorize=vectorize)
syst[sites] = np.eye(lat.norbs)
ranges = syst.finalized().site_ranges
expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)]
- assert ranges == expected
- # poison system with a single site with no norbs defined.
- # Also catch the deprecation warning.
- with warnings.catch_warnings():
- warnings.simplefilter("ignore")
- syst[kwant.lattice.chain(norbs=None)(0)] = 1
- ranges = syst.finalized().site_ranges
- assert ranges == None
-
-
-def test_hamiltonian_evaluation():
+ assert np.array_equal(ranges, expected)
+ if not vectorize: # vectorized systems *must* have all norbs
+ # poison system with a single site with no norbs defined.
+ # Also catch the deprecation warning.
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore")
+ syst[kwant.lattice.chain(norbs=None)(0)] = 1
+ ranges = syst.finalized().site_ranges
+ assert ranges is None
+
+
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_hamiltonian_evaluation(vectorize):
def f_onsite(site):
return site.tag[0]
@@ -505,14 +562,26 @@ def test_hamiltonian_evaluation():
a, b = a.tag, b.tag
return complex(a[0] + b[0], a[1] - b[1])
+ def f_onsite_vectorized(sites):
+ return sites.tags[:, 0]
+
+ def f_hopping_vectorized(a, b):
+ a, b = a.tags, b.tags
+ return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1])
+
tags = [(0, 0), (1, 1), (2, 2), (3, 3)]
edges = [(0, 1), (0, 2), (0, 3), (1, 2)]
- syst = builder.Builder()
+ syst = builder.Builder(vectorize=vectorize)
fam = builder.SimpleSiteFamily(norbs=1)
sites = [fam(*tag) for tag in tags]
- syst[(fam(*tag) for tag in tags)] = f_onsite
- syst[((fam(*tags[i]), fam(*tags[j])) for (i, j) in edges)] = f_hopping
+ hoppings = [(sites[i], sites[j]) for i, j in edges]
+ if vectorize:
+ syst[sites] = f_onsite_vectorized
+ syst[hoppings] = f_hopping_vectorized
+ else:
+ syst[sites] = f_onsite
+ syst[hoppings] = f_hopping
fsyst = syst.finalized()
assert fsyst.graph.num_nodes == len(tags)
@@ -521,12 +590,16 @@ def test_hamiltonian_evaluation():
for i in range(len(tags)):
site = fsyst.sites[i]
assert site in sites
- assert fsyst.hamiltonian(i, i) == syst[site](site)
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ assert fsyst.hamiltonian(i, i) == f_onsite(site)
for t, h in fsyst.graph:
tsite = fsyst.sites[t]
hsite = fsyst.sites[h]
- assert fsyst.hamiltonian(t, h) == syst[tsite, hsite](tsite, hsite)
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ assert fsyst.hamiltonian(t, h) == f_hopping(tsite, hsite)
# test when user-function raises errors
def onsite_raises(site):
@@ -539,13 +612,18 @@ def test_hamiltonian_evaluation():
a, b = hop
# exceptions are converted to kwant.UserCodeError and we add our message
with raises(kwant.UserCodeError) as exc_info:
- fsyst.hamiltonian(a, a)
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ fsyst.hamiltonian(a, a)
msg = 'Error occurred in user-supplied value function "onsite_raises"'
assert msg in str(exc_info.value)
for hop in [(a, b), (b, a)]:
with raises(kwant.UserCodeError) as exc_info:
- fsyst.hamiltonian(*hop)
+ with warnings.catch_warnings():
+ # Ignore deprecation warnings for 'hamiltonian'
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ fsyst.hamiltonian(*hop)
msg = ('Error occurred in user-supplied '
'value function "hopping_raises"')
assert msg in str(exc_info.value)
@@ -555,7 +633,7 @@ def test_hamiltonian_evaluation():
syst[new_hop[0]] = onsite_raises
syst[new_hop] = hopping_raises
fsyst = syst.finalized()
- hop = tuple(map(fsyst.sites.index, new_hop))
+ hop = tuple(map(fsyst.id_by_site.__getitem__, new_hop))
test_raising(fsyst, hop)
# test with infinite system
@@ -563,10 +641,118 @@ def test_hamiltonian_evaluation():
for k, v in it.chain(syst.site_value_pairs(), syst.hopping_value_pairs()):
inf_syst[k] = v
inf_fsyst = inf_syst.finalized()
- hop = tuple(map(inf_fsyst.sites.index, new_hop))
+ hop = tuple(map(inf_fsyst.id_by_site.__getitem__, new_hop))
test_raising(inf_fsyst, hop)
+def test_vectorized_hamiltonian_evaluation():
+
+ def onsite(site):
+ return site.tag[0]
+
+ def vectorized_onsite(sites):
+ return sites.tags[:, 0]
+
+ def hopping(to_site, from_site):
+ a, b = to_site.tag, from_site.tag
+ return a[0] + b[0] + 1j * (a[1] - b[1])
+
+ def vectorized_hopping(to_sites, from_sites):
+ a, b = to_sites.tags, from_sites.tags
+ return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1])
+
+ tags = [(0, 0), (1, 1), (2, 2), (3, 3)]
+ edges = [(0, 1), (0, 2), (0, 3), (1, 2)]
+
+ fam = builder.SimpleSiteFamily(norbs=1)
+ sites = [fam(*tag) for tag in tags]
+ hops = [(fam(*tags[i]), fam(*tags[j])) for (i, j) in edges]
+
+ syst_simple = builder.Builder(vectorize=False)
+ syst_simple[sites] = onsite
+ syst_simple[hops] = hopping
+ fsyst_simple = syst_simple.finalized()
+
+ syst_vectorized = builder.Builder(vectorize=True)
+ syst_vectorized[sites] = vectorized_onsite
+ syst_vectorized[hops] = vectorized_hopping
+ fsyst_vectorized = syst_vectorized.finalized()
+
+ assert fsyst_vectorized.graph.num_nodes == len(tags)
+ assert fsyst_vectorized.graph.num_edges == 2 * len(edges)
+ assert len(fsyst_vectorized.site_arrays) == 1
+ assert fsyst_vectorized.site_arrays[0] == system.SiteArray(fam, tags)
+
+ assert np.allclose(
+ fsyst_simple.hamiltonian_submatrix(),
+ fsyst_vectorized.hamiltonian_submatrix(),
+ )
+
+ for i in range(len(tags)):
+ site = fsyst_vectorized.sites[i]
+ assert site in sites
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ assert (
+ fsyst_vectorized.hamiltonian(i, i)
+ == fsyst_simple.hamiltonian(i, i))
+
+ for t, h in fsyst_vectorized.graph:
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ assert (
+ fsyst_vectorized.hamiltonian(t, h)
+ == fsyst_simple.hamiltonian(t, h))
+
+ # Test infinite system, including hoppings that go both ways into
+ # the next cell
+ lat = kwant.lattice.square(norbs=1)
+
+ syst_vectorized = builder.Builder(kwant.TranslationalSymmetry((-1, 0)),
+ vectorize=True)
+ syst_vectorized[lat(0, 0)] = 4
+ syst_vectorized[lat(0, 1)] = 5
+ syst_vectorized[lat(0, 2)] = vectorized_onsite
+ syst_vectorized[(lat(1, 0), lat(0, 0))] = 1j
+ syst_vectorized[(lat(2, 1), lat(1, 1))] = vectorized_hopping
+ fsyst_vectorized = syst_vectorized.finalized()
+
+ syst_simple = builder.Builder(kwant.TranslationalSymmetry((-1, 0)),
+ vectorize=False)
+ syst_simple[lat(0, 0)] = 4
+ syst_simple[lat(0, 1)] = 5
+ syst_simple[lat(0, 2)] = onsite
+ syst_simple[(lat(1, 0), lat(0, 0))] = 1j
+ syst_simple[(lat(2, 1), lat(1, 1))] = hopping
+ fsyst_simple = syst_simple.finalized()
+
+ assert np.allclose(
+ fsyst_vectorized.hamiltonian_submatrix(),
+ fsyst_simple.hamiltonian_submatrix(),
+ )
+ assert np.allclose(
+ fsyst_vectorized.cell_hamiltonian(),
+ fsyst_simple.cell_hamiltonian(),
+ )
+ assert np.allclose(
+ fsyst_vectorized.inter_cell_hopping(),
+ fsyst_simple.inter_cell_hopping(),
+ )
+
+
+def test_vectorized_requires_norbs():
+
+ # Catch deprecation warning for lack of norbs
+ with warnings.catch_warnings():
+ warnings.simplefilter("ignore")
+ fam = builder.SimpleSiteFamily()
+
+ syst = builder.Builder(vectorize=True)
+ syst[fam(0, 0)] = 1
+
+ raises(ValueError, syst.finalized)
+
+
def test_dangling():
def make_system():
# 1
@@ -1204,15 +1390,22 @@ def test_discrete_symmetries():
# We need to keep testing 'args', but we don't want to see
# all the deprecation warnings in the test logs
@pytest.mark.filterwarnings("ignore:.*'args' parameter")
-def test_argument_passing():
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_argument_passing(vectorize):
chain = kwant.lattice.chain(norbs=1)
# Test for passing parameters to hamiltonian matrix elements
def onsite(site, p1, p2):
- return p1 + p2
+ if vectorize:
+ return (p1 + p2) * np.ones(len(site))
+ else:
+ return p1 + p2
def hopping(site1, site2, p1, p2):
- return p1 - p2
+ if vectorize:
+ return (p1 - p2) * np.ones(len(site1))
+ else:
+ return p1 - p2
def gen_fill_syst(onsite, hopping, syst):
syst[(chain(i) for i in range(3))] = onsite
@@ -1221,8 +1414,9 @@ def test_argument_passing():
fill_syst = ft.partial(gen_fill_syst, onsite, hopping)
- syst = fill_syst(kwant.Builder())
- inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,))))
+ syst = fill_syst(kwant.Builder(vectorize=vectorize))
+ inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,)),
+ vectorize=vectorize))
tests = (
syst.hamiltonian_submatrix,
@@ -1238,28 +1432,34 @@ def test_argument_passing():
# test that mixing 'args' and 'params' raises TypeError
with raises(TypeError):
- syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1))
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1))
+
with raises(TypeError):
- inf_syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1))
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ inf_syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1))
# Missing parameters raises TypeError
with raises(TypeError):
- syst.hamiltonian(0, 0, params=dict(p1=2))
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ syst.hamiltonian(0, 0, params=dict(p1=2))
+
with raises(TypeError):
- syst.hamiltonian_submatrix(params=dict(p1=2))
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
+ syst.hamiltonian_submatrix(params=dict(p1=2))
# test that passing parameters without default values works, and that
# passing parameters with default values fails
- def onsite(site, p1, p2):
- return p1 + p2
-
- def hopping(site, site2, p1, p2):
- return p1 - p2
fill_syst = ft.partial(gen_fill_syst, onsite, hopping)
- syst = fill_syst(kwant.Builder())
- inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,))))
+ syst = fill_syst(kwant.Builder(vectorize=vectorize))
+ inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,)),
+ vectorize=vectorize))
tests = (
syst.hamiltonian_submatrix,
@@ -1275,12 +1475,18 @@ def test_argument_passing():
# Some common, some different args for value functions
def onsite2(site, a, b):
- return site.pos + a + b
+ if vectorize:
+ return site.positions()[:, 0] + a + b
+ else:
+ return site.pos[0] + a + b
def hopping2(site1, site2, a, c, b):
- return a + b + c
+ if vectorize:
+ return (a + b + c) * np.ones(len(site1))
+ else:
+ return a + b + c
- syst = kwant.Builder()
+ syst = kwant.Builder(vectorize=vectorize)
syst[(chain(i) for i in range(3))] = onsite2
syst[((chain(i), chain(i + 1)) for i in range(2))] = hopping2
fsyst = syst.finalized()
@@ -1391,6 +1597,7 @@ def test_subs():
lead = lead.finalized()
assert lead.parameters == {'lead_a', 'lead_b', 'lead_c'}
+
def test_attach_stores_padding():
lat = kwant.lattice.chain(norbs=1)
syst = kwant.Builder()
diff --git a/kwant/tests/test_operator.py b/kwant/tests/test_operator.py
index bb8a185e6675356d90552ee3d21b06d241fda83c..c5111b75d88033a74654eb5aa8e0661ee2dc51e9 100644
--- a/kwant/tests/test_operator.py
+++ b/kwant/tests/test_operator.py
@@ -445,7 +445,7 @@ def test_arg_passing(A):
lat1 = kwant.lattice.chain(norbs=1)
syst = kwant.Builder()
- syst[lat1(0)] = syst[lat1(1)] = lambda s0, a, b: s0.pos + a + b
+ syst[lat1(0)] = syst[lat1(1)] = lambda s0, a, b: s0.pos[0] + a + b
syst[lat1.neighbors()] = lambda s0, s1, a, b: a - b
fsyst = syst.finalized()
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index 9bba152bc881013cd57bc949f46fca3244aa4153..3a4034776ee8492221910d0b00dd5a8e1aa97cf1 100644
--- a/kwant/tests/test_system.py
+++ b/kwant/tests/test_system.py
@@ -1,4 +1,4 @@
-# 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
@@ -8,6 +8,7 @@
import pickle
import copy
+import pytest
from pytest import raises
import numpy as np
from scipy import sparse
@@ -15,9 +16,11 @@ import kwant
from kwant._common import ensure_rng
-def test_hamiltonian_submatrix():
- syst = kwant.Builder()
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_hamiltonian_submatrix(vectorize):
+ syst = kwant.Builder(vectorize=vectorize)
chain = kwant.lattice.chain(norbs=1)
+ chain2 = kwant.lattice.chain(norbs=2)
for i in range(3):
syst[chain(i)] = 0.5 * i
for i in range(2):
@@ -27,7 +30,11 @@ def test_hamiltonian_submatrix():
mat = syst2.hamiltonian_submatrix()
assert mat.shape == (3, 3)
# Sorting is required due to unknown compression order of builder.
- perm = np.argsort([os[0] for os in syst2.onsites])
+ if vectorize:
+ _, (site_offsets, _) = syst2.subgraphs[0]
+ else:
+ site_offsets = [os[0] for os in syst2.onsites]
+ perm = np.argsort(site_offsets)
mat_should_be = np.array([[0, 1j, 0], [-1j, 0.5, 2j], [0, -2j, 1]])
mat = mat[perm, :]
@@ -41,19 +48,12 @@ def test_hamiltonian_submatrix():
mat = mat[:, perm]
np.testing.assert_array_equal(mat, mat_should_be)
- mat = syst2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]])
- np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3])
-
- mat = syst2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]], sparse=True)
- mat = mat.toarray()
- np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3])
-
# Test for correct treatment of matrix input.
- syst = kwant.Builder()
- syst[chain(0)] = np.array([[0, 1j], [-1j, 0]])
+ syst = kwant.Builder(vectorize=vectorize)
+ syst[chain2(0)] = np.array([[0, 1j], [-1j, 0]])
syst[chain(1)] = np.array([[1]])
syst[chain(2)] = np.array([[2]])
- syst[chain(1), chain(0)] = np.array([[1, 2j]])
+ syst[chain(1), chain2(0)] = np.array([[1, 2j]])
syst[chain(2), chain(1)] = np.array([[3j]])
syst2 = syst.finalized()
mat_dense = syst2.hamiltonian_submatrix()
@@ -62,9 +62,10 @@ def test_hamiltonian_submatrix():
# Test precalculation of modes.
rng = ensure_rng(5)
- lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
- lead[chain(0)] = np.zeros((2, 2))
- lead[chain(0), chain(1)] = rng.randn(2, 2)
+ lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)),
+ vectorize=vectorize)
+ lead[chain2(0)] = np.zeros((2, 2))
+ lead[chain2(0), chain2(1)] = rng.randn(2, 2)
syst.attach_lead(lead)
syst2 = syst.finalized()
smatrix = kwant.smatrix(syst2, .1).data
@@ -74,19 +75,26 @@ def test_hamiltonian_submatrix():
raises(ValueError, kwant.solvers.default.greens_function, syst3, 0.2)
# Test for shape errors.
- syst[chain(0), chain(2)] = np.array([[1, 2]])
+ syst[chain2(0), chain(2)] = np.array([[1, 2]])
syst2 = syst.finalized()
raises(ValueError, syst2.hamiltonian_submatrix)
raises(ValueError, syst2.hamiltonian_submatrix, sparse=True)
- syst[chain(0), chain(2)] = 1
+ syst[chain2(0), chain(2)] = 1
syst2 = syst.finalized()
raises(ValueError, syst2.hamiltonian_submatrix)
raises(ValueError, syst2.hamiltonian_submatrix, sparse=True)
+ if vectorize: # non-vectorized systems don't check this at finalization
+ # Add another hopping of the same type but with a different
+ # (and still incompatible) shape.
+ syst[chain2(0), chain(1)] = np.array([[1, 2]])
+ raises(ValueError, syst.finalized)
-def test_pickling():
- syst = kwant.Builder()
- lead = kwant.Builder(symmetry=kwant.TranslationalSymmetry([1.]))
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_pickling(vectorize):
+ syst = kwant.Builder(vectorize=vectorize)
+ lead = kwant.Builder(symmetry=kwant.TranslationalSymmetry([1.]),
+ vectorize=vectorize)
lat = kwant.lattice.chain(norbs=1)
syst[lat(0)] = syst[lat(1)] = 0
syst[lat(0), lat(1)] = 1
diff --git a/kwant/wraparound.py b/kwant/wraparound.py
index 08efec8c668567acb374b12bfbed7777f0571bb5..e5688bab42919e664124fd73d8da0bdf7bb52ae6 100644
--- a/kwant/wraparound.py
+++ b/kwant/wraparound.py
@@ -183,6 +183,9 @@ def wraparound(builder, keep=None, *, coordinate_names='xyz'):
f.__signature__ = inspect.Signature(params.values())
return f
+ if builder.vectorize:
+ raise TypeError("'wraparound' does not work with vectorized Builders.")
+
try:
momenta = ['k_{}'.format(coordinate_names[i])
for i in range(len(builder.symmetry.periods))]