Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • kwant/kwant
  • jbweston/kwant
  • anton-akhmerov/kwant
  • cwg/kwant
  • Mathieu/kwant
  • slavoutich/kwant
  • pacome/kwant
  • behrmann/kwant
  • michaelwimmer/kwant
  • albeercik/kwant
  • eunjongkim/kwant
  • basnijholt/kwant
  • r-j-skolasinski/kwant
  • sahmed95/kwant
  • pablopiskunow/kwant
  • mare/kwant
  • dvarjas/kwant
  • Paul/kwant
  • bbuijtendorp/kwant
  • tkloss/kwant
  • torosdahl/kwant
  • kel85uk/kwant
  • kpoyhonen/kwant
  • Fromeworld/kwant
  • quaeritis/kwant
  • marwahaha/kwant
  • fernandodfufrpe/kwant
  • oly/kwant
  • jiamingh/kwant
  • mehdi2369/kwant
  • ValFadeev/kwant
  • Kostas/kwant
  • chelseabaptiste03/kwant
33 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 2427 additions and 580 deletions
kwant/builder.py
View file @ 44e42c27
# 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
......@@ -13,189 +13,29 @@ import collections
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
from .physics import DiscreteSymmetry, magnetic_gauge
from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
interleave, deprecate_args)
interleave, deprecate_args, memoize)
__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead']
__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 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.
......@@ -301,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.
......@@ -323,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):
......@@ -406,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
......@@ -564,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.
......@@ -572,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.
......@@ -601,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):
......@@ -632,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.
"""
......@@ -663,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.
"""
......@@ -791,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
......@@ -821,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
-----
......@@ -899,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:
......@@ -909,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 = {}
......@@ -992,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
......@@ -1446,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)``.
......@@ -1456,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
------
......@@ -1608,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.
......@@ -1618,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
------
......@@ -1671,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
......@@ -1691,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,
......@@ -1787,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.')
......@@ -1807,6 +1719,176 @@ class Builder:
inter_cell_hopping = cell_hamiltonian = precalculated = \
_require_system
def _add_peierls_phase(syst, peierls_parameter):
@memoize
def wrap_hopping(hop):
def f(*args):
a, b, *args, phi = args
return hop(a, b, *args) * phi(a, b)
params = ('_site0', '_site1')
if callable(hop):
params += get_parameters(hop)[2:] # cut off both site parameters
params += (peierls_parameter,)
f.__signature__ = inspect.Signature(
inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
for name in params)
return f
const_hoppings = {}
def const_hopping(a, b, phi):
return const_hoppings[a, b] * phi(a, b)
# fix the signature
params = ('_site0', '_site1', peierls_parameter)
const_hopping.__signature__ = inspect.Signature(
inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
for name in params)
ret = Builder(
syst.symmetry,
# No time reversal, as magnetic fields break it
conservation_law=syst.conservation_law,
particle_hole=syst.particle_hole,
chiral=syst.chiral,
)
for i, lead in enumerate(syst.leads):
ret.leads.append(
BuilderLead(
_add_peierls_phase(
lead.builder,
peierls_parameter + '_lead{}'.format(i),
),
lead.interface,
lead.padding,
)
)
for site, val in syst.site_value_pairs():
ret[site] = val
for hop, val in syst.hopping_value_pairs():
if callable(val):
ret[hop] = wrap_hopping(val)
else:
const_hoppings[hop] = val
ret[hop] = const_hopping
return ret
def add_peierls_phase(syst, peierls_parameter='phi', fix_gauge=True):
"""Add a Peierls phase parameter to a Builder.
Parameters
----------
syst : kwant.Builder
The system to which to add a Peierls phase.
peierls_parameter : str
The name of the Peierls phase parameter to value functions.
fix_gauge : bool (default: True)
If True, fix the gauge of the system and return it also
Returns
-------
syst : kwant.builder.Finitesystem or kwant.builder.InfiniteSystem
A system where all hoppings are given by value functions
that have an additional parameter 'phi' that should be
the Peierls phase returned by `kwant.physics.magnetic_gauge`.
Any leads have similar parameters named 'phi_lead0',
'phi_lead1' etc.
gauge : callable
Only returned if 'fix_gauge' is True. Called with magnetic
field(s), and returns a parameter dictionary that can be
passed to the system as 'params' (see the example below).
Examples
--------
>>> import numpy as np
>>> import kwant
>>>
>>> syst = make_system()
>>> lead = make_lead()
>>> syst.attach_lead(lead)
>>> syst.attach_lead(lead.reversed())
>>>
>>> fsyst, phase = kwant.builder.add_peierls_phase(syst)
>>>
>>> def B_syst(pos):
>>> return np.exp(-np.sum(pos * pos))
>>>
>>> kwant.smatrix(fsyst, parameters=dict(t=-1, **phase(B_syst, 0, 0)))
"""
def wrap_gauge(gauge):
phase_names = [peierls_parameter]
phase_names += [peierls_parameter + '_lead{}'.format(i)
for i in range(len(syst.leads))]
doc = textwrap.dedent("""\
Return the Peierls phase for a magnetic field configuration.
Parameters
----------
syst_field : scalar, vector or callable
The magnetic field to apply to the scattering region.
If callable, takes a position and returns the
magnetic field at that position. Can be a scalar if
the system is 1D or 2D, otherwise must be a vector.
Magnetic field is expressed in units :math:`φ₀ / l²`,
where :math:`φ₀` is the magnetic flux quantum and
:math:`l` is the unit of length.
*lead_fields : scalar, vector or callable
The magnetic fields to apply to each of the leads, in
the same format as 'syst_field'. In addition, if a callable
is provided, then the magnetic field must have the symmetry
of the associated lead.
tol : float, default: 1E-8
The tolerance to which to calculate the flux through each
hopping loop in the system.
average : bool, default: False
If True, estimate the magnetic flux through each hopping loop
in the system by evaluating the magnetic field at a single
position inside the loop and multiplying it by the area of the
loop. If False, then ``scipy.integrate.quad`` is used to integrate
the magnetic field. This parameter is only used when 'syst_field'
or 'lead_fields' are callable.
Returns
-------
phases : dict of callables
To be passed directly as the parameters of a system wrapped with
'kwant.builder.add_peierls_phase'.
""")
@wraps(gauge)
def f(*args, **kwargs):
phases = gauge(*args, **kwargs)
# When there are no leads 'gauge' returns a single callable
if not syst.leads:
phases = (phases,)
return dict(zip(phase_names, phases))
f.__doc__ = doc
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:
return ret, wrap_gauge(magnetic_gauge(ret))
else:
return ret
################ Finalized systems
......@@ -1946,6 +2028,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 = system._normalize_matrix_blocks(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.
......@@ -1978,6 +2189,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.
......@@ -1986,11 +2255,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``.
"""
......@@ -2004,57 +2273,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
......@@ -2063,7 +2285,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
......@@ -2083,7 +2305,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
......@@ -2096,8 +2318,480 @@ 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 _sort_term(term, value):
term = np.asarray(term)
if not callable(value):
value = system._normalize_matrix_blocks(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):
......@@ -2106,10 +2800,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``.
......@@ -2133,69 +2827,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
......@@ -2203,41 +2840,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:
......@@ -2264,7 +2880,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
......@@ -2275,6 +2891,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
......@@ -2283,5 +3018,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))
kwant/continuum/__init__.py
View file @ 44e42c27
......@@ -10,11 +10,12 @@ try:
from .discretizer import discretize, discretize_symbolic, build_discretized
from ._common import sympify, lambdify
from ._common import momentum_operators, position_operators
from .landau_levels import to_landau_basis, discretize_landau, LandauLattice
except ImportError as error:
msg = ("'kwant.continuum' is not available because one or more of its "
"dependencies is not installed.")
raise ImportError(msg) from error
__all__ = ['discretize', 'discretize_symbolic', 'build_discretized',
'to_landau_basis', 'discretize_landau', 'LandauLattice',
'sympify', 'lambdify', 'momentum_operators', 'position_operators']
kwant/continuum/discretizer.py
View file @ 44e42c27
# Copyright 2011-2017 Kwant authors.
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
......@@ -20,7 +20,7 @@ from sympy.printing.lambdarepr import LambdaPrinter
from sympy.printing.precedence import precedence
from sympy.core.function import AppliedUndef
from .. import builder, lattice
from .. import builder, lattice, system
from .. import KwantDeprecationWarning
from .._common import reraise_warnings
from ._common import (sympify, gcd, position_operators, momentum_operators,
......@@ -63,7 +63,7 @@ class _DiscretizedBuilder(builder.Builder):
for key, val in itertools.chain(self.site_value_pairs(),
self.hopping_value_pairs()):
if isinstance(key, builder.Site):
if isinstance(key, system.Site):
result.append("# Onsite element:\n")
else:
a, b = key
......
kwant/continuum/landau_levels.py 0 → 100644
View file @ 44e42c27
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
# http://kwant-project.org/license. A list of Kwant authors can be found in
# the file AUTHORS.rst at the top-level directory of this distribution and at
# http://kwant-project.org/authors.
from keyword import iskeyword
import functools
import operator
import collections
from math import sqrt
import tinyarray as ta
import numpy as np
import sympy
import kwant.lattice
import kwant.builder
import kwant.continuum
import kwant.continuum._common
import kwant.continuum.discretizer
from kwant.continuum import momentum_operators, position_operators
__all__ = ["to_landau_basis", "discretize_landau"]
coordinate_vectors = dict(zip("xyz", np.eye(3)))
ladder_lower, ladder_raise = sympy.symbols(r"a a^\dagger", commutative=False)
def to_landau_basis(hamiltonian, momenta=None):
r"""Replace two momenta by Landau level ladder operators.
Replaces:
k_0 -> sqrt(B/2) * (a + a^\dagger)
k_1 -> 1j * sqrt(B/2) * (a - a^\dagger)
Parameters
----------
hamiltonian : str or sympy expression
The Hamiltonian to which to apply the Landau level transformation.
momenta : sequence of str (optional)
The momenta to replace with Landau level ladder operators. If not
provided, 'k_x' and 'k_y' are used
Returns
-------
hamiltonian : sympy expression
momenta : sequence of sympy atoms
The momentum operators that have been replaced by ladder operators.
normal_coordinate : sympy atom
The remaining position coordinate. May or may not be present
in 'hamiltonian'.
"""
hamiltonian = kwant.continuum.sympify(hamiltonian)
momenta = _normalize_momenta(momenta)
normal_coordinate = _find_normal_coordinate(hamiltonian, momenta)
# Substitute ladder operators for Landau level momenta
B = sympy.symbols("B")
hamiltonian = hamiltonian.subs(
{
momenta[0]: sympy.sqrt(abs(B) / 2) * (ladder_raise + ladder_lower),
momenta[1]: sympy.I
* sympy.sqrt(abs(B) / 2)
* (ladder_lower - ladder_raise),
}
)
return hamiltonian, momenta, normal_coordinate
def discretize_landau(hamiltonian, N, momenta=None, grid_spacing=1):
"""Discretize a Hamiltonian in a basis of Landau levels.
Parameters
----------
hamiltonian : str or sympy expression
N : positive integer
The number of Landau levels in the basis.
momenta : sequence of str (optional)
The momenta defining the plane perpendicular to the magnetic field.
If not provided, "k_x" and "k_y" are used.
grid_spacing : float, default: 1
The grid spacing to use when discretizing the normal direction
(parallel to the magnetic field).
Returns
-------
builder : `~kwant.builder.Builder`
'hamiltonian' discretized in a basis of Landau levels in the plane
defined by 'momenta'. If a third coordinate is present in 'hamiltonian',
then this also has a translational symmetry in that coordinate direction.
The builder has a parameter 'B' in addition to any other parameters
present in the provided 'hamiltonian'.
Notes
-----
The units of magnetic field are :math:`ϕ₀ / 2 π a²` with :math:`ϕ₀ = h/e`
the magnetic flux quantum and :math:`a` the unit length.
"""
if N <= 0:
raise ValueError("N must be positive")
hamiltonian, momenta, normal_coordinate = to_landau_basis(hamiltonian, momenta)
# Discretize in normal direction and split terms for onsites/hoppings into
# monomials in ladder operators.
tb_hamiltonian, _ = kwant.continuum.discretize_symbolic(
hamiltonian, coords=[normal_coordinate.name]
)
tb_hamiltonian = {
key: kwant.continuum._common.monomials(value, gens=(ladder_lower, ladder_raise))
for key, value in tb_hamiltonian.items()
}
# Replace ladder operator monomials by tuple of integers:
# e.g. a^\dagger a^2 -> (+1, -2).
tb_hamiltonian = {
outer_key: {
_ladder_term(inner_key): inner_value
for inner_key, inner_value in outer_value.items()
}
for outer_key, outer_value in tb_hamiltonian.items()
}
# Construct map from LandauLattice HoppingKinds to a sequence of pairs
# that encode the ladder operators and multiplying expression.
tb_hoppings = collections.defaultdict(list)
for outer_key, outer_value in tb_hamiltonian.items():
for inner_key, inner_value in outer_value.items():
tb_hoppings[(*outer_key, sum(inner_key))] += [(inner_key, inner_value)]
# Extract the number of orbitals on each site/hopping
random_element = next(iter(tb_hoppings.values()))[0][1]
norbs = 1 if isinstance(random_element, sympy.Expr) else random_element.shape[0]
tb_onsite = tb_hoppings.pop((0, 0), None)
# Construct Builder
if _has_coordinate(normal_coordinate, hamiltonian):
sym = kwant.lattice.TranslationalSymmetry([grid_spacing, 0])
else:
sym = kwant.builder.NoSymmetry()
lat = LandauLattice(grid_spacing, norbs=norbs)
syst = kwant.Builder(sym)
# Add onsites
landau_sites = (lat(0, j) for j in range(N))
if tb_onsite is None:
syst[landau_sites] = ta.zeros((norbs, norbs))
else:
syst[landau_sites] = _builder_value(
tb_onsite, normal_coordinate.name, grid_spacing, is_onsite=True
)
# Add zero hoppings between adjacent Landau levels.
# Necessary to be able to use the Landau level builder
# to populate another builder using builder.fill().
syst[kwant.builder.HoppingKind((0, 1), lat)] = ta.zeros((norbs, norbs))
# Add the hoppings from the Hamiltonian
for hopping, parts in tb_hoppings.items():
syst[kwant.builder.HoppingKind(hopping, lat)] = _builder_value(
parts, normal_coordinate.name, grid_spacing, is_onsite=False
)
return syst
# This has to subclass lattice so that it will work with TranslationalSymmetry.
class LandauLattice(kwant.lattice.Monatomic):
"""
A `~kwant.lattice.Monatomic` lattice with a Landau level index per site.
Site tags (see `~kwant.system.SiteFamily`) are pairs of integers, where
the first integer describes the real space position and the second the
Landau level index.
The real space Bravais lattice is one dimensional, oriented parallel
to the magnetic field. The Landau level index represents the harmonic
oscillator basis states used for the Landau quantization in the plane
normal to the field.
Parameters
----------
grid_spacing : float
Real space lattice spacing (parallel to the magnetic field).
offset : float (optional)
Displacement of the lattice origin from the real space
coordinates origin.
"""
def __init__(self, grid_spacing, offset=None, name="", norbs=None):
if offset is not None:
offset = [offset, 0]
# The second vector and second coordinate do not matter (they are
# not used in pos())
super().__init__([[grid_spacing, 0], [0, 1]], offset, name, norbs)
def pos(self, tag):
return ta.array((self.prim_vecs[0, 0] * tag[0] + self.offset[0],))
def landau_index(self, tag):
return tag[-1]
def _builder_value(terms, normal_coordinate, grid_spacing, is_onsite):
"""Construct an onsite/hopping value function from a list of terms
Parameters
----------
terms : list
Each element is a pair (ladder_term, sympy expression/matrix).
ladder_term is a tuple of integers that encodes a string of
Landau raising/lowering operators and the sympy expression
is the rest
normal_coordinate : str
grid_spacing : float
The grid spacing in the normal direction
is_onsite : bool
True if we are constructing an onsite value function
"""
ladder_term_symbols = [sympy.Symbol(_ladder_term_name(lt)) for lt, _ in terms]
# Construct a single expression from the terms, multiplying each part
# by the placeholder that represents the prefactor from the ladder operator term.
(ladder, (_, part)), *rest = zip(ladder_term_symbols, terms)
expr = ladder * part
for ladder, (_, part) in rest:
expr += ladder * part
expr = expr.subs(
{kwant.continuum.discretizer._displacements[normal_coordinate]: grid_spacing}
)
# Construct the return string and temporary variable names
# for function calls.
return_string, map_func_calls, const_symbols, _cache = kwant.continuum.discretizer._return_string(
expr, coords=[normal_coordinate]
)
# Remove the ladder term placeholders, as these are not parameters
const_symbols = set(const_symbols)
for ladder_term in ladder_term_symbols:
const_symbols.discard(ladder_term)
# Construct the argument part of signature. Arguments
# consist of any constants and functions in the return string.
arg_names = set.union(
{s.name for s in const_symbols}, {str(k.func) for k in map_func_calls}
)
arg_names.discard("Abs") # Abs function is not a parameter
for arg_name in arg_names:
if not (arg_name.isidentifier() and not iskeyword(arg_name)):
raise ValueError(
"Invalid name in used symbols: {}\n"
"Names of symbols used in Hamiltonian "
"must be valid Python identifiers and "
"may not be keywords".format(arg_name)
)
arg_names = ", ".join(sorted(arg_names))
# Construct site part of the function signature
site_string = "from_site" if is_onsite else "to_site, from_site"
# Construct function signature
if arg_names:
function_header = "def _({}, {}):".format(site_string, arg_names)
else:
function_header = "def _({}):".format(site_string)
# Construct function body
function_body = []
if "B" not in arg_names:
# B is not a parameter for terms with no ladder operators but we still
# need something to pass to _evaluate_ladder_term
function_body.append("B = +1")
function_body.extend(
[
"{}, = from_site.pos".format(normal_coordinate),
"_ll_index = from_site.family.landau_index(from_site.tag)",
]
)
# To get the correct hopping if B < 0, we need to set the Hermitian
# conjugate of the ladder operator matrix element, which swaps the
# from_site and to_site Landau level indices.
if not is_onsite:
function_body.extend(
["if B < 0:",
" _ll_index = to_site.family.landau_index(to_site.tag)"
]
)
function_body.extend(
"{} = _evaluate_ladder_term({}, _ll_index, B)".format(symb.name, lt)
for symb, (lt, _) in zip(ladder_term_symbols, terms)
)
function_body.extend(
"{} = {}".format(v, kwant.continuum.discretizer._print_sympy(k))
for k, v in map_func_calls.items()
)
function_body.append(return_string)
func_code = "\n ".join([function_header] + function_body)
# Add "I" to namespace just in case sympy again would miss to replace it
# with Python's 1j as it was the case with SymPy 1.2 when "I" was argument
# of some function.
namespace = {
"pi": np.pi,
"I": 1j,
"_evaluate_ladder_term": _evaluate_ladder_term,
"Abs": abs,
}
namespace.update(_cache)
# Construct full source, including cached arrays
source = []
for k, v in _cache.items():
source.append("{} = (\n{})\n".format(k, repr(np.array(v))))
source.append(func_code)
exec(func_code, namespace)
f = namespace["_"]
f._source = "".join(source)
return f
def _ladder_term(operator_string):
r"""Return a tuple of integers representing a string of ladder operators
Parameters
----------
operator_string : Sympy expression
Monomial in ladder operators, e.g. a^\dagger a^2 a^\dagger.
Returns
-------
ladder_term : tuple of int
e.g. a^\dagger a^2 -> (+1, -2)
"""
ret = []
for factor in operator_string.as_ordered_factors():
ladder_op, exponent = factor.as_base_exp()
if ladder_op == ladder_lower:
sign = -1
elif ladder_op == ladder_raise:
sign = +1
else:
sign = 0
ret.append(sign * int(exponent))
return tuple(ret)
def _ladder_term_name(ladder_term):
"""
Parameters
----------
ladder_term : tuple of int
Returns
-------
ladder_term_name : str
"""
def ns(i):
if i >= 0:
return str(i)
else:
return "_" + str(-i)
return "_ladder_{}".format("_".join(map(ns, ladder_term)))
def _evaluate_ladder_term(ladder_term, n, B):
r"""Evaluates the prefactor for a ladder operator on a landau level.
Example: a^\dagger a^2 -> (n - 1) * sqrt(n)
Parameters
----------
ladder_term : tuple of int
Represents a string of ladder operators. Positive
integers represent powers of the raising operator,
negative integers powers of the lowering operator.
n : non-negative int
Landau level index on which to act with ladder_term.
B : float
Magnetic field with sign
Returns
-------
ladder_term_prefactor : float
"""
assert n >= 0
# For negative B we swap a and a^\dagger.
if B < 0:
ladder_term = tuple(-i for i in ladder_term)
ret = 1
for m in reversed(ladder_term):
if m > 0:
factors = range(n + 1, n + m + 1)
elif m < 0:
factors = range(n + m + 1, n + 1)
if n == 0:
return 0 # a|0> = 0
else:
factors = (1,)
ret *= sqrt(functools.reduce(operator.mul, factors))
n += m
return ret
def _normalize_momenta(momenta=None):
"""Return Landau level momenta as Sympy atoms
Parameters
----------
momenta : None or pair of int or pair of str
The momenta to choose. If None then 'k_x' and 'k_y'
are chosen. If integers, then these are the indices
of the momenta: 0 → k_x, 1 → k_y, 2 → k_z. If strings,
then these name the momenta.
Returns
-------
The specified momenta as sympy atoms.
"""
# Specify which momenta to substitute for the Landau level basis.
if momenta is None:
# Use k_x and k_y by default
momenta = momentum_operators[:2]
else:
if len(momenta) != 2:
raise ValueError("Two momenta must be specified.")
k_names = [k.name for k in momentum_operators]
if all([type(i) is int for i in momenta]) and all(
[i >= 0 and i < 3 for i in momenta]
):
momenta = [momentum_operators[i] for i in momenta]
elif all([isinstance(momentum, str) for momentum in momenta]) and all(
[momentum in k_names for momentum in momenta]
):
momenta = [
momentum_operators[k_names.index(momentum)] for momentum in momenta
]
else:
raise ValueError("Momenta must all be integers or strings.")
return tuple(momenta)
def _find_normal_coordinate(hamiltonian, momenta):
discrete_momentum = next(
momentum for momentum in momentum_operators if momentum not in momenta
)
normal_coordinate = position_operators[momentum_operators.index(discrete_momentum)]
return normal_coordinate
def _has_coordinate(coord, expr):
momentum = momentum_operators[position_operators.index(coord)]
atoms = set(expr.atoms())
return coord in atoms or momentum in atoms
kwant/continuum/tests/test_discretizer.py
View file @ 44e42c27
......@@ -548,14 +548,14 @@ def test_grid(ham, grid_spacing, grid):
@pytest.mark.parametrize('ham, grid_offset, offset, norbs', [
('k_x', None, 0, None),
('k_x', None, 0, 1),
('k_x', None, 0, 1),
('k_x * eye(2)', None, 0, 2),
('k_x', (0,), 0, None),
('k_x', (1,), 1, None),
('k_x + k_y', None, (0, 0), None),
('k_x + k_y', (0, 0), (0, 0), None),
('k_x + k_y', (1, 2), (1, 2), None),
('k_x', (0,), 0, 1),
('k_x', (1,), 1, 1),
('k_x + k_y', None, (0, 0), 1),
('k_x + k_y', (0, 0), (0, 0), 1),
('k_x + k_y', (1, 2), (1, 2), 1),
])
def test_grid_input(ham, grid_offset, offset, norbs):
# build appriopriate grid
......@@ -579,7 +579,7 @@ def test_grid_input(ham, grid_offset, offset, norbs):
def test_grid_offset_passed_to_functions():
V = lambda x: x
grid = Monatomic([[1, ]], offset=[0.5, ])
grid = Monatomic([[1, ]], offset=[0.5, ], norbs=1)
tb = discretize('V(x)', 'x', grid=grid)
onsite = tb[tb.lattice(0)]
bools = [np.allclose(onsite(tb.lattice(i), V), V(tb.lattice(i).pos))
......@@ -588,11 +588,11 @@ def test_grid_offset_passed_to_functions():
@pytest.mark.parametrize("ham, coords, grid", [
("k_x", None, Monatomic([[1, 0]])),
("k_x", 'xy', Monatomic([[1, 0]])),
("k_x", None, Monatomic([[1, 0]], norbs=1)),
("k_x", 'xy', Monatomic([[1, 0]], norbs=1)),
("k_x", None, Monatomic([[1, ]], norbs=2)),
("k_x * eye(2)", None, Monatomic([[1, ]], norbs=1)),
("k_x+k_y", None, Monatomic([[1, 0], [1, 1]])),
("k_x+k_y", None, Monatomic([[1, 0], [1, 1]], norbs=1)),
])
def test_grid_constraints(ham, coords, grid):
with pytest.raises(ValueError):
......@@ -606,7 +606,7 @@ def test_check_symbol_names(name):
def test_rectangular_grid():
lat = Monatomic([[1, 0], [0, 2]])
lat = Monatomic([[1, 0], [0, 2]], norbs=1)
tb = discretize("V(x, y)", 'xy', grid=lat)
assert np.allclose(tb[lat(0, 0)](lat(1, 0), lambda x, y: x), 1)
......
kwant/continuum/tests/test_landau_levels.py 0 → 100644
View file @ 44e42c27
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
# http://kwant-project.org/license. A list of Kwant authors can be found in
# the file AUTHORS.rst at the top-level directory of this distribution and at
# http://kwant-project.org/authors.
from math import sqrt
import numpy as np
import sympy
import pytest
import itertools
import kwant.builder
import kwant.lattice
from .._common import position_operators, momentum_operators, sympify
from ..landau_levels import (
to_landau_basis,
discretize_landau,
_ladder_term,
_evaluate_ladder_term,
ladder_lower,
ladder_raise,
LandauLattice,
)
x, y, z = position_operators
k_x, k_y, k_z = momentum_operators
B = sympy.symbols("B")
V = sympy.symbols("V", cls=sympy.Function)
a, a_dag = ladder_lower, ladder_raise
def test_ladder_term():
assert _ladder_term(a ** 2 * a_dag) == (-2, 1)
assert _ladder_term(a_dag ** 5 * a ** 3 * a_dag) == (5, -3, 1)
# non-ladder operators give a shift of 0
assert _ladder_term(B * a ** 2) == (0, -2)
def test_evaluate_ladder_term():
assert np.isclose(_evaluate_ladder_term((+1, -1), 1, +1), 1)
assert np.isclose(_evaluate_ladder_term((+1, -1), 2, +1), 2)
assert np.isclose(_evaluate_ladder_term((-2, +3, -2), 5, +1), 4 * 5 * 6 * sqrt(5))
# annihilating |0> is always 0
assert _evaluate_ladder_term((-1,), 0, +1) == 0
assert _evaluate_ladder_term((-2,), 1, +1) == 0
assert _evaluate_ladder_term((-3,), 1, +1) == 0
assert _evaluate_ladder_term((+3, -2), 1, +1) == 0
assert _evaluate_ladder_term((-3, -2, +3), 1, +1) == 0
# negative B swaps creation and annihilation operators
assert _evaluate_ladder_term((+1, -1), 2, +1) == _evaluate_ladder_term(
(-1, +1), 2, -1
)
assert _evaluate_ladder_term((-2, +3, -2), 5, +1) == _evaluate_ladder_term(
(+2, -3, +2), 5, -1
)
def test_to_landau_basis():
# test basic usage
ham, momenta, normal_coord = to_landau_basis("k_x**2 + k_y**2")
assert sympy.expand(ham) == abs(B) * a * a_dag + abs(B) * a_dag * a
assert momenta == (k_x, k_y)
assert normal_coord == z
# test that hamiltonian can be specified as a sympy expression
ham, momenta, normal_coord = to_landau_basis(sympify("k_x**2 + k_y**2"))
assert sympy.expand(ham) == abs(B) * a * a_dag + abs(B) * a_dag * a
assert momenta == (k_x, k_y)
assert normal_coord == z
# test that
ham, momenta, normal_coord = to_landau_basis("k_x**2 + k_y**2 + k_z**2 + V(z)")
assert sympy.expand(ham) == (
abs(B) * a * a_dag + abs(B) * a_dag * a + k_z ** 2 + V(z)
)
assert momenta == (k_x, k_y)
assert normal_coord == z
# test for momenta explicitly specified
ham, momenta, normal_coord = to_landau_basis(
"k_x**2 + k_y**2 + k_z**2 + k_x*k_y", momenta=("k_z", "k_y")
)
assert sympy.expand(ham) == (
abs(B) * a * a_dag
+ abs(B) * a_dag * a
+ k_x ** 2
+ sympy.I * sympy.sqrt(abs(B) / 2) * k_x * a
- sympy.I * sympy.sqrt(abs(B) / 2) * k_x * a_dag
)
assert momenta == (k_z, k_y)
assert normal_coord == x
def test_discretize_landau():
n_levels = 10
magnetic_field = 1 / 3 # a suitably arbitrary value
# Ensure that we can handle products of ladder operators by truncating
# several levels higher than the number of levels we actually want.
a = np.diag(np.sqrt(np.arange(1, n_levels + 5)), k=1)
a_dag = a.conjugate().transpose()
k_x = sqrt(magnetic_field / 2) * (a + a_dag)
k_y = 1j * sqrt(magnetic_field / 2) * (a - a_dag)
sigma_0 = np.eye(2)
sigma_x = np.array([[0, 1], [1, 0]])
sigma_y = np.array([[0, -1j], [1j, 0]])
# test that errors are raised on invalid input
with pytest.raises(ValueError):
discretize_landau("k_x**2 + k_y**2", N=0)
with pytest.raises(ValueError):
discretize_landau("k_x**2 + k_y**2", N=-1)
# test a basic Hamiltonian with no normal coordinate dependence
syst = discretize_landau("k_x**2 + k_y**2", N=n_levels)
lat = LandauLattice(1, norbs=1)
assert isinstance(syst.symmetry, kwant.builder.NoSymmetry)
syst = syst.finalized()
assert set(syst.sites) == {lat(0, j) for j in range(n_levels)}
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=0)), np.zeros((n_levels, n_levels))
)
should_be = k_x @ k_x + k_y @ k_y
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
should_be[:n_levels, :n_levels],
)
# test negative magnetic field
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=-magnetic_field)),
should_be[:n_levels, :n_levels],
)
# test hamiltonian with no onsite elements
syst = discretize_landau("k_x", N=n_levels)
syst = syst.finalized()
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
k_x[:n_levels, :n_levels],
)
# test a basic Hamiltonian with normal coordinate dependence
grid = 1 / 7 # a suitably arbitrary value
syst = discretize_landau(
"k_x**2 + k_y**2 + k_z**2 + V(z)", N=n_levels, grid_spacing=grid
)
assert isinstance(syst.symmetry, kwant.lattice.TranslationalSymmetry)
syst = syst.finalized()
zero_potential = syst.cell_hamiltonian(params=dict(B=magnetic_field, V=lambda z: 0))
with_potential = syst.cell_hamiltonian(params=dict(B=magnetic_field, V=lambda z: 1))
# extra +2/grid**2 from onsite kinetic energy
assert np.allclose(
zero_potential,
should_be[:n_levels, :n_levels] + (2 / grid ** 2) * np.eye(n_levels),
)
# V(z) just adds the same energy to each onsite
assert np.allclose(with_potential - zero_potential, np.eye(n_levels))
# hopping matrix does not exchange landau levels
assert np.allclose(
syst.inter_cell_hopping(params=dict(B=magnetic_field, V=lambda z: 0)),
-np.eye(n_levels) / grid ** 2,
)
# test a Hamiltonian with mixing between Landau levels
# and spatial degrees of freedom.
syst = discretize_landau("k_x**2 + k_y**2 + k_x*k_z", N=n_levels)
syst = syst.finalized()
assert np.allclose(
syst.inter_cell_hopping(params=dict(B=magnetic_field)),
-1j * k_x[:n_levels, :n_levels] / 2,
)
# test a Hamiltonian with extra degrees of freedom
syst = discretize_landau("sigma_0 * k_x**2 + sigma_x * k_y**2", N=n_levels)
syst = syst.finalized()
assert syst.sites[0].family.norbs == 2
should_be = np.kron(k_x @ k_x, sigma_0) + np.kron(k_y @ k_y, sigma_x)
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
should_be[: 2 * n_levels, : 2 * n_levels],
)
# test a linear Hamiltonian
syst = discretize_landau("sigma_y * k_x - sigma_x * k_y", N=n_levels)
syst = syst.finalized()
should_be = np.kron(k_x, sigma_y) - np.kron(k_y, sigma_x)
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
should_be[: 2 * n_levels, : 2 * n_levels],
)
assert np.allclose(
syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
syst.hamiltonian_submatrix(params=dict(B=-magnetic_field)),
)
def test_analytical_spectrum():
hamiltonian = """(k_x**2 + k_y**2) * sigma_0 +
alpha * (k_x * sigma_y - k_y * sigma_x) +
g * B * sigma_z"""
def exact_Es(n, B, alpha, g):
# See e.g. R. Winkler (2003), section 8.4.1
sign_B = np.sign(B)
B = np.abs(B)
Ep = 2*B*(n+1) - 0.5*np.sqrt((2*B - sign_B*2*g*B)**2 + 8*B*alpha**2*(n+1))
Em = 2*B*n + 0.5*np.sqrt((2*B - sign_B*2*g*B)**2 + 8*B*alpha**2*n)
return Ep, Em
N = 20
syst = discretize_landau(hamiltonian, N)
syst = syst.finalized()
params = dict(alpha=0.07, g=0.04)
for _ in range(5):
B = 0.01 + np.random.rand()*3
params['B'] = B
exact = [exact_Es(n, B, params['alpha'], params['g']) for n in range(N)]
# We only check the first N levels - the SOI couples adjacent levels,
# so the higher ones are not necessarily determined accurately in the
# discretization
exact = np.sort([energy for energies in exact for energy in energies])[:N]
ll_spect = np.linalg.eigvalsh(syst.hamiltonian_submatrix(params=params))[:len(exact)]
assert np.allclose(ll_spect, exact)
def test_fill():
def shape(site, lower, upper):
(z, )= site.pos
return lower <= z < upper
hamiltonian = "k_x**2 + k_y**2 + k_z**2"
N = 6
template = discretize_landau(hamiltonian, N)
syst = kwant.Builder()
width = 4
syst.fill(template, lambda site: shape(site, 0, width), (0, ));
correct_tags = [(coordinate, ll_index) for coordinate, ll_index
in itertools.product(range(width), range(N))]
syst_tags = [site.tag for site in syst.sites()]
assert len(syst_tags) == len(correct_tags)
assert all(tag in correct_tags for tag in syst_tags)
kwant/graph/__init__.py
View file @ 44e42c27
......@@ -9,8 +9,7 @@
"""Functionality for graphs"""
# Merge the public interface of all submodules.
__all__ = []
for module in ['core', 'defs']:
exec('from . import {0}'.format(module))
exec('from .{0} import *'.format(module))
exec('__all__.extend({0}.__all__)'.format(module))
from .core import *
from .defs import *
__all__ = [core.__all__ + defs.__all__]
kwant/kpm.py
View file @ 44e42c27
# -*- coding: utf-8 -*-
# Copyright 2011-2016 Kwant authors.
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
......@@ -731,23 +731,25 @@ class Correlator:
g_kernel[0] /= 2
mu_kernel = np.outer(g_kernel, g_kernel) * self.moments_matrix
e = (self.energies - self._b) / self._a
e_scaled = (self.energies - self._b) / self._a
# arrays for faster calculation
sqrt_e = np.sqrt(1 - e ** 2)
arccos_e = np.arccos(e)
m_array = np.arange(n_moments)
def _integral_factor(e):
# arrays for faster calculation
sqrt_e = np.sqrt(1 - e ** 2)
arccos_e = np.arccos(e)
exp_n = np.exp(1j * np.outer(arccos_e, np.arange(n_moments)))
t_n = np.real(exp_n)
exp_n = np.exp(1j * arccos_e * m_array)
t_n = np.real(exp_n)
e_plus = (np.outer(e, np.ones(n_moments)) -
1j * np.outer(sqrt_e, np.arange(n_moments)))
e_plus = e_plus * exp_n
e_plus = (e - 1j * sqrt_e * m_array)
e_plus = e_plus * exp_n
big_gamma = e_plus[:, None, :] * t_n[:, :, None]
big_gamma += big_gamma.conj().swapaxes(1, 2)
self._integral_factor = np.tensordot(mu_kernel, big_gamma.T)
big_gamma = e_plus[None, :] * t_n[:, None]
big_gamma += big_gamma.conj().T
return np.tensordot(mu_kernel, big_gamma.T)
self._integral_factor = np.array([_integral_factor(e)
for e in e_scaled]).T
def conductivity(hamiltonian, alpha='x', beta='x', positions=None, **kwargs):
......@@ -926,7 +928,7 @@ def RandomVectors(syst, where=None, rng=None):
----------
syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
If a system is passed, it should contain no leads.
where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
where : sequence of `int` or `~kwant.system.Site`, or callable, optional
Spatial range of the vectors produced. If ``syst`` is a
`~kwant.system.FiniteSystem`, where behaves as in
`~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
......@@ -948,7 +950,7 @@ class LocalVectors:
----------
syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
If a system is passed, it should contain no leads.
where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
where : sequence of `int` or `~kwant.system.Site`, or callable, optional
Spatial range of the vectors produced. If ``syst`` is a
`~kwant.system.FiniteSystem`, where behaves as in
`~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
......
kwant/lattice.py
View file @ 44e42c27
# Copyright 2011-2013 Kwant authors.
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
......@@ -13,7 +13,7 @@ from math import sqrt
from itertools import product
import numpy as np
import tinyarray as ta
from . import builder
from . import builder, system
from .linalg import lll
from ._common import ensure_isinstance
......@@ -70,7 +70,7 @@ class Polyatomic:
A Bravais lattice with an arbitrary number of sites in the basis.
Contains `Monatomic` sublattices. Note that an instance of ``Polyatomic`` is
not itself a `~kwant.builder.SiteFamily`, only its sublattices are.
not itself a `~kwant.system.SiteFamily`, only its sublattices are.
Parameters
----------
......@@ -171,7 +171,7 @@ class Polyatomic:
>>> syst[lat.neighbors()] = 1
"""
def shape_sites(symmetry=None):
Site = builder.Site
Site = system.Site
if symmetry is None:
symmetry = builder.NoSymmetry()
......@@ -306,7 +306,7 @@ class Polyatomic:
"""
# This algorithm is not designed to be fast and can be improved,
# however there is no real need.
Site = builder.Site
Site = system.Site
sls = self.sublattices
shortest_hopping = sls[0].n_closest(
sls[0].pos(([0] * sls[0].lattice_dim)), 2)[-1]
......@@ -396,12 +396,12 @@ def short_array_str(array):
return full[1 : -1]
class Monatomic(builder.SiteFamily, Polyatomic):
class Monatomic(system.SiteFamily, Polyatomic):
"""
A Bravais lattice with a single site in the basis.
Instances of this class provide the `~kwant.builder.SiteFamily` interface.
Site tags (see `~kwant.builder.SiteFamily`) are sequences of integers and
Instances of this class provide the `~kwant.system.SiteFamily` interface.
Site tags (see `~kwant.system.SiteFamily`) are sequences of integers and
describe the lattice coordinates of a site.
``Monatomic`` instances are used as site families on their own or as
......@@ -470,6 +470,13 @@ class Monatomic(builder.SiteFamily, Polyatomic):
def __str__(self):
return self.cached_str
def normalize_tags(self, tags):
tags = np.asarray(tags, int)
tags.flags.writeable = False
if tags.shape[1] != self.lattice_dim:
raise ValueError("Dimensionality mismatch.")
return tags
def normalize_tag(self, tag):
tag = ta.array(tag, int)
if len(tag) != self.lattice_dim:
......@@ -495,6 +502,10 @@ class Monatomic(builder.SiteFamily, Polyatomic):
"""
return ta.array(self.n_closest(pos)[0])
def positions(self, tags):
"""Return the real-space positions of the sites with the given tags."""
return tags @ self._prim_vecs + self.offset
def pos(self, tag):
"""Return the real-space position of the site with a given tag."""
return ta.dot(tag, self._prim_vecs) + self.offset
......@@ -687,26 +698,48 @@ class TranslationalSymmetry(builder.Symmetry):
def which(self, site):
det_x_inv_m_part, det_m = self._get_site_family_data(site.family)[-2:]
result = ta.dot(det_x_inv_m_part, site.tag) // det_m
if isinstance(site, system.Site):
result = ta.dot(det_x_inv_m_part, site.tag) // det_m
elif isinstance(site, system.SiteArray):
result = np.dot(det_x_inv_m_part, site.tags.transpose()) // det_m
else:
raise TypeError("'site' must be a Site or a SiteArray")
return -result if self.is_reversed else result
def act(self, element, a, b=None):
element = ta.array(element)
if element.dtype is not int:
is_site = isinstance(a, system.Site)
# Tinyarray for small arrays (single site) else numpy
array_mod = ta if is_site else np
element = array_mod.array(element)
if not np.issubdtype(element.dtype, np.integer):
raise ValueError("group element must be a tuple of integers")
if (len(element.shape) == 2 and is_site):
raise ValueError("must provide a single group element when "
"acting on single sites.")
if (len(element.shape) == 1 and not is_site):
raise ValueError("must provide a sequence of group elements "
"when acting on site arrays.")
m_part = self._get_site_family_data(a.family)[0]
try:
delta = ta.dot(m_part, element)
delta = array_mod.dot(m_part, element)
except ValueError:
msg = 'Expecting a {0}-tuple group element, but got `{1}` instead.'
raise ValueError(msg.format(self.num_directions, element))
if self.is_reversed:
delta = -delta
if b is None:
return builder.Site(a.family, a.tag + delta, True)
if is_site:
return system.Site(a.family, a.tag + delta, True)
else:
return system.SiteArray(a.family, a.tags + delta.transpose())
elif b.family == a.family:
return (builder.Site(a.family, a.tag + delta, True),
builder.Site(b.family, b.tag + delta, True))
if is_site:
return (system.Site(a.family, a.tag + delta, True),
system.Site(b.family, b.tag + delta, True))
else:
return (system.SiteArray(a.family, a.tags + delta.transpose()),
system.SiteArray(b.family, b.tags + delta.transpose()))
else:
m_part = self._get_site_family_data(b.family)[0]
try:
......@@ -717,8 +750,12 @@ class TranslationalSymmetry(builder.Symmetry):
raise ValueError(msg.format(self.num_directions, element))
if self.is_reversed:
delta2 = -delta2
return (builder.Site(a.family, a.tag + delta, True),
builder.Site(b.family, b.tag + delta2, True))
if is_site:
return (system.Site(a.family, a.tag + delta, True),
system.Site(b.family, b.tag + delta2, True))
else:
return (system.SiteArray(a.family, a.tags + delta.transpose()),
system.SiteArray(b.family, b.tags + delta2.transpose()))
def reversed(self):
"""Return a reversed copy of the symmetry.
......
kwant/linalg/__init__.py
View file @ 44e42c27
......@@ -10,7 +10,10 @@ __all__ = ['lapack']
from . import lapack
# Merge the public interface of the other submodules.
for module in ['decomp_lu', 'decomp_ev', 'decomp_schur']:
exec('from . import {0}'.format(module))
exec('from .{0} import *'.format(module))
exec('__all__.extend({0}.__all__)'.format(module))
from .decomp_lu import *
from .decomp_schur import *
from .decomp_ev import *
__all__.extend([decomp_lu.__all__,
decomp_ev.__all__,
decomp_schur.__all__])
kwant/linalg/tests/test_mumps.py
View file @ 44e42c27
......@@ -65,7 +65,7 @@ def test_schur_complement_with_dense():
def test_error_minus_9(r=10):
"""Test if MUMPSError -9 is properly caught by increasing memory"""
graphene = honeycomb()
graphene = honeycomb(norbs=1)
a, b = graphene.sublattices
def circle(pos):
......
kwant/operator.pxd
View file @ 44e42c27
......@@ -4,15 +4,24 @@ from .graph.defs import gint_dtype
cdef gint _bisect(gint[:] a, gint x)
cdef int _is_herm_conj(complex[:, :] a, complex[:, :] b,
double atol=*, double rtol=*) except -1
cdef int _is_hermitian(
complex[:, :] a, double atol=*, double rtol=*
) except -1
cdef int _is_hermitian_3d(
complex[:, :, :] a, double atol=*, double rtol=*
) except -1
cdef _select(gint[:, :] arr, gint[:] indexes)
cdef int _check_onsite(complex[:, :] M, gint norbs,
int check_hermiticity) except -1
cdef int _check_ham(complex[:, :] H, ham, args, params,
gint a, gint a_norbs, gint b, gint b_norbs,
int check_hermiticity) except -1
cdef int _check_onsites(complex[:, :, :] M, gint norbs,
int check_hermiticity) except -1
cdef int _check_hams(complex[:, :, :] H, gint to_norbs, gint from_norbs,
int check_hermiticity) except -1
cdef void _get_orbs(gint[:, :] site_ranges, gint site,
gint *start_orb, gint *norbs)
......@@ -28,7 +37,7 @@ cdef class BlockSparseMatrix:
cdef class _LocalOperator:
cdef public int check_hermiticity, sum
cdef public object syst, onsite, _onsite_param_names
cdef public object syst, onsite, _onsite_param_names, _terms
cdef public gint[:, :] where, _site_ranges
cdef public BlockSparseMatrix _bound_onsite, _bound_hamiltonian
......
kwant/operator.pyx
View file @ 44e42c27
# 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
......@@ -10,9 +10,10 @@
__all__ = ['Density', 'Current', 'Source']
import cython
from operator import itemgetter
import itertools
import functools as ft
import collections
import warnings
import numpy as np
import tinyarray as ta
......@@ -20,11 +21,19 @@ from scipy.sparse import coo_matrix
from libc cimport math
cdef extern from "complex.h":
double cabs(double complex)
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, is_vectorized, Site, SiteArray, _normalize_matrix_blocks
)
from ._common import (
UserCodeError, KwantDeprecationWarning, get_parameters, deprecate_args
)
################ Generic Utility functions
......@@ -45,32 +54,75 @@ cdef gint _bisect(gint[:] a, gint x):
@cython.boundscheck(False)
@cython.wraparound(False)
cdef int _is_herm_conj(complex[:, :] a, complex[:, :] b,
double atol=1e-300, double rtol=1e-13) except -1:
"Return True if `a` is the Hermitian conjugate of `b`."
assert a.shape[0] == b.shape[1]
assert a.shape[1] == b.shape[0]
cdef int _is_hermitian(
complex[:, :] a, double atol=1e-300, double rtol=1e-13
) except -1:
"Return True if 'a' is Hermitian"
if a.shape[0] != a.shape[1]:
return False
# compute max(a)
cdef double tmp, max_a = 0
cdef gint i, j
cdef gint i, j, k
for i in range(a.shape[0]):
for j in range(a.shape[1]):
tmp = a[i, j].real * a[i, j].real + a[i, j].imag * a[i, j].imag
tmp = cabs(a[i, j])
if tmp > max_a:
max_a = tmp
max_a = math.sqrt(max_a)
cdef double tol = rtol * max_a + atol
cdef complex ctmp
for i in range(a.shape[0]):
for j in range(a.shape[1]):
ctmp = a[i, j] - b[j, i].conjugate()
tmp = ctmp.real * ctmp.real + ctmp.imag * ctmp.imag
for j in range(i, a.shape[1]):
tmp = cabs(a[i, j] - a[j, i].conjugate())
if tmp > tol:
return False
return True
@cython.boundscheck(False)
@cython.wraparound(False)
cdef int _is_hermitian_3d(
complex[:, :, :] a, double atol=1e-300, double rtol=1e-13
) except -1:
"Return True if 'a' is Hermitian"
if a.shape[1] != a.shape[2]:
return False
# compute max(a)
cdef double tmp, max_a = 0
cdef gint i, j, k
for k in range(a.shape[0]):
for i in range(a.shape[1]):
for j in range(a.shape[2]):
tmp = cabs(a[k, i, j])
if tmp > max_a:
max_a = tmp
max_a = math.sqrt(max_a)
cdef double tol = rtol * max_a + atol
for k in range(a.shape[0]):
for i in range(a.shape[1]):
for j in range(i, a.shape[2]):
tmp = cabs(a[k, i, j] - a[k, j, i].conjugate())
if tmp > tol:
return False
return True
@cython.boundscheck(False)
@cython.wraparound(False)
cdef _select(gint[:, :] arr, gint[:] indexes):
ret = np.empty((indexes.shape[0], arr.shape[1]), dtype=gint_dtype)
cdef gint[:, :] ret_view = ret
cdef gint i, j
for i in range(indexes.shape[0]):
for j in range(arr.shape[1]):
ret_view[i, j] = arr[indexes[i], j]
return ret
################ Helper functions
......@@ -87,22 +139,28 @@ cdef int _check_onsite(complex[:, :] M, gint norbs,
raise UserCodeError('Onsite matrix is not square')
if M.shape[0] != norbs:
raise UserCodeError(_shape_msg.format('Onsite'))
if check_hermiticity and not _is_herm_conj(M, M):
if check_hermiticity and not _is_hermitian(M):
raise ValueError(_herm_msg.format('Onsite'))
return 0
cdef int _check_ham(complex[:, :] H, ham, args, params,
gint a, gint a_norbs, gint b, gint b_norbs,
int check_hermiticity) except -1:
"Check Hamiltonian matrix for correct shape and hermiticity."
if H.shape[0] != a_norbs and H.shape[1] != b_norbs:
cdef int _check_onsites(complex[:, :, :] M, gint norbs,
int check_hermiticity) except -1:
"Check onsite matrix for correct shape and hermiticity."
if M.shape[1] != M.shape[2]:
raise UserCodeError('Onsite matrix is not square')
if M.shape[1] != norbs:
raise UserCodeError(_shape_msg.format('Onsite'))
if check_hermiticity and not _is_hermitian_3d(M):
raise ValueError(_herm_msg.format('Onsite'))
return 0
cdef int _check_hams(complex[:, :, :] H, gint to_norbs, gint from_norbs,
int check_hermiticity) except -1:
if H.shape[1] != to_norbs or H.shape[2] != from_norbs:
raise UserCodeError(_shape_msg.format('Hamiltonian'))
if check_hermiticity:
# call the "partner" element if we are not on the diagonal
H_conj = H if a == b else ta.matrix(ham(b, a, *args, params=params),
complex)
if not _is_herm_conj(H_conj, H):
if check_hermiticity and not _is_hermitian_3d(H):
raise ValueError(_herm_msg.format('Hamiltonian'))
return 0
......@@ -148,7 +206,7 @@ def _get_all_orbs(gint[:, :] where, gint[:, :] site_ranges):
def _get_tot_norbs(syst):
cdef gint _unused, tot_norbs
is_infinite_system = isinstance(syst, InfiniteSystem)
is_infinite_system = is_infinite(syst)
n_sites = syst.cell_size if is_infinite_system else syst.graph.num_nodes
_get_orbs(np.asarray(syst.site_ranges, dtype=gint_dtype),
n_sites, &tot_norbs, &_unused)
......@@ -159,34 +217,41 @@ def _normalize_site_where(syst, where):
"""Normalize the format of `where` when `where` contains sites.
If `where` is None, then all sites in the system are returned.
If it is a general iterator then it is expanded into an array. If `syst`
is a finalized Builder then `where` should contain `Site` objects,
If it is a general sequence then it is expanded into an array. If `syst`
is a finalized Builder then `where` may contain `Site` objects,
otherwise it should contain integers.
"""
if where is None:
size = (syst.cell_size
if isinstance(syst, InfiniteSystem) else syst.graph.num_nodes)
_where = list(range(size))
if is_infinite(syst):
where = list(range(syst.cell_size))
else:
where = list(range(syst.graph.num_nodes))
elif callable(where):
try:
_where = [syst.id_by_site[a] for a in filter(where, syst.sites)]
where = [syst.id_by_site[s] for s in filter(where, syst.sites)]
except AttributeError:
_where = list(filter(where, range(syst.graph.num_nodes)))
if is_infinite(syst):
where = [s for s in range(syst.cell_size) if where(s)]
else:
where = [s for s in range(syst.graph.num_nodes) if where(s)]
else:
try:
_where = list(syst.id_by_site[s] for s in where)
except AttributeError:
_where = list(where)
if any(w < 0 or w >= syst.graph.num_nodes for w in _where):
raise ValueError('`where` contains sites that are not in the '
'system.')
if isinstance(where[0], Site):
try:
where = [syst.id_by_site[s] for s in where]
except AttributeError:
raise TypeError("'where' contains Sites, but the system is not "
"a finalized Builder.")
where = np.asarray(where, dtype=gint_dtype).reshape(-1, 1)
if isinstance(syst, InfiniteSystem):
if any(w >= syst.cell_size for w in _where):
raise ValueError('Only sites in the fundamental domain may be '
'specified using `where`.')
if is_infinite(syst) and np.any(where >= syst.cell_size):
raise ValueError('Only sites in the fundamental domain may be '
'specified using `where`.')
if np.any(np.logical_or(where < 0, where >= syst.graph.num_nodes)):
raise ValueError('`where` contains sites that are not in the '
'system.')
return np.asarray(_where, dtype=gint_dtype).reshape(len(_where), 1)
return where
def _normalize_hopping_where(syst, where):
......@@ -194,61 +259,124 @@ def _normalize_hopping_where(syst, where):
If `where` is None, then all hoppings in the system are returned.
If it is a general iterator then it is expanded into an array. If `syst` is
a finalized Builder then `where` should contain pairs of `Site` objects,
a finalized Builder then `where` may contain pairs of `Site` objects,
otherwise it should contain pairs of integers.
"""
if where is None:
# we cannot extract the hoppings in the same order as they are in the
# graph while simultaneously excluding all inter-cell hoppings
if isinstance(syst, InfiniteSystem):
if is_infinite(syst):
raise ValueError('`where` must be provided when calculating '
'current in an InfiniteSystem.')
_where = list(syst.graph)
where = list(syst.graph)
elif callable(where):
if hasattr(syst, "sites"):
def idx_where(hop):
def idxwhere(hop):
a, b = hop
return where(syst.sites[a], syst.sites[b])
_where = list(filter(idx_where, syst.graph))
where = list(filter(idxwhere, syst.graph))
else:
_where = list(filter(lambda h: where(*h), syst.graph))
where = list(filter(lambda h: where(*h), syst.graph))
else:
if isinstance(where[0][0], Site):
try:
where = list((syst.id_by_site[a], syst.id_by_site[b])
for a, b in where)
except AttributeError:
raise TypeError("'where' contains Sites, but the system is not "
"a finalized Builder.")
# NOTE: if we ever have operators that contain elements that are
# not in the system graph, then we should modify this check
try:
_where = list((syst.id_by_site[a], syst.id_by_site[b])
for a, b in where)
except AttributeError:
_where = list(where)
# NOTE: if we ever have operators that contain elements that are
# not in the system graph, then we should modify this check
error = ValueError('`where` contains hoppings that are not '
'in the system.')
if any(not syst.graph.has_edge(*w) for w in where):
raise ValueError('`where` contains hoppings that are not in the '
'system.')
raise error
# If where contains: negative integers, or integers > # of sites
except (NodeDoesNotExistError, DisabledFeatureError):
raise error
if isinstance(syst, InfiniteSystem):
if any(a > syst.cell_size or b > syst.cell_size for a, b in _where):
raise ValueError('Only intra-cell hoppings may be specified '
'using `where`.')
where = np.asarray(where, dtype=gint_dtype)
return np.asarray(_where, dtype=gint_dtype)
if is_infinite(syst) and np.any(where > syst.cell_size):
raise ValueError('Only intra-cell hoppings may be specified '
'using `where`.')
return where
## These two classes are here to avoid using closures, as these will
## These four classes are here to avoid using closures, as these will
## break pickle support. These are only used inside '_normalize_onsite'.
def _raise_user_error(exc, func, vectorized):
msg = [
'Error occurred in user-supplied onsite function "{0}".',
'Did you remember to vectorize "{0}"?' if vectorized else '',
'See the upper part of the above backtrace for more information.',
]
msg = '\n'.join(line for line in msg if line).format(func.__name__)
raise UserCodeError(msg) from exc
class _FunctionalOnsite:
def __init__(self, onsite, sites):
def __init__(self, onsite, sites, site_ranges):
self.onsite = onsite
self.sites = sites
self.site_ranges = site_ranges
def __call__(self, site_id, *args):
return self.onsite(self.sites[site_id], *args)
def __call__(self, site_range, site_offsets, *args):
sites = self.sites
offset = self.site_ranges[site_range][0]
try:
ret = [self.onsite(sites[offset + i], *args) for i in site_offsets]
except Exception as exc:
_raise_user_error(exc, self.onsite, vectorized=False)
return _normalize_matrix_blocks(ret, len(site_offsets))
class _DictOnsite(_FunctionalOnsite):
class _VectorizedFunctionalOnsite:
def __call__(self, site_id, *args):
return self.onsite[self.sites[site_id].family]
def __init__(self, onsite, site_arrays):
self.onsite = onsite
self.site_arrays = site_arrays
def __call__(self, site_range, site_offsets, *args):
site_array = self.site_arrays[site_range]
tags = site_array.tags[site_offsets]
sites = SiteArray(site_array.family, tags)
try:
ret = self.onsite(sites, *args)
except Exception as exc:
_raise_user_error(exc, self.onsite, vectorized=True)
return _normalize_matrix_blocks(ret, len(site_offsets))
class _FunctionalOnsiteNoTransform:
def __init__(self, onsite, site_ranges):
self.onsite = onsite
self.site_ranges = site_ranges
def __call__(self, site_range, site_offsets, *args):
site_ids = self.site_ranges[site_range][0] + site_offsets
try:
ret = [self.onsite(id, *args) for id in site_ids]
except Exception as exc:
_raise_user_error(exc, self.onsite, vectorized=False)
return _normalize_matrix_blocks(ret, len(site_offsets))
class _DictOnsite:
def __init__(self, onsite, range_family_map):
self.onsite = onsite
self.range_family_map = range_family_map
def __call__(self, site_range, site_offsets, *args):
fam = self.range_family_map[site_range]
ret = [self.onsite[fam]] * len(site_offsets)
return _normalize_matrix_blocks(ret, len(site_offsets))
def _normalize_onsite(syst, onsite, check_hermiticity):
......@@ -262,21 +390,29 @@ def _normalize_onsite(syst, onsite, check_hermiticity):
if callable(onsite):
# make 'onsite' compatible with hamiltonian value functions
param_names = get_parameters(onsite)[1:]
try:
_onsite = _FunctionalOnsite(onsite, syst.sites)
except AttributeError:
_onsite = onsite
elif isinstance(onsite, collections.abc.Mapping):
if not hasattr(syst, 'sites'):
raise TypeError('Provide `onsite` as a value or a function for '
'systems that are not finalized Builders.')
if is_vectorized(syst):
_onsite = _VectorizedFunctionalOnsite(onsite, syst.site_arrays)
elif hasattr(syst, "sites"): # probably a non-vectorized finalized Builder
_onsite = _FunctionalOnsite(onsite, syst.sites, syst.site_ranges)
else: # generic old-style system, therefore *not* vectorized.
_onsite = _FunctionalOnsiteNoTransform(onsite, syst.site_ranges)
elif isinstance(onsite, collections.abc.Mapping):
# onsites known; immediately check for correct shape and hermiticity
for fam, _onsite in onsite.items():
_onsite = ta.matrix(_onsite, complex)
_check_onsite(_onsite, fam.norbs, check_hermiticity)
_onsite = _DictOnsite(onsite, syst.sites)
if is_vectorized(syst):
range_family_map = [arr.family for arr in syst.site_arrays]
elif not hasattr(syst, 'sites'):
raise TypeError('Provide `onsite` as a value or a function for '
'systems that are not finalized Builders.')
else:
# The last entry in 'site_ranges' is just an end marker, so we remove it
range_family_map = [syst.sites[r[0]].family for r in syst.site_ranges[:-1]]
_onsite = _DictOnsite(onsite, range_family_map)
else:
# single onsite; immediately check for correct shape and hermiticity
_onsite = ta.matrix(onsite, complex)
......@@ -287,7 +423,7 @@ def _normalize_onsite(syst, onsite, check_hermiticity):
# bottleneck, then we can add a code path for scalar onsites
max_norbs = max(norbs for (_, norbs, _) in syst.site_ranges)
_onsite = _onsite[0, 0] * ta.identity(max_norbs, complex)
elif len(set(map(itemgetter(1), syst.site_ranges[:-1]))) == 1:
elif len(set(norbs for _, norbs, _ in syst.site_ranges[:-1])) == 1:
# we have the same number of orbitals everywhere
norbs = syst.site_ranges[0][1]
if _onsite.shape[0] != norbs:
......@@ -302,6 +438,101 @@ def _normalize_onsite(syst, onsite, check_hermiticity):
return _onsite, param_names
def _make_onsite_or_hopping_terms(site_ranges, where):
terms = dict()
cdef gint[:] site_offsets_ = np.asarray(site_ranges, dtype=gint_dtype)[:, 0]
cdef gint i
if where.shape[1] == 1: # onsite
for i in range(where.shape[0]):
a = _bisect(site_offsets_, where[i, 0]) - 1
terms.setdefault((a, a), []).append(i)
else: # hopping
for i in range(where.shape[0]):
a = _bisect(site_offsets_, where[i, 0]) - 1
b = _bisect(site_offsets_, where[i, 1]) - 1
terms.setdefault((a, b), []).append(i)
return [(a, None, b) for a, b in terms.items()]
def _vectorized_make_onsite_terms(syst, where):
assert is_vectorized(syst)
assert where.shape[1] == 1
site_offsets = [r[0] for r in syst.site_ranges]
terms = {}
term_by_id = syst._onsite_term_by_site_id
for i in range(where.shape[0]):
w = where[i, 0]
terms.setdefault(term_by_id[w], []).append(i)
ret = []
for term_id, which in terms.items():
term = syst.terms[term_id]
((term_sa, _), (term_sites, _)) = syst.subgraphs[term.subgraph]
term_sites += site_offsets[term_sa]
which = np.asarray(which, dtype=gint_dtype)
sites = _select(where, which).reshape(-1)
selector = np.searchsorted(term_sites, sites)
ret.append((term_id, selector, which))
return ret
def _vectorized_make_hopping_terms(syst, where):
assert is_vectorized(syst)
assert where.shape[1] == 2
site_offsets = [r[0] for r in syst.site_ranges]
terms = {}
term_by_id = syst._hopping_term_by_edge_id
for i in range(where.shape[0]):
a, b = where[i, 0], where[i, 1]
edge = syst.graph.first_edge_id(a, b)
terms.setdefault(term_by_id[edge], []).append(i)
ret = []
dtype = np.dtype([('f0', int), ('f1', int)])
for term_id, which in terms.items():
herm_conj = term_id < 0
if herm_conj:
real_term_id = -term_id - 1
else:
real_term_id = term_id
which = np.asarray(which, dtype=gint_dtype)
# Select out the hoppings and reverse them if we are
# dealing with Hermitian conjugate hoppings
hops = _select(where, which)
if herm_conj:
hops = hops[:, ::-1]
# Force contiguous array to handle herm conj case also.
# Needs to be contiguous to cast to compound dtype
hops = np.ascontiguousarray(hops, dtype=int)
hops = hops.view(dtype).reshape(-1)
term = syst.terms[real_term_id]
# Get array of pairs
((to_sa, from_sa), term_hoppings) = syst.subgraphs[term.subgraph]
term_hoppings = term_hoppings.transpose() + (site_offsets[to_sa], site_offsets[from_sa])
term_hoppings = term_hoppings.view(dtype).reshape(-1)
selector = np.recarray.searchsorted(term_hoppings, hops)
ret.append((term_id, selector, which))
return ret
def _make_matrix_elements(evaluate_term, terms):
matrix_elements = []
for (term_id, term_selector, which) in terms:
which = np.asarray(which, dtype=gint_dtype)
data = evaluate_term(term_id, term_selector, which)
matrix_elements.append((which, data))
return matrix_elements
cdef class BlockSparseMatrix:
"""A sparse matrix stored as dense blocks.
......@@ -316,11 +547,10 @@ cdef class BlockSparseMatrix:
in the sparse matrix: ``(row_offset, col_offset)``.
block_shapes : gint[:, :]
``Nx2`` array: the shapes of each block, ``(n_rows, n_cols)``.
f : callable
evaluates matrix blocks. Has signature ``(a, n_rows, b, n_cols)``
where all the arguments are integers and
``a`` and ``b`` are the contents of ``where``. This function
must return a matrix of shape ``(n_rows, n_cols)``.
matrix_elements : sequence of pairs (where_indices, data)
'data' is a 3D complex array; a vector of matrices.
'where_indices' is a 1D array of indices for 'where';
'data[i]' should be placed at 'where[where_indices[i]]'.
Attributes
----------
......@@ -339,7 +569,7 @@ cdef class BlockSparseMatrix:
@cython.boundscheck(False)
@cython.wraparound(False)
def __init__(self, gint[:, :] where, gint[:, :] block_offsets,
gint[:, :] block_shapes, f):
gint[:, :] block_shapes, matrix_elements):
if (block_offsets.shape[0] != where.shape[0] or
block_shapes.shape[0] != where.shape[0]):
raise ValueError('Arrays should be the same length along '
......@@ -354,20 +584,19 @@ cdef class BlockSparseMatrix:
data_size += block_shapes[w, 0] * block_shapes[w, 1]
### Populate data array
self.data = np.empty((data_size,), dtype=complex)
cdef complex[:, :] mat
cdef gint i, j, off, a, b, a_norbs, b_norbs
for w in range(where.shape[0]):
off = self.data_offsets[w]
a_norbs = self.block_shapes[w, 0]
b_norbs = self.block_shapes[w, 1]
a = where[w, 0]
b = a if where.shape[1] == 1 else where[w, 1]
# call the function that gives the matrix
mat = f(a, a_norbs, b, b_norbs)
# Copy data
for i in range(a_norbs):
for j in range(b_norbs):
self.data[off + i * b_norbs + j] = mat[i, j]
cdef complex[:, :, :] data
cdef gint[:] where_indexes
cdef gint i, j, k, off, a, b, a_norbs, b_norbs
for where_indexes, data in matrix_elements:
for i in range(where_indexes.shape[0]):
w = where_indexes[i]
off = self.data_offsets[w]
a_norbs = self.block_shapes[w, 0]
b_norbs = self.block_shapes[w, 1]
# Copy data
for j in range(a_norbs):
for k in range(b_norbs):
self.data[off + j * b_norbs + k] = data[i, j, k]
cdef complex* get(self, gint block_idx):
return <complex*> &self.data[0] + self.data_offsets[block_idx]
......@@ -436,6 +665,11 @@ cdef class _LocalOperator:
'Declare the number of orbitals using the '
'`norbs` keyword argument when constructing '
'the site families (lattices).')
# TODO: Update this when it becomes clear how ND systems will be
# implemented.
if is_vectorized(syst) and is_infinite(syst):
raise TypeError('Vectorized infinite systems cannot yet be '
'used with operators.')
self.syst = syst
self.onsite, self._onsite_param_names = _normalize_onsite(
......@@ -447,6 +681,40 @@ cdef class _LocalOperator:
self._bound_onsite = None
self._bound_hamiltonian = None
# Here we pre-compute the datastructures that will enable us to evaluate
# the Hamiltonian and onsite functions in a vectorized way. Essentially
# we group the sites/hoppings in 'where' by what term of 'syst' they are
# in (for vectorized systems), or by the site family(s) (for
# non-vectorized systems). If the system is vectorized we store a list
# of triples:
#
# (term_id, term_selector, which)
#
# otherwise
#
# ((to_site_range, from_site_range), None, which)
#
# Where:
#
# 'term_id' → integer: term index, may be negative (indicates herm. conj.)
# 'term_selector' → 1D integer array: selects which elements from the
# subgraph of term number 'term_id' should be evaluated.
# 'which' → 1D integer array: selects which elements of 'where' this
# vectorized evaluation corresponds to.
# 'to/from_site_range' → integer: the site ranges that the elements of
# 'where' indexed by 'which' correspond to.
#
# Note that all sites/hoppings indicated by 'which' belong to the *same*
# pair of site families by construction. This is what allows for
# vectorized evaluation.
if is_vectorized(syst):
if self.where.shape[1] == 1:
self._terms = _vectorized_make_onsite_terms(syst, where)
else:
self._terms = _vectorized_make_hopping_terms(syst, where)
else:
self._terms = _make_onsite_or_hopping_terms(self._site_ranges, where)
@cython.embedsignature
def __call__(self, bra, ket=None, args=(), *, params=None):
r"""Return the matrix elements of the operator.
......@@ -525,12 +793,6 @@ cdef class _LocalOperator:
elif ket.shape != (tot_norbs,):
raise ValueError('ket vector is incorrect shape')
where = np.asarray(self.where)
where.setflags(write=False)
if self.where.shape[1] == 1:
# if `where` just contains sites, then we want a strictly 1D array
where = where.reshape(-1)
result = np.zeros((self.where.shape[0],), dtype=complex)
self._operate(out_data=result, bra=bra, ket=ket, args=args,
params=params, op=MAT_ELS)
......@@ -648,9 +910,10 @@ cdef class _LocalOperator:
"""Evaluate the onsite matrices on all elements of `where`"""
assert callable(self.onsite)
assert not (args and params)
matrix = ta.matrix
onsite = self.onsite
check_hermiticity = self.check_hermiticity
syst = self.syst
_is_vectorized = is_vectorized(syst)
if params:
try:
......@@ -662,34 +925,93 @@ cdef class _LocalOperator:
', '.join(map('"{}"'.format, missing)))
raise TypeError(''.join(msg))
def get_onsite(a, a_norbs, b, b_norbs):
mat = matrix(onsite(a, *args), complex)
_check_onsite(mat, a_norbs, check_hermiticity)
return mat
# Evaluate many onsites at once. See _LocalOperator.__init__
# for an explanation the parameters.
def eval_onsite(term_id, term_selector, which):
if _is_vectorized:
if term_id >= 0:
(sr, _), _ = syst.subgraphs[syst.terms[term_id].subgraph]
else:
(_, sr), _ = syst.subgraphs[syst.terms[-term_id - 1].subgraph]
else:
sr, _ = term_id
start_site, norbs, _ = self.syst.site_ranges[sr]
# All sites selected by 'which' are part of the same site family.
site_offsets = _select(self.where, which)[:, 0] - start_site
data = self.onsite(sr, site_offsets, *args)
_check_onsites(data, norbs, self.check_hermiticity)
return data
matrix_elements = _make_matrix_elements(eval_onsite, self._terms)
offsets, norbs = _get_all_orbs(self.where, self._site_ranges)
return BlockSparseMatrix(self.where, offsets, norbs, get_onsite)
return BlockSparseMatrix(self.where, offsets, norbs, matrix_elements)
cdef BlockSparseMatrix _eval_hamiltonian(self, args, params):
"""Evaluate the Hamiltonian on all elements of `where`."""
matrix = ta.matrix
hamiltonian = self.syst.hamiltonian
where = self.where
syst = self.syst
is_onsite = self.where.shape[1] == 1
check_hermiticity = self.check_hermiticity
def get_ham(a, a_norbs, b, b_norbs):
mat = matrix(hamiltonian(a, b, *args, params=params), complex)
_check_ham(mat, hamiltonian, args, params,
a, a_norbs, b, b_norbs, check_hermiticity)
return mat
if is_vectorized(self.syst):
# Evaluate many Hamiltonian elements at once.
# See _LocalOperator.__init__ for an explanation the parameters.
def eval_hamiltonian(term_id, term_selector, which):
herm_conj = term_id < 0
assert not is_onsite or (is_onsite and not herm_conj) # onsite terms are never hermitian conjugate
if herm_conj:
term_id = -term_id - 1
data = syst.hamiltonian_term(term_id, term_selector,
args=args, params=params)
if herm_conj:
data = data.conjugate().transpose(0, 2, 1)
# Checks for data consistency
(to_sr, from_sr), _ = syst.subgraphs[syst.terms[term_id].subgraph]
to_norbs = syst.site_ranges[to_sr][1]
from_norbs = syst.site_ranges[from_sr][1]
if herm_conj:
to_norbs, from_norbs = from_norbs, to_norbs
_check_hams(data, to_norbs, from_norbs, is_onsite and check_hermiticity)
return data
offsets, norbs = _get_all_orbs(self.where, self._site_ranges)
return BlockSparseMatrix(self.where, offsets, norbs, get_ham)
else:
# Evaluate many Hamiltonian elements at once.
# See _LocalOperator.__init__ for an explanation the parameters.
def eval_hamiltonian(term_id, term_selector, which):
if is_onsite:
data = [
syst.hamiltonian(where[i, 0], where[i, 0], *args, params=params)
for i in which
]
else:
data = [
syst.hamiltonian(where[i, 0], where[i, 1], *args, params=params)
for i in which
]
data = _normalize_matrix_blocks(data, len(which))
# Checks for data consistency
(to_sr, from_sr) = term_id
to_norbs = syst.site_ranges[to_sr][1]
from_norbs = syst.site_ranges[from_sr][1]
_check_hams(data, to_norbs, from_norbs, is_onsite and check_hermiticity)
return data
matrix_elements = _make_matrix_elements(eval_hamiltonian, self._terms)
offsets, norbs = _get_all_orbs(where, self._site_ranges)
return BlockSparseMatrix(where, offsets, norbs, matrix_elements)
def __getstate__(self):
return (
(self.check_hermiticity, self.sum),
(self.syst, self.onsite, self._onsite_param_names),
tuple(map(np.asarray, (self.where, self._site_ranges))),
(self._terms,),
(self._bound_onsite, self._bound_hamiltonian),
)
......@@ -697,6 +1019,7 @@ cdef class _LocalOperator:
((self.check_hermiticity, self.sum),
(self.syst, self.onsite, self._onsite_param_names),
(self.where, self._site_ranges),
(self._terms,),
(self._bound_onsite, self._bound_hamiltonian),
) = state
......@@ -716,10 +1039,10 @@ cdef class Density(_LocalOperator):
maps from site families to square matrices. If a function is given it
must take the same arguments as the onsite Hamiltonian functions of the
system.
where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
where : sequence of `int` or `~kwant.system.Site`, or callable, optional
Where to evaluate the operator. If ``syst`` is not a finalized Builder,
then this should be a sequence of integers. If a function is provided,
it should take a single `int` or `~kwant.builder.Site` (if ``syst`` is
it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
a finalized builder) and return True or False. If not provided, the
operator will be calculated over all sites in the system.
check_hermiticity: bool
......@@ -865,11 +1188,11 @@ cdef class Current(_LocalOperator):
matrices (scalars are allowed if the site family has 1 orbital per
site). If a function is given it must take the same arguments as the
onsite Hamiltonian functions of the system.
where : sequence of pairs of `int` or `~kwant.builder.Site`, or callable, optional
where : sequence of pairs of `int` or `~kwant.system.Site`, or callable, optional
Where to evaluate the operator. If ``syst`` is not a finalized Builder,
then this should be a sequence of pairs of integers. If a function is
provided, it should take a pair of integers or a pair of
`~kwant.builder.Site` (if ``syst`` is a finalized builder) and return
`~kwant.system.Site` (if ``syst`` is a finalized builder) and return
True or False. If not provided, the operator will be calculated over
all hoppings in the system.
check_hermiticity : bool
......@@ -997,10 +1320,10 @@ cdef class Source(_LocalOperator):
matrices (scalars are allowed if the site family has 1 orbital per
site). If a function is given it must take the same arguments as the
onsite Hamiltonian functions of the system.
where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
where : sequence of `int` or `~kwant.system.Site`, or callable, optional
Where to evaluate the operator. If ``syst`` is not a finalized Builder,
then this should be a sequence of integers. If a function is provided,
it should take a single `int` or `~kwant.builder.Site` (if ``syst`` is
it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
a finalized builder) and return True or False. If not provided, the
operator will be calculated over all sites in the system.
check_hermiticity : bool
......
kwant/physics/__init__.py
View file @ 44e42c27
......@@ -8,8 +8,14 @@
"""Physics-related algorithms"""
# Merge the public interface of all submodules.
__all__ = []
for module in ['leads', 'dispersion', 'noise', 'symmetry', 'gauge']:
exec('from . import {0}'.format(module))
exec('from .{0} import *'.format(module))
exec('__all__.extend({0}.__all__)'.format(module))
from .leads import *
from .dispersion import *
from .noise import *
from .symmetry import *
from .gauge import *
__all__ = [leads.__all__
+ dispersion.__all__
+ noise.__all__
+ symmetry.__all__
+ gauge.__all__]
kwant/physics/gauge.py
View file @ 44e42c27
# Copyright 2011-2018 Kwant authors.
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
......@@ -303,7 +303,7 @@ def _loops_in_infinite(syst):
extended_sites : callable : int -> Site
Given a site index in the extended system consisting of
two unit cells, returns the associated high-level
`kwant.builder.Site`.
`kwant.system.Site`.
"""
assert isinstance(syst, system.InfiniteSystem)
_check_infinite_syst(syst)
......@@ -375,7 +375,7 @@ def _loops_in_composite(syst):
-1 if the site is part of the reduced scattering region (see notes).
extended_sites : callable : int -> Site
Given a site index in the extended scattering region (see notes),
returns the associated high-level `kwant.builder.Site`.
returns the associated high-level `kwant.system.Site`.
Notes
-----
......@@ -401,7 +401,7 @@ def _loops_in_composite(syst):
# Get distance matrix for the extended scattering region,
# a function that maps sites to their lead patches (-1 for sites
# in the reduced scattering region), and a function that maps sites
# to high-level 'kwant.builder.Site' objects.
# to high-level 'kwant.system.Site' objects.
distance_matrix, which_patch, extended_sites =\
_extended_scattering_region(syst)
......@@ -439,7 +439,7 @@ def _extended_scattering_region(syst):
-1 if the site is part of the reduced scattering region.
extended_sites : callable : int -> Site
Given a site index in the extended scattering region, returns
the associated high-level `kwant.builder.Site`.
the associated high-level `kwant.system.Site`.
Notes
-----
......@@ -1027,7 +1027,7 @@ class magnetic_gauge:
The first callable computes the Peierls phase in the scattering
region and the remaining callables compute the Peierls phases
in each of the leads. Each callable takes a pair of
`~kwant.builder.Site` (a hopping) and returns a unit complex
`~kwant.system.Site` (a hopping) and returns a unit complex
number (Peierls phase) that multiplies that hopping.
"""
return self._peierls(syst_field, *lead_fields, tol=tol, average=False)
kwant/physics/leads.py
View file @ 44e42c27
......@@ -173,11 +173,11 @@ class StabilizedModes:
Translation eigenvectors divided by the corresponding eigenvalues.
nmodes : int
Number of left-moving (or right-moving) modes.
sqrt_hop : numpy array or None
Part of the SVD of `h_hop`, or None if the latter is invertible.
sqrt_hop : numpy array
Part of the SVD of `h_hop`.
"""
def __init__(self, vecs, vecslmbdainv, nmodes, sqrt_hop=None):
def __init__(self, vecs, vecslmbdainv, nmodes, sqrt_hop):
kwargs = locals()
kwargs.pop('self')
self.__dict__.update(kwargs)
......
kwant/physics/tests/test_dispersion.py
View file @ 44e42c27
......@@ -17,7 +17,7 @@ from math import pi, cos, sin
def make_lead():
syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
lat = kwant.lattice.square()
lat = kwant.lattice.square(norbs=1)
syst[[lat(0, 0), lat(0, 1)]] = 3
syst[lat(0, 1), lat(0, 0)] = -1
syst[((lat(1, y), lat(0, y)) for y in range(2))] = -1
......@@ -35,7 +35,7 @@ def test_band_energies(N=5):
def test_same_as_lead():
syst = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
lat = kwant.lattice.chain()
lat = kwant.lattice.chain(norbs=1)
syst[lat(0)] = 0
syst[lat(0), lat(1)] = complex(cos(0.2), sin(0.2))
......@@ -49,7 +49,7 @@ def test_same_as_lead():
def test_raise_nonhermitian():
syst = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
lat = kwant.lattice.chain()
lat = kwant.lattice.chain(norbs=1)
syst[lat(0)] = 1j
syst[lat(0), lat(1)] = complex(cos(0.2), sin(0.2))
syst = syst.finalized()
......@@ -58,7 +58,7 @@ def test_raise_nonhermitian():
def test_band_velocities():
syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
lat = kwant.lattice.square()
lat = kwant.lattice.square(norbs=1)
syst[lat(0, 0)] = 1
syst[lat(0, 1)] = 3
syst[lat(1, 0), lat(0, 0)] = -1
......@@ -75,7 +75,7 @@ def test_band_velocities():
def test_band_velocity_derivative():
syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
lat = kwant.lattice.square()
lat = kwant.lattice.square(norbs=1)
syst[lat(0, 0)] = 1
syst[lat(0, 1)] = 3
syst[lat(1, 0), lat(0, 0)] = -1
......
kwant/physics/tests/test_leads.py
View file @ 44e42c27
......@@ -267,7 +267,7 @@ def test_modes():
def test_modes_bearded_ribbon():
# Check if bearded graphene ribbons work.
lat = kwant.lattice.honeycomb()
lat = kwant.lattice.honeycomb(norbs=1)
syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
syst[lat.shape((lambda pos: -20 < pos[1] < 20),
(0, 0))] = 0.3
......@@ -417,7 +417,7 @@ def test_zero_hopping():
def make_clean_lead(W, E, t):
syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
lat = kwant.lattice.square()
lat = kwant.lattice.square(norbs=1)
syst[(lat(0, j) for j in range(W))] = E
syst[lat.neighbors()] = -t
return syst.finalized()
......
kwant/physics/tests/test_noise.py
View file @ 44e42c27
......@@ -14,7 +14,7 @@ from kwant.physics import two_terminal_shotnoise
from kwant._common import ensure_rng
n = 5
chain = kwant.lattice.chain()
chain = kwant.lattice.chain(norbs=n)
def twoterminal_system():
rng = ensure_rng(11)
......
kwant/plotter.py
View file @ 44e42c27
# -*- coding: utf-8 -*-
# Copyright 2011-2018 Kwant authors.
# Copyright 2011-2019 Kwant authors.
#
# This file is part of Kwant. It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
......@@ -402,7 +402,7 @@ def sys_leads_sites(sys, num_lead_cells=2):
Returns
-------
sites : list of (site, lead_number, copy_number) tuples
A site is a `~kwant.builder.Site` instance if the system is not finalized,
A site is a `~kwant.system.Site` instance if the system is not finalized,
and an integer otherwise. For system sites `lead_number` is `None` and
`copy_number` is `0`, for leads both are integers.
lead_cells : list of slices
......@@ -427,7 +427,7 @@ def sys_leads_sites(sys, num_lead_cells=2):
lead.builder.sites() for i in
range(num_lead_cells)))
lead_cells.append(slice(start, len(sites)))
elif isinstance(syst, system.FiniteSystem):
elif system.is_finite(syst):
sites = [(i, None, 0) for i in range(syst.graph.num_nodes)]
for leadnr, lead in enumerate(syst.leads):
start = len(sites)
......@@ -528,7 +528,7 @@ def sys_leads_hoppings(sys, num_lead_cells=2):
Returns
-------
hoppings : list of (hopping, lead_number, copy_number) tuples
A site is a `~kwant.builder.Site` instance if the system is not finalized,
A site is a `~kwant.system.Site` instance if the system is not finalized,
and an integer otherwise. For system sites `lead_number` is `None` and
`copy_number` is `0`, for leads both are integers.
lead_cells : list of slices
......@@ -735,17 +735,21 @@ def plot(sys, num_lead_cells=2, unit='nn',
symbols specifications (only for kwant.system.FiniteSystem).
site_size : number, function, array, or `None`
Relative (linear) size of the site symbol.
An array may not be used when 'syst' is a kwant.Builder.
site_color : ``matplotlib`` color description, function, array, or `None`
A color used for plotting a site in the system. If a colormap is used,
it should be a function returning single floats or a one-dimensional
array of floats. By default sites are colored by their site family,
using the current matplotlib color cycle.
An array of colors may not be used when 'syst' is a kwant.Builder.
site_edgecolor : ``matplotlib`` color description, function, array, or `None`
Color used for plotting the edges of the site symbols. Only
valid matplotlib color descriptions are allowed (and no
combination of floats and colormap as for site_color).
An array of colors may not be used when 'syst' is a kwant.Builder.
site_lw : number, function, array, or `None`
Linewidth of the site symbol edges.
An array may not be used when 'syst' is a kwant.Builder.
hop_color : ``matplotlib`` color description or a function
Same as `site_color`, but for hoppings. A function is passed two sites
in this case. (arrays are not allowed in this case).
......@@ -808,7 +812,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
- The meaning of "site" depends on whether the system to be plotted is a
builder or a low level system. For builders, a site is a
kwant.builder.Site object. For low level systems, a site is an integer
kwant.system.Site object. For low level systems, a site is an integer
-- the site number.
- color and symbol definitions may be tuples, but not lists or arrays.
......@@ -911,7 +915,11 @@ def plot(sys, num_lead_cells=2, unit='nn',
raise ValueError('Invalid value of unit argument.')
# make all specs proper: either constant or lists/np.arrays:
def make_proper_site_spec(spec, fancy_indexing=False):
def make_proper_site_spec(spec_name, spec, fancy_indexing=False):
if _p.isarray(spec) and isinstance(syst, builder.Builder):
raise TypeError('{} cannot be an array when plotting'
' a Builder; use a function instead.'
.format(spec_name))
if callable(spec):
spec = [spec(i[0]) for i in sites if i[1] is None]
if (fancy_indexing and _p.isarray(spec)
......@@ -933,7 +941,8 @@ def plot(sys, num_lead_cells=2, unit='nn',
spec = np.asarray(spec, dtype='object')
return spec
site_symbol = make_proper_site_spec(site_symbol)
site_symbol = make_proper_site_spec('site_symbol', site_symbol)
if site_symbol is None: site_symbol = defaults['site_symbol'][dim]
# separate different symbols (not done in 3D, the separation
# would mess up sorting)
......@@ -952,7 +961,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
if site_color is None:
cycle = _color_cycle()
if isinstance(syst, (builder.FiniteSystem, builder.InfiniteSystem)):
if builder.is_system(syst):
# Skipping the leads for brevity.
families = sorted({site.family for site in syst.sites})
color_mapping = dict(zip(families, cycle))
......@@ -967,10 +976,10 @@ def plot(sys, num_lead_cells=2, unit='nn',
# Unknown finalized system, no sites access.
site_color = defaults['site_color'][dim]
site_size = make_proper_site_spec(site_size, fancy_indexing)
site_color = make_proper_site_spec(site_color, fancy_indexing)
site_edgecolor = make_proper_site_spec(site_edgecolor, fancy_indexing)
site_lw = make_proper_site_spec(site_lw, fancy_indexing)
site_size = make_proper_site_spec('site_size', site_size, fancy_indexing)
site_color = make_proper_site_spec('site_color', site_color, fancy_indexing)
site_edgecolor = make_proper_site_spec('site_edgecolor', site_edgecolor, fancy_indexing)
site_lw = make_proper_site_spec('site_lw', site_lw, fancy_indexing)
hop_color = make_proper_hop_spec(hop_color)
hop_lw = make_proper_hop_spec(hop_lw)
......@@ -1282,7 +1291,7 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None,
if callable(value):
value = [value(site[0]) for site in sites]
else:
if not isinstance(syst, system.FiniteSystem):
if not system.is_finite(syst):
raise ValueError('List of values is only allowed as input '
'for finalized systems.')
value = np.array(value)
......@@ -1398,7 +1407,7 @@ def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
"for bands()")
syst = sys # for naming consistency inside function bodies
_common.ensure_isinstance(syst, system.InfiniteSystem)
_common.ensure_isinstance(syst, (system.InfiniteSystem, system.InfiniteVectorizedSystem))
momenta = np.array(momenta)
if momenta.ndim != 1:
......@@ -1474,7 +1483,7 @@ def spectrum(syst, x, y=None, params=None, mask=None, file=None,
if y is not None and not _p.has3d:
raise RuntimeError("Installed matplotlib does not support 3d plotting")
if isinstance(syst, system.FiniteSystem):
if system.is_finite(syst):
def ham(**kwargs):
return syst.hamiltonian_submatrix(params=kwargs, sparse=False)
elif callable(syst):
......@@ -1742,7 +1751,7 @@ def interpolate_current(syst, current, relwidth=None, abswidth=None, n=9):
the extents of `field`: ((x0, x1), (y0, y1), ...)
"""
if not isinstance(syst, builder.FiniteSystem):
if not builder.is_finite_system(syst):
raise TypeError("The system needs to be finalized.")
if len(current) != syst.graph.num_edges:
......@@ -1835,7 +1844,7 @@ def interpolate_density(syst, density, relwidth=None, abswidth=None, n=9,
the extents of ``field``: ((x0, x1), (y0, y1), ...)
"""
if not isinstance(syst, builder.FiniteSystem):
if not builder.is_finite_system(syst):
raise TypeError("The system needs to be finalized.")
if len(density) != len(syst.sites):
......