Coverage for kwant/builder.py : 95%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Copyright 2011-2019 Kwant authors.
2#
3# This file is part of Kwant. It is subject to the license terms in the file
4# LICENSE.rst found in the top-level directory of this distribution and at
5# https://kwant-project.org/license. A list of Kwant authors can be found in
6# the file AUTHORS.rst at the top-level directory of this distribution and at
7# https://kwant-project.org/authors.
9import abc
10import warnings
11import operator
12import collections
13import copy
14from functools import wraps, update_wrapper
15from itertools import islice, chain
16import textwrap
17import bisect
18import numbers
19import inspect
20import tinyarray as ta
21import numpy as np
22from scipy import sparse
23from . import system, graph, KwantDeprecationWarning, UserCodeError
24from .system import Site, SiteArray, SiteFamily, Symmetry, NoSymmetry
25from .linalg import lll
26from .operator import Density
27from .physics import DiscreteSymmetry, magnetic_gauge
28from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
29 interleave, deprecate_args, memoize)
32__all__ = ['Builder', 'HoppingKind', 'Lead',
33 'BuilderLead', 'SelfEnergyLead', 'ModesLead', 'add_peierls_phase']
36def validate_hopping(hopping):
37 """Verify that the argument is a valid hopping."""
39 # This check is essential to maintaining the requirement that hoppings must
40 # be tuples. Without it, list "hoppings" would work in some cases
41 # (e.g. with Builder.__contains__).
42 if not isinstance(hopping, tuple):
43 raise TypeError("{0} object is not a valid key.".format(
44 type(hopping).__name__))
46 # The following check is not strictly necessary (there would be an error
47 # anyway), but the error message would be confusing.
48 try:
49 a, b = hopping
50 except:
51 raise IndexError("Only length 2 tuples (=hoppings) are valid keys.")
53 # This check is essential for Builder.__contains__ - without it a builder
54 # would simply "not contain" such invalid hoppings. In other cases it
55 # provides a nicer error message.
56 for site in hopping:
57 if not isinstance(site, Site):
58 raise TypeError("Hopping elements must be Site objects, not {0}."
59 .format(type(site).__name__))
61 # This again is an essential check. Without it, builders would accept loop
62 # hoppings.
63 if a == b:
64 raise ValueError("A hopping connects the following site to itself:\n"
65 "{0}".format(a))
68################ Hopping kinds
70class HoppingKind(tuple):
71 """A pattern for matching hoppings.
73 An alias exists for this common name: ``kwant.HoppingKind``.
75 A hopping ``(a, b)`` matches precisely when the site family of ``a`` equals
76 `family_a` and that of ``b`` equals `family_b` and ``(a.tag - b.tag)`` is
77 equal to `delta`. In other words, the matching hoppings have the form:
78 ``(family_a(x + delta), family_b(x))``
80 Parameters
81 ----------
82 delta : Sequence of integers
83 The sequence is interpreted as a vector with integer elements.
84 family_a : `~kwant.system.SiteFamily`
85 family_b : `~kwant.system.SiteFamily` or ``None`` (default)
86 The default value means: use the same family as `family_a`.
88 Notes
89 -----
90 A ``HoppingKind`` is a callable object: When called with a
91 `~kwant.builder.Builder` as sole argument, an instance of this class will
92 return an iterator over all possible matching hoppings whose sites are
93 already present in the system. The hoppings do *not* have to be already
94 present in the system. For example::
96 kind = kwant.builder.HoppingKind((1, 0), lat)
97 syst[kind(syst)] = 1
99 Because a `~kwant.builder.Builder` can be indexed with functions or
100 iterables of functions, ``HoppingKind`` instances (or any non-tuple
101 iterables of them, e.g. a list) can be used directly as "wildcards" when
102 setting or deleting hoppings::
104 kinds = [kwant.builder.HoppingKind(v, lat) for v in [(1, 0), (0, 1)]]
105 syst[kinds] = 1
106 """
107 __slots__ = ()
109 delta = property(operator.itemgetter(0),
110 doc="The difference between the tags of the hopping's sites")
111 family_a = property(operator.itemgetter(1),
112 doc="The family of the first site in the hopping")
113 family_b = property(operator.itemgetter(2),
114 doc="The family of the second site in the hopping")
116 def __new__(cls, delta, family_a, family_b=None):
117 delta = ta.array(delta, int)
118 ensure_isinstance(family_a, SiteFamily)
119 if family_b is None:
120 family_b = family_a
121 else:
122 ensure_isinstance(family_b, SiteFamily)
123 family_b = family_b
125 try:
126 Site(family_b, family_a.normalize_tag(delta) - delta)
127 except Exception as e:
128 same_fams = family_b is family_a
129 msg = (str(family_a),
130 'and {} are'.format(family_b) if not same_fams else ' is',
131 'not compatible with delta={}'.format(delta),
132 )
133 raise ValueError(' '.join(msg)) from e
135 return tuple.__new__(cls, (delta, family_a, family_b))
137 def __call__(self, builder):
138 delta = self.delta
139 family_a = self.family_a
140 family_b = self.family_b
141 H = builder.H
142 symtofd = builder.symmetry.to_fd
144 for a in H:
145 if a.family != family_a:
146 continue
147 b = Site(family_b, a.tag - delta, True)
148 if symtofd(b) in H:
149 yield a, b
151 def __repr__(self):
152 return '{0}({1}, {2}{3})'.format(
153 self.__class__.__name__, repr(tuple(self.delta)),
154 repr(self.family_a),
155 ', ' + repr(self.family_b) if self.family_a != self.family_b else '')
157 def __str__(self):
158 return '{0}({1}, {2}{3})'.format(
159 self.__class__.__name__, tuple(self.delta),
160 self.family_a,
161 ', ' + str(self.family_b) if self.family_a != self.family_b else '')
165################ Support for Hermitian conjugation
167def herm_conj(value):
168 """
169 Calculate the hermitian conjugate of a python object.
171 If the object is neither a complex number nor a matrix, the original value
172 is returned. In the context of this module, this is the correct behavior
173 for functions.
174 """
175 if hasattr(value, 'conjugate'): 175 ↛ 186line 175 didn't jump to line 186, because the condition on line 175 was never false
176 value = value.conjugate()
177 if hasattr(value, 'shape'):
178 if len(value.shape) > 2:
179 is_ta = isinstance(value, (
180 ta.ndarray_int, ta.ndarray_float, ta.ndarray_complex))
181 value = np.swapaxes(value, -1, -2)
182 if is_ta: 182 ↛ 183line 182 didn't jump to line 183, because the condition on line 182 was never true
183 value = ta.array(value)
184 else:
185 value = value.transpose()
186 return value
189class HermConjOfFunc:
190 """Proxy returning the hermitian conjugate of the original result."""
191 __slots__ = ('function')
193 def __init__(self, function):
194 self.function = function
196 def __call__(self, i, j, *args, **kwargs):
197 return herm_conj(self.function(j, i, *args, **kwargs))
199 @property
200 def __signature__(self):
201 return inspect.signature(self.function)
204################ Leads
206class Lead(metaclass=abc.ABCMeta):
207 """Abstract base class for leads that can be attached to a `Builder`.
209 To attach a lead to a builder, append it to the builder's `~Builder.leads`
210 instance variable. See the documentation of `kwant.builder` for the
211 concrete classes of leads derived from this one.
213 Attributes
214 ----------
215 interface : sequence of sites
217 """
219 @abc.abstractmethod
220 def finalized(self):
221 """Return a finalized version of the lead.
223 Returns
224 -------
225 finalized_lead
227 Notes
228 -----
229 The finalized lead must be an object that can be used as a lead
230 in a `kwant.system.FiniteSystem`, i.e. an instance of
231 `kwant.system.InfiniteSystem`. Typically it will be a finalized
232 builder: `kwant.builder.InfiniteSystem`.
234 The order of sites for the finalized lead must be the one specified in
235 `interface`.
236 """
237 pass
240class BuilderLead(Lead):
241 """A lead made from a `Builder` with a spatial symmetry.
243 Parameters
244 ----------
245 builder : `Builder`
246 The tight-binding system of a lead. It has to possess appropriate
247 symmetry, and it may not contain hoppings between further than
248 neighboring images of the fundamental domain.
249 interface : sequence of `~kwant.system.Site` instances
250 Sequence of sites in the scattering region to which the lead is
251 attached.
253 Attributes
254 ----------
255 builder : `Builder`
256 The tight-binding system of a lead.
257 interface : list of `~kwant.system.Site` instances
258 A sorted list of interface sites.
259 padding : list of `~kwant.system.Site` instances
260 A sorted list of sites that originate from the lead, have the same
261 onsite Hamiltonian, and are connected by the same hoppings as the lead
262 sites.
264 Notes
265 -----
266 The hopping from the scattering region to the lead is assumed to be equal to
267 the hopping from a lead unit cell to the next one in the direction of the
268 symmetry vector (i.e. the lead is 'leaving' the system and starts with a
269 hopping).
271 Every system has an attribute `leads`, which stores a list of
272 `BuilderLead` objects with all the information about the leads that are
273 attached.
274 """
275 def __init__(self, builder, interface, padding=None):
276 self.builder = builder
277 self.interface = sorted(interface)
278 self.padding = sorted(padding) if padding is not None else []
280 def finalized(self):
281 """Return a `kwant.builder.InfiniteSystem` corresponding to the
282 compressed lead.
284 The order of interface sites is kept during finalization.
285 """
286 if self.builder.vectorize:
287 return InfiniteVectorizedSystem(self.builder, self.interface)
288 else:
289 return InfiniteSystem(self.builder, self.interface)
292def _ensure_signature(func):
293 """
294 Ensure that a modes/selfenergy function has a keyword-only parameter
295 ``params``, or takes ``**kwargs`` by potentially wrapping it.
296 """
297 parameters = inspect.signature(func).parameters
298 has_params = bool(parameters.get('params'))
299 has_kwargs = any(p.kind == inspect.Parameter.VAR_KEYWORD
300 for p in parameters.values())
301 if has_params or has_kwargs:
302 return func
304 # function conforming to old API: needs wrapping
305 @deprecate_args
306 def wrapper(energy, args=(), *, params=None):
307 return func(energy, args)
309 return wrapper
312class SelfEnergyLead(Lead):
313 """A general lead defined by its self energy.
315 Parameters
316 ----------
317 selfenergy_func : function
318 Has the same signature as `selfenergy` (without the ``self``
319 parameter) and returns the self energy matrix for the interface sites.
320 interface : sequence of `~kwant.system.Site` instances
321 parameters : sequence of strings
322 The parameters on which the lead depends.
323 """
324 def __init__(self, selfenergy_func, interface, parameters):
325 self.interface = tuple(interface)
326 # we changed the API of 'selfenergy_func' to have a keyword-only
327 # parameter 'params', but we still need to support the old API
328 # XXX: remove this when releasing Kwant 2.0
329 self.selfenergy_func = _ensure_signature(selfenergy_func)
330 self.parameters = frozenset(parameters)
332 def finalized(self):
333 """Trivial finalization: the object is returned itself."""
334 return self
336 @deprecate_args
337 def selfenergy(self, energy, args=(), *, params=None):
338 return self.selfenergy_func(energy, args, params=params)
341class ModesLead(Lead):
342 """A general lead defined by its modes wave functions.
344 Parameters
345 ----------
346 modes_func : function
347 Has the same signature as `modes` (without the ``self`` parameter)
348 and returns the modes of the lead as a tuple of
349 `~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes`.
350 interface :
351 sequence of `~kwant.system.Site` instances
352 parameters : sequence of strings
353 The parameters on which the lead depends.
354 """
355 def __init__(self, modes_func, interface, parameters):
356 self.interface = tuple(interface)
357 # we changed the API of 'selfenergy_func' to have a keyword-only
358 # parameter 'params', but we still need to support the old API
359 # XXX: remove this when releasing Kwant 2.0
360 self.modes_func = _ensure_signature(modes_func)
361 self.parameters = frozenset(parameters)
363 def finalized(self):
364 """Trivial finalization: the object is returned itself."""
365 return self
367 @deprecate_args
368 def modes(self, energy, args=(), *, params=None):
369 return self.modes_func(energy, args, params=params)
371 @deprecate_args
372 def selfenergy(self, energy, args=(), *, params=None):
373 stabilized = self.modes(energy, args, params=params)[1]
374 return stabilized.selfenergy()
378################ Builder class
380# A marker, meaning for hopping (i, j): this value is given by the Hermitian
381# conjugate the value of the hopping (j, i). Used by Builder and System.
382class Other:
383 pass
386def edges(seq):
387 result = interleave(seq)
388 next(result) # Skip the special loop edge.
389 return result
392def _site_ranges(sites):
393 """Return a sequence of ranges for ``sites``.
395 Here, "ranges" are defined as sequences of sites that have the same site
396 family. Because site families now have a fixed number of orbitals,
397 this coincides with the definition given in `~kwant.system.System`.
398 """
399 # we shall start a new range of different `SiteFamily`s separately,
400 # even if they happen to contain the same number of orbitals.
401 total_norbs = 0
402 current_fam = None
403 site_ranges = []
404 for idx, fam in enumerate(s.family for s in sites):
405 if not fam.norbs:
406 # can't provide site_ranges if norbs not given
407 return None
408 if fam != current_fam: # start a new run
409 current_fam = fam
410 current_norbs = fam.norbs
411 site_ranges.append((idx, current_norbs, total_norbs))
412 total_norbs += current_norbs
413 # add sentinel to the end
414 site_ranges.append((len(sites), 0, total_norbs))
415 return site_ranges
418class _Substituted:
419 """Proxy that renames function parameters."""
421 def __init__(self, func, params):
422 self.func = func
423 self.params = params
424 update_wrapper(self, func)
426 def __eq__(self, other):
427 if not isinstance(other, _Substituted): 427 ↛ 428line 427 didn't jump to line 428, because the condition on line 427 was never true
428 return False
429 return (self.func == other.func and self.params == other.params)
431 def __hash__(self):
432 return hash((self.func, self.params))
434 @property
435 def __signature__(self):
436 return inspect.Signature(
437 [inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
438 for name in self.params])
440 def __call__(self, *args):
441 return self.func(*args)
444def _substitute_params(func, subs):
445 """Substitute 'params' from 'subs' into 'func'."""
446 assert callable(func)
448 if isinstance(func, _Substituted):
449 old_params = func.params
450 old_func = func.func
451 else:
452 old_params = get_parameters(func)
453 old_func = func
455 params = tuple(subs.get(p, p) for p in old_params)
457 duplicates = [p for p, count in collections.Counter(params).items()
458 if count > 1]
459 if duplicates:
460 msg = ('Cannot rename parameters ',
461 ','.join('"{}"'.format(d) for d in duplicates),
462 ': parameters with the same name exist')
463 raise ValueError(''.join(msg))
465 if params == old_params:
466 return func
467 else:
468 return _Substituted(old_func, params)
471class Builder:
472 """A tight binding system defined on a graph.
474 An alias exists for this common name: ``kwant.Builder``.
476 This is one of the central types in Kwant. It is used to construct tight
477 binding systems in a flexible way.
479 The nodes of the graph are `~kwant.system.Site` instances. The edges,
480 i.e. the hoppings, are pairs (2-tuples) of sites. Each node and each
481 edge has a value associated with it. The values associated with nodes
482 are interpreted as on-site Hamiltonians, the ones associated with edges
483 as hopping integrals.
485 To make the graph accessible in a way that is natural within the Python
486 language it is exposed as a *mapping* (much like a built-in Python
487 dictionary). Keys are sites or hoppings. Values are 2d arrays
488 (e.g. NumPy or Tinyarray) or numbers (interpreted as 1 by 1 matrices).
490 Parameters
491 ----------
492 symmetry : `~kwant.system.Symmetry` or `None`
493 The spatial symmetry of the system.
494 conservation_law : 2D array, dictionary, function, or `None`
495 An onsite operator with integer eigenvalues that commutes with the
496 Hamiltonian. The ascending order of eigenvalues corresponds to the
497 selected ordering of the Hamiltonian subblocks. If a dict is
498 given, it maps from site families to such matrices. If a function is
499 given it must take the same arguments as the onsite Hamiltonian
500 functions of the system and return the onsite matrix.
501 time_reversal : scalar, 2D array, dictionary, function, or `None`
502 The unitary part of the onsite time-reversal symmetry operator.
503 Same format as that of `conservation_law`.
504 particle_hole : scalar, 2D array, dictionary, function, or `None`
505 The unitary part of the onsite particle-hole symmetry operator.
506 Same format as that of `conservation_law`.
507 chiral : 2D array, dictionary, function or `None`
508 The unitary part of the onsite chiral symmetry operator.
509 Same format as that of `conservation_law`.
510 vectorize : bool, default: False
511 If True then the finalized Builder will evaluate its Hamiltonian in
512 a vectorized manner. This requires that all value functions take
513 `~kwant.system.SiteArray` as first/second parameter rather than
514 `~kwant.system.Site`, and return an array of values. The returned
515 array must have the same length as the provided SiteArray(s), and
516 may contain either scalars (i.e. a vector of values) or matrices
517 (i.e. a 3d array of values).
519 Notes
520 -----
522 Values can be also functions that receive the site or the hopping (passed
523 to the function as two sites) and possibly additional arguments and are
524 expected to return a valid value. This allows to define systems quickly,
525 to modify them without reconstructing, and to save memory for many-orbital
526 models.
528 In addition to simple keys (single sites and hoppings) more powerful keys
529 are possible as well that allow to manipulate multiple sites/hoppings in a
530 single operation. Such keys are internally expanded into a sequence of
531 simple keys by using the method `Builder.expand`. For example,
532 ``syst[general_key] = value`` is equivalent to ::
534 for simple_key in syst.expand(general_key):
535 syst[simple_key] = value
537 Builder instances automatically ensure that every hopping is Hermitian, so
538 that if ``builder[a, b]`` has been set, there is no need to set
539 ``builder[b, a]``.
541 Builder instances can be made to automatically respect a `~kwant.system.Symmetry`
542 that is passed to them during creation. The behavior of builders with a symmetry
543 is slightly more sophisticated: all keys are mapped to the fundamental
544 domain of the symmetry before storing them. This may produce confusing
545 results when neighbors of a site are queried.
547 The method `attach_lead` *works* only if the sites affected by them have
548 tags which are sequences of integers. It *makes sense* only when these
549 sites live on a regular lattice, like the ones provided by `kwant.lattice`.
551 Attaching a lead manually (without the use of `~Builder.attach_lead`)
552 amounts to creating a `Lead` object and appending it to the list of leads
553 accessbile as the `~Builder.leads` attribute.
555 `conservation_law`, `time_reversal`, `particle_hole`, and `chiral`
556 affect the basis in which scattering modes derived from the builder
557 are expressed - see `~kwant.physics.DiscreteSymmetry` for details.
559 .. warning::
561 If functions are used to set values in a builder with a symmetry, then
562 they must satisfy the same symmetry. There is (currently) no check and
563 wrong results will be the consequence of a misbehaving function.
565 Attributes
566 ----------
567 leads : list of `Lead` instances
568 The leads that are attached to the system.
570 Examples
571 --------
572 Define a site.
574 >>> builder[site] = value
576 Print the value of a site.
578 >>> print(builder[site])
580 Define a hopping.
582 >>> builder[site1, site2] = value
584 Delete a site.
586 >>> del builder[site3]
588 Detach the last lead. (This does not remove the sites that were added to
589 the scattering region by `~Builder.attach_lead`.)
591 >>> del builder.leads[-1]
593 """
595 def __init__(self, symmetry=None, *, conservation_law=None, time_reversal=None,
596 particle_hole=None, chiral=None, vectorize=False):
597 if symmetry is None:
598 symmetry = NoSymmetry()
599 else:
600 ensure_isinstance(symmetry, Symmetry)
601 self.symmetry = symmetry
602 self.conservation_law = conservation_law
603 self.time_reversal = time_reversal
604 self.particle_hole = particle_hole
605 self.chiral = chiral
606 self.vectorize = vectorize
607 self.leads = []
608 self.H = {}
610 #### Note on H ####
611 #
612 # This dictionary stores a directed graph optimized for efficient querying
613 # and modification. The nodes are instances of `Site`.
614 #
615 # Each edge, specified by a ``(tail, head)`` pair of nodes, holds an object
616 # as a value. Likewise, each tail which occurs in the graph also holds a
617 # value. (Nodes which only occur as heads are not required to have
618 # values.) Every tail node has to be in the fundamental domain of the
619 # builder's symmetry.
620 #
621 # For a given `tail` site, H[tail] is a list alternately storing
622 # heads and values. (The heads occupy even locations followed by the
623 # values at odd locations.) Each pair of entries thus describes a single
624 # directed edge of the graph.
625 #
626 # The first pair of entries in each list is special: it always
627 # corresponds to a loop edge. (The head is equal to the tail.) This
628 # special edge has two purposes: It is used to store the value
629 # associated with the tail node itself, and it is necessary for the
630 # method getkey_tail which helps to conserve memory by storing equal
631 # node label only once.
633 def _get_edge(self, tail, head):
634 for h, value in edges(self.H[tail]):
635 if h == head:
636 return value
638 # (tail, head) is not present in the system, but tail is.
639 if head in self.H: 639 ↛ 640line 639 didn't jump to line 640, because the condition on line 639 was never true
640 raise KeyError((tail, head))
641 else:
642 # If already head is missing, we only report this. This way the
643 # behavior is symmetric with regard to tail and head.
644 raise KeyError(head)
646 def _set_edge(self, tail, head, value):
647 hvhv = self.H[tail]
648 heads = hvhv[2::2]
649 if head in heads:
650 i = 2 + 2 * heads.index(head)
651 hvhv[i] = head
652 hvhv[i + 1] = value
653 else:
654 hvhv.append(head)
655 hvhv.append(value)
657 def _del_edge(self, tail, head):
658 hvhv = self.H[tail]
659 heads = hvhv[2::2]
661 try:
662 i = 2 + 2 * heads.index(head)
663 except ValueError:
664 # (tail, head) is not present in the system, but tail is.
665 if head in self.H: 665 ↛ 666line 665 didn't jump to line 666, because the condition on line 665 was never true
666 raise KeyError((tail, head))
667 else:
668 # If already head is missing, we only report this. This way the
669 # behavior is symmetric with regard to tail and head.
670 raise KeyError(head)
672 del hvhv[i : i + 2]
674 def _out_neighbors(self, tail):
675 hvhv = self.H[tail]
676 return islice(hvhv, 2, None, 2)
678 def _out_degree(self, tail):
679 hvhv = self.H[tail]
680 return len(hvhv) // 2 - 1
682 def __copy__(self):
683 """Shallow copy"""
684 result = object.__new__(Builder)
685 result.symmetry = self.symmetry
686 result.conservation_law = self.conservation_law
687 result.time_reversal = self.time_reversal
688 result.particle_hole = self.particle_hole
689 result.chiral = self.chiral
690 result.vectorize = self.vectorize
691 result.leads = self.leads
692 result.H = self.H
693 return result
695 # TODO: write a test for this method.
696 def reversed(self):
697 """Return a shallow copy of the builder with the symmetry reversed.
699 This method can be used to attach the same infinite system as lead from
700 two opposite sides. It requires a builder to which an infinite
701 symmetry is associated.
702 """
703 if self.leads: 703 ↛ 704line 703 didn't jump to line 704, because the condition on line 703 was never true
704 raise ValueError('System to be reversed may not have leads.')
705 result = copy.copy(self)
706 # if we don't assign a new list we will inadvertantly add leads to
707 # the reversed system if we add leads to *this* system
708 # (because we only shallow copy)
709 result.leads = []
710 result.symmetry = self.symmetry.reversed()
711 return result
713 def __bool__(self):
714 return bool(self.H)
716 def expand(self, key):
717 """
718 Expand a general key into an iterator over simple keys.
720 Parameters
721 ----------
722 key : builder key (see notes)
723 The key to be expanded
725 Notes
726 -----
727 Keys are (recursively):
728 * Simple keys: sites or 2-tuples of sites (=hoppings).
729 * Any (non-tuple) iterable of keys, e.g. a list or a generator
730 expression.
731 * Any function that returns a key when passed a builder as sole
732 argument, e.g. a `HoppingKind` instance or the function returned
733 by `~kwant.lattice.Polyatomic.shape`.
735 This method is internally used to expand the keys when getting or
736 deleting items of a builder (i.e. ``syst[key] = value`` or ``del
737 syst[key]``).
739 """
740 itr = iter((key,))
741 iter_stack = [None]
742 while iter_stack:
743 for key in itr:
744 while callable(key):
745 key = key(self)
746 if isinstance(key, tuple):
747 # Site instances are also tuples.
748 yield key
749 else:
750 iter_stack.append(itr)
751 try:
752 itr = iter(key)
753 except TypeError:
754 raise TypeError("{0} object is not a valid key."
755 .format(type(key).__name__))
756 break
757 else:
758 itr = iter_stack.pop()
760 def __getitem__(self, key):
761 """Get the value of a single site or hopping."""
762 if isinstance(key, Site):
763 site = self.symmetry.to_fd(key)
764 return self.H[site][1]
766 sym = self.symmetry
767 validate_hopping(key)
768 a, b = sym.to_fd(*key)
769 value = self._get_edge(a, b)
770 if value is Other:
771 if not sym.in_fd(b):
772 b, a = sym.to_fd(b, a)
773 assert not sym.in_fd(a)
774 value = self._get_edge(b, a)
775 if callable(value):
776 assert not isinstance(value, HermConjOfFunc)
777 value = HermConjOfFunc(value)
778 else:
779 value = herm_conj(value)
780 return value
782 def __contains__(self, key):
783 """Tell whether the system contains a site or hopping."""
784 if isinstance(key, Site):
785 key = self.symmetry.to_fd(key)
786 return key in self.H
788 validate_hopping(key)
789 a, b = self.symmetry.to_fd(*key)
790 hvhv = self.H.get(a, ())
791 return b in islice(hvhv, 2, None, 2)
793 def _set_site(self, site, value):
794 """Set a single site."""
795 if not isinstance(site, Site):
796 raise TypeError('Expecting a site, got {0} instead.'.format(type(site).__name__))
797 site = self.symmetry.to_fd(site)
798 hvhv = self.H.setdefault(site, [])
799 if hvhv:
800 hvhv[1] = value
801 else:
802 hvhv[:] = [site, value]
804 def _set_hopping(self, hopping, value):
805 """Set a single hopping."""
806 sym = self.symmetry
807 validate_hopping(hopping)
808 a, b = sym.to_fd(*hopping)
810 if sym.in_fd(b):
811 # Make sure that we do not waste space by storing multiple instances
812 # of identical sites.
813 a2 = a = self.H[a][0]
814 b2 = b = self.H[b][0]
815 else:
816 b2, a2 = sym.to_fd(b, a)
817 assert not sym.in_fd(a2)
818 if b2 not in self.H:
819 raise KeyError(b)
821 # It's important that we have verified that b2 belongs to the system.
822 # Otherwise, we risk ending up with a half-added hopping.
823 if isinstance(value, HermConjOfFunc):
824 # Avoid nested HermConjOfFunc instances.
825 self._set_edge(a, b, Other) # May raise KeyError(a).
826 self._set_edge(b2, a2, value.function) # Must succeed.
827 else:
828 self._set_edge(a, b, value) # May raise KeyError(a).
829 self._set_edge(b2, a2, Other) # Must succeed.
831 def __setitem__(self, key, value):
832 """Set a single site/hopping or a bunch of them."""
833 func = None
834 for sh in self.expand(key):
835 if func is None:
836 func = (self._set_site if isinstance(sh, Site)
837 else self._set_hopping)
838 func(sh, value)
840 def _del_site(self, site):
841 """Delete a single site and all associated hoppings."""
842 if not isinstance(site, Site):
843 raise TypeError('Expecting a site, got {0} instead.'.format(
844 type(site).__name__))
846 tfd = self.symmetry.to_fd
847 site = tfd(site)
849 out_neighbors = self._out_neighbors(site)
851 for neighbor in out_neighbors:
852 if neighbor in self.H:
853 self._del_edge(neighbor, site)
854 else:
855 assert not self.symmetry.in_fd(neighbor)
856 self._del_edge(*tfd(neighbor, site))
858 del self.H[site]
860 def _del_hopping(self, hopping):
861 """Delete a single hopping."""
862 sym = self.symmetry
863 validate_hopping(hopping)
864 a, b = sym.to_fd(*hopping)
865 self._del_edge(a, b)
867 if sym.in_fd(b):
868 self._del_edge(b, a)
869 else:
870 b, a = sym.to_fd(b, a)
871 assert not sym.in_fd(a)
872 self._del_edge(b, a)
874 def __delitem__(self, key):
875 """Delete a single site/hopping or bunch of them."""
876 func = None
877 for sh in self.expand(key):
878 if func is None:
879 func = (self._del_site if isinstance(sh, Site)
880 else self._del_hopping)
881 func(sh)
883 def eradicate_dangling(self):
884 """Keep deleting dangling sites until none are left.
886 Sites are considered as dangling when less than two hoppings
887 lead to them.
888 """
889 sites = list(site for site in self.H
890 if self._out_degree(site) < 2)
892 for site in sites:
893 # We could have deleted the site already if there was e.g. an
894 # isolated hopping.
895 if site not in self.H:
896 continue
898 while self._out_degree(site) < 2:
899 neighbors = tuple(self._out_neighbors(site))
900 if not neighbors:
901 del self.H[site]
902 break
904 (neighbor,) = neighbors
905 if neighbor in self.H:
906 # neighbor is in the fundamental domain
907 edge = (neighbor, site)
908 else:
909 edge = self.symmetry.to_fd(neighbor, site)
910 neighbor = self.symmetry.to_fd(neighbor)
911 self._del_edge(*edge)
913 del self.H[site]
914 # Neighbor could become dangling now.
915 site = neighbor
917 def __iter__(self):
918 """Return an iterator over all sites and hoppings."""
919 return chain(self.H, self.hoppings())
921 def sites(self):
922 """Return a read-only set over all sites.
924 The sites that are returned belong to the fundamental domain of the
925 `Builder` symmetry, and are not necessarily the ones that were set
926 initially (but always the equivalent ones).
927 """
928 try:
929 return self.H.keys()
930 except AttributeError:
931 return frozenset(self.H)
933 def site_value_pairs(self):
934 """Return an iterator over all (site, value) pairs."""
935 for site, hvhv in self.H.items():
936 yield site, hvhv[1]
938 def hoppings(self):
939 """Return an iterator over all Builder hoppings.
941 The hoppings that are returned belong to the fundamental domain of the
942 `Builder` symmetry, and are not necessarily the ones that were set
943 initially (but always the equivalent ones).
944 """
945 for tail, hvhv in self.H.items():
946 for head, value in edges(hvhv):
947 if value is Other:
948 continue
949 yield tail, head
951 def hopping_value_pairs(self):
952 """Return an iterator over all (hopping, value) pairs."""
953 for tail, hvhv in self.H.items():
954 for head, value in edges(hvhv):
955 if value is Other:
956 continue
957 yield (tail, head), value
959 def dangling(self):
960 """Return an iterator over all dangling sites.
962 Sites are considered as dangling when less than two hoppings
963 lead to them.
964 """
965 for site in self.H:
966 if self._out_degree(site) < 2:
967 yield site
969 def degree(self, site):
970 """Return the number of neighbors of a site."""
971 if not isinstance(site, Site):
972 raise TypeError('Expecting a site, got {0} instead.'.format(
973 type(site).__name__))
974 site = self.symmetry.to_fd(site)
975 return self._out_degree(site)
977 def neighbors(self, site):
978 """Return an iterator over all neighbors of a site.
980 Technical note: This method respects the symmetry of the builder,
981 i.e. the returned sites are really connected to the given site (and not
982 to its image in the fundamental domain).
983 """
984 if not isinstance(site, Site):
985 raise TypeError('Expecting a site, got {0} instead.'.format(
986 type(site).__name__))
987 sym = self.symmetry
988 if isinstance(sym, NoSymmetry):
989 # Optimization for common case.
990 yield from self._out_neighbors(site)
991 return
992 shift = sym.which(site)
993 site = sym.act(-shift, site)
994 for neighbor in self._out_neighbors(site):
995 yield sym.act(shift, neighbor)
997 def closest(self, pos):
998 """Return the site that is closest to the given position.
1000 This function takes into account the symmetry of the builder. It is
1001 assumed that the symmetry is a translational symmetry.
1003 This function executes in a time proportional to the number of sites,
1004 so it is not efficient for large builders. It is especially slow for
1005 builders with a symmetry, but such systems often contain only a limited
1006 number of sites.
1008 """
1009 errmsg = ("Builder.closest() requires site families that provide "
1010 "pos().\nThe following one does not:\n")
1011 sym = self.symmetry
1012 n = sym.num_directions
1014 if n:
1015 # Determine basis in real space from first site. (The result from
1016 # any site would do.)
1017 I = ta.identity(n, int)
1018 site = next(iter(self.H))
1019 space_basis = [sym.act(element, site).pos - site.pos
1020 for element in I]
1021 space_basis, transf = lll.lll(space_basis)
1022 transf = ta.array(transf.T, int)
1024 tag_basis_cache = {}
1025 dist = float('inf')
1026 result = None
1027 for site in self.H:
1028 try:
1029 site_pos = site.pos
1030 except AttributeError:
1031 raise AttributeError(errmsg + str(site.family))
1032 if n:
1033 fam = site.family
1034 tag_basis = tag_basis_cache.get(fam)
1035 if tag_basis is None:
1036 zero_site = Site(fam, ta.zeros(len(site.tag), int))
1037 tag_basis = [sym.act(element, zero_site).tag
1038 for element in I]
1039 tag_basis = ta.dot(transf, tag_basis)
1040 tag_basis_cache[fam] = tag_basis
1041 shift = lll.cvp(pos - site_pos, space_basis, 1)[0]
1042 site = Site(fam, ta.dot(shift, tag_basis) + site.tag)
1043 site_pos = site.pos
1044 d = site_pos - pos
1045 d = ta.dot(d, d)
1046 if d < dist:
1047 dist = d
1048 result = site
1049 return result
1051 def update(self, other):
1052 """Update builder from `other`.
1054 All sites and hoppings of `other`, together with their values, are
1055 written to `self`, overwriting already existing sites and hoppings.
1056 The leads of `other` are appended to the leads of the system being
1057 updated.
1059 This method requires that both builders share the same symmetry.
1060 """
1061 if (not self.symmetry.has_subgroup(other.symmetry) 1061 ↛ 1063line 1061 didn't jump to line 1063, because the condition on line 1061 was never true
1062 or not other.symmetry.has_subgroup(self.symmetry)):
1063 raise ValueError("Both builders involved in update() must have "
1064 "equal symmetries.")
1065 for site, value in other.site_value_pairs():
1066 self[site] = value
1067 for hop, value in other.hopping_value_pairs():
1068 self[hop] = value
1069 self.leads.extend(other.leads)
1071 def __iadd__(self, other):
1072 warnings.warn("The += operator of builders is deprecated. Use "
1073 "'Builder.update()' instead.", KwantDeprecationWarning,
1074 stacklevel=2)
1075 self.update(other)
1076 return self
1078 def substituted(self, **subs):
1079 """Return a copy of this Builder with modified parameter names.
1080 """
1081 # Construct the a copy of the system with new value functions.
1082 if self.leads:
1083 raise ValueError("For simplicity, 'subsituted' is limited "
1084 "to builders without leads. Use 'substituted' "
1085 "before attaching leads to avoid this error.")
1087 # Get value *functions* only
1088 onsites = list(set(
1089 onsite for _, onsite in self.site_value_pairs()
1090 if callable(onsite)))
1091 hoppings = list(set(
1092 hop for _, hop in self.hopping_value_pairs()
1093 if callable(hop)))
1095 flatten = chain.from_iterable
1097 # Get parameter names to be substituted for each function,
1098 # without the 'site' parameter(s)
1099 onsite_params = [get_parameters(v)[1:] for v in onsites]
1100 hopping_params = [get_parameters(v)[2:] for v in hoppings]
1102 system_params = set(flatten(chain(onsite_params, hopping_params)))
1103 nonexistant_params = set(subs.keys()).difference(system_params)
1104 if nonexistant_params:
1105 msg = ('Parameters {} are not used by any onsite or hopping '
1106 'value function in this system.'
1107 ).format(nonexistant_params)
1108 warnings.warn(msg, RuntimeWarning, stacklevel=2)
1110 # Precompute map from old onsite/hopping value functions to ones
1111 # with substituted parameters.
1112 value_map = {value: _substitute_params(value, subs)
1113 for value in chain(onsites, hoppings)}
1115 result = copy.copy(self)
1116 # if we don't assign a new list we will inadvertantly add leads to
1117 # the reversed system if we add leads to *this* system
1118 # (because we only shallow copy)
1119 result.leads = []
1120 # Copy the 'H' dictionary, mapping old values to new ones using
1121 # 'value_map'. If a value does not appear in the map then it means
1122 # that the old value should be used.
1123 result.H = {}
1124 for tail, hvhv in self.H.items():
1125 result.H[tail] = list(flatten(
1126 (head, value_map.get(value, value))
1127 for head, value in interleave(hvhv)))
1129 return result
1131 def fill(self, template, shape, start, *, max_sites=10**7):
1132 """Populate builder using another one as a template.
1134 Starting from one or multiple sites, traverse the graph of the template
1135 builder and copy sites and hoppings to the target builder. The
1136 traversal stops at sites that are already present in the target and on
1137 sites that are not inside the provided shape.
1139 This function takes into account translational symmetry. As such,
1140 typically the template will have a higher symmetry than the target.
1142 Newly added sites are connected by hoppings to sites that were already
1143 present. This facilitates construction of a system by a series of
1144 calls to 'fill'.
1146 Parameters
1147 ----------
1148 template : `Builder` instance
1149 The builder used as the template. The symmetry of the target builder
1150 must be a subgroup of the symmetry of the template.
1151 shape : callable
1152 A boolean function of site returning whether the site should be
1153 added to the target builder or not. The shape must be compatible
1154 with the symmetry of the target builder.
1155 start : `~kwant.system.Site` or iterable thereof or iterable of numbers
1156 The site(s) at which the the flood-fill starts. If start is an
1157 iterable of numbers, the starting site will be
1158 ``template.closest(start)``.
1159 max_sites : positive number
1160 The maximal number of sites that may be added before
1161 ``RuntimeError`` is raised. Used to prevent using up all memory.
1163 Returns
1164 -------
1165 added_sites : list of `~kwant.system.Site` that were added to the system.
1167 Raises
1168 ------
1169 ValueError
1170 If the symmetry of the target isn't a subgroup of the template
1171 symmetry.
1172 RuntimeError
1173 If more than `max_sites` sites are to be added. The target builder
1174 will be left in an unusable state.
1176 Notes
1177 -----
1178 This function uses a flood-fill algorithm. If the template builder
1179 consists of disconnected parts, the fill will stop at their boundaries.
1181 """
1182 if not max_sites > 0:
1183 raise ValueError("max_sites must be positive.")
1185 to_fd = self.symmetry.to_fd
1186 H = self.H
1187 templ_sym = template.symmetry
1189 # Check that symmetries are commensurate.
1190 if not templ_sym.has_subgroup(self.symmetry): 1190 ↛ 1191line 1190 didn't jump to line 1191, because the condition on line 1190 was never true
1191 raise ValueError("Builder symmetry is not a subgroup of the "
1192 "template symmetry")
1194 if isinstance(start, Site):
1195 start = [start]
1196 else:
1197 start = list(start)
1198 if start and not isinstance(start[0], Site):
1199 start = [template.closest(start)]
1201 if any(s not in template for s in start):
1202 warnings.warn("fill(): Some of the starting sites are "
1203 "not in the template builder.",
1204 RuntimeWarning, stacklevel=2)
1205 start = [s for s in start if s in template]
1206 if not start:
1207 return []
1209 try:
1210 # "Active" are sites (mapped to the target's FD) that have been
1211 # verified to lie inside the shape, have been added to the target
1212 # (with `None` as value), but yet without their hoppings.
1213 active = set()
1214 congested = True
1215 for s in start:
1216 s = to_fd(s)
1217 if s not in H:
1218 congested = False
1219 if shape(s):
1220 active.add(s)
1221 H.setdefault(s, [s, None])
1223 if not active:
1224 if congested:
1225 warnings.warn("fill(): The target builder already contains "
1226 "all starting sites.", RuntimeWarning,
1227 stacklevel=2)
1228 else:
1229 warnings.warn("fill(): None of the starting sites is in "
1230 "the desired shape",
1231 RuntimeWarning, stacklevel=2)
1232 return []
1234 done = []
1235 old_active = set()
1236 new_active = set()
1238 # Flood-fill on the graph. We work site by site, writing all the
1239 # outgoing edges.
1240 while active:
1241 old_active.update(active)
1243 for tail in active:
1244 done.append(tail)
1245 if len(done) > max_sites:
1246 raise RuntimeError("Maximal number of sites (max_sites "
1247 "parameter of fill()) added.")
1249 # Make an iterator over head-value-pairs.
1250 shift = templ_sym.which(tail)
1251 templ_hvhv = template.H[templ_sym.act(-shift, tail)]
1252 templ_hvhv = iter(templ_hvhv)
1253 templ_hvhv = iter(zip(templ_hvhv, templ_hvhv))
1255 hvhv = H[tail]
1256 hvhv[1] = next(templ_hvhv)[1]
1257 old_heads = hvhv[2::2]
1259 # The remaining pairs are the heads and their associated
1260 # values.
1261 for head, value in templ_hvhv:
1262 head = templ_sym.act(shift, head)
1263 head_fd = to_fd(head)
1265 if (head_fd not in old_active
1266 and head_fd not in new_active):
1267 # The 'head' site has not been filled yet.
1268 if head_fd in H:
1269 # The 'head' site exists. (It doesn't matter
1270 # whether it's in the shape or not.) Fill the
1271 # incoming edge as well to balance the hopping.
1272 other_value = template._get_edge(
1273 *templ_sym.to_fd(head, tail))
1274 self._set_edge(*to_fd(head, tail)
1275 + (other_value,))
1276 else:
1277 if not shape(head_fd):
1278 # There is no site at 'head' and it's
1279 # outside the shape.
1280 continue
1281 new_active.add(head_fd)
1282 H.setdefault(head_fd, [head_fd, None])
1284 # Fill the outgoing edge.
1285 if head in old_heads: 1285 ↛ 1286line 1285 didn't jump to line 1286, because the condition on line 1285 was never true
1286 i = 2 + 2 * old_heads.index(head)
1287 hvhv[i] = head
1288 hvhv[i + 1] = value
1289 else:
1290 hvhv.extend((head, value))
1292 old_active, active, new_active = active, new_active, old_active
1293 new_active.clear()
1294 except Exception as e:
1295 # The graph has unbalanced edges: delete it.
1296 self.H = {}
1297 # Re-raise the exception with an additional message.
1298 msg = ("All sites of this builder have been deleted because an "
1299 "exception\noccurred during the execution of fill(): "
1300 "see above.")
1301 raise RuntimeError(msg) from e
1303 return done
1305 def attach_lead(self, lead_builder, origin=None, add_cells=0):
1306 """Attach a lead to the builder, possibly adding missing sites.
1308 This method first adds sites from 'lead_builder' until the interface
1309 where the lead will attach is "smooth". Then it appends the
1310 'lead_builder' and the interface sites as a
1311 `~kwant.builder.BuilderLead` to the 'leads' of this builder.
1313 Parameters
1314 ----------
1315 lead_builder : `Builder` with 1D translational symmetry
1316 Builder of the lead which has to be attached.
1317 origin : `~kwant.system.Site`
1318 The site which should belong to a domain where the lead should
1319 begin. It is used to attach a lead inside the system, e.g. to an
1320 inner radius of a ring.
1321 add_cells : int
1322 Number of complete unit cells of the lead to be added to the system
1323 *after* the missing sites have been added.
1325 Returns
1326 -------
1327 added_sites : list of `~kwant.system.Site`
1329 Raises
1330 ------
1331 ValueError
1332 If `lead_builder` does not have proper symmetry, has hoppings with
1333 range of more than one lead unit cell, or if it is not completely
1334 interrupted by the system.
1336 Notes
1337 -----
1338 This method is not fool-proof, i.e. if it raises an error, there is
1339 no guarantee that the system stayed unaltered.
1341 The system must "interrupt" the lead that is being attached. This means
1342 that for each site in the lead that has a hopping to a neighbnoring
1343 unit cell there must be at least one site in the system that is an
1344 image of the lead site under the action of the lead's translational
1345 symmetry. In order to interrupt the lead, the system must contain
1346 sites from the same site family as the sites in the lead. Below are
1347 three examples of leads being attached to systems::
1349 Successful Successful Unsuccessful
1350 ---------- ---------- ------------
1351 Lead System Lead System Lead System
1352 x- … o x … x- …
1353 | | | |
1354 x- … o-o x- … o-o x- … o-o
1355 | | | | | | | | |
1356 x- … o-o-o x- … o-o-o x- … o-o
1358 The second case succeeds, as even though the top site has no image in
1359 the system, because the top site has no hoppings to sites in other unit
1360 cells.
1362 Sites may be added to the system when the lead is attached, so that the
1363 interface to the lead is "smooth". Below we show the system after
1364 having attached a lead. The 'x' symbols in the system indicate the
1365 added sites::
1367 Lead System Lead System
1368 x- … x-x-o x … x
1369 | | | | | |
1370 x- … x-o-o x- … x-o-o
1371 | | | | | | | |
1372 x- … o-o-o x- … o-o-o
1373 """
1374 if self.symmetry.num_directions: 1374 ↛ 1375line 1374 didn't jump to line 1375, because the condition on line 1374 was never true
1375 raise ValueError("Leads can only be attached to finite systems.")
1377 if add_cells < 0 or int(add_cells) != add_cells: 1377 ↛ 1378line 1377 didn't jump to line 1378, because the condition on line 1377 was never true
1378 raise ValueError('add_cells must be an integer >= 0.')
1380 if self.vectorize != lead_builder.vectorize: 1380 ↛ 1381line 1380 didn't jump to line 1381, because the condition on line 1380 was never true
1381 raise ValueError(
1382 "Sites of the lead were added to the scattering "
1383 "region, but only one of these systems is vectorized. "
1384 "Vectorize both systems to allow attaching leads."
1385 )
1387 sym = lead_builder.symmetry
1388 H = lead_builder.H
1390 if sym.num_directions != 1: 1390 ↛ 1391line 1390 didn't jump to line 1391, because the condition on line 1390 was never true
1391 raise ValueError('Only builders with a 1D symmetry are allowed.')
1393 try:
1394 hop_range = max(abs(sym.which(hopping[1])[0])
1395 for hopping in lead_builder.hoppings())
1396 except ValueError: # if there are no hoppings max() will raise
1397 hop_range = 0
1399 if hop_range > 1:
1400 # Automatically increase the period, potentially warn the user.
1401 new_lead = Builder(
1402 sym.subgroup((hop_range,)),
1403 conservation_law=lead_builder.conservation_law,
1404 time_reversal=lead_builder.time_reversal,
1405 particle_hole=lead_builder.particle_hole,
1406 chiral=lead_builder.chiral,
1407 vectorize=lead_builder.vectorize,
1408 )
1409 with reraise_warnings():
1410 new_lead.fill(lead_builder, lambda site: True,
1411 lead_builder.sites(), max_sites=float('inf'))
1412 lead_builder = new_lead
1413 sym = lead_builder.symmetry
1414 H = lead_builder.H
1416 if not H:
1417 raise ValueError('Lead to be attached contains no sites.')
1419 # Check if site families of the lead are present in the system (catches
1420 # a common and a hard to find bug).
1421 families = set(site.family for site in H)
1422 lead_only_families = families.copy()
1423 for site in self.H: 1423 ↛ 1428line 1423 didn't jump to line 1428, because the loop on line 1423 didn't complete
1424 lead_only_families.discard(site.family)
1425 if not lead_only_families:
1426 break
1427 else:
1428 msg = ('Sites with site families {0} do not appear in the system, '
1429 'hence the system does not interrupt the lead.')
1430 raise ValueError(msg.format(tuple(lead_only_families)))
1432 all_doms = set()
1433 for site in self.H:
1434 if site.family not in families:
1435 continue
1436 ge = sym.which(site)
1437 if sym.act(-ge, site) in H:
1438 all_doms.add(ge[0])
1440 if origin is not None:
1441 orig_dom = sym.which(origin)[0]
1442 all_doms = [dom for dom in all_doms if dom <= orig_dom]
1443 if len(all_doms) == 0:
1444 raise ValueError('Builder does not intersect with the lead,'
1445 ' this lead cannot be attached.')
1446 max_dom = max(all_doms) + add_cells
1447 min_dom = min(all_doms)
1448 del all_doms
1450 def shape(site):
1451 domain, = sym.which(site)
1452 if domain < min_dom: 1452 ↛ 1453line 1452 didn't jump to line 1453, because the condition on line 1452 was never true
1453 raise ValueError('Builder does not interrupt the lead,'
1454 ' this lead cannot be attached.')
1455 return domain <= max_dom + 1
1457 # We start flood-fill from the first domain that doesn't belong to the
1458 # system (this one is guaranteed to contain a complete unit cell of the
1459 # lead). After flood-fill we remove that domain.
1460 start = {sym.act((max_dom + 1,), site) for site in H}
1461 with reraise_warnings():
1462 all_added = self.fill(lead_builder, shape, start,
1463 max_sites=float('inf'))
1464 all_added = [site for site in all_added if site not in start]
1465 del self[start]
1467 # Calculate the interface.
1468 interface = set()
1469 for site in H:
1470 for neighbor in lead_builder.neighbors(site):
1471 neighbor = sym.act((max_dom + 1,), neighbor)
1472 if sym.which(neighbor)[0] == max_dom:
1473 interface.add(neighbor)
1475 self.leads.append(BuilderLead(lead_builder, interface, all_added))
1476 return all_added
1478 def finalized(self):
1479 """Return a finalized (=usable with solvers) copy of the system.
1481 Returns
1482 -------
1483 finalized_system : `kwant.builder.FiniteSystem`
1484 If there is no symmetry.
1485 finalized_system : `kwant.builder.InfiniteSystem`
1486 If a symmetry is present.
1488 Notes
1489 -----
1490 This method does not modify the Builder instance for which it is
1491 called.
1493 Upon finalization, it is implicitly assumed that **every** function
1494 assigned as a value to a builder with a symmetry possesses the same
1495 symmetry.
1497 Attached leads are also finalized and will be present in the finalized
1498 system to be returned.
1500 Currently, only Builder instances without or with a 1D translational
1501 `~kwant.system.Symmetry` can be finalized.
1502 """
1503 if self.symmetry.num_directions == 0:
1504 if self.vectorize:
1505 return FiniteVectorizedSystem(self)
1506 else:
1507 return FiniteSystem(self)
1508 elif self.symmetry.num_directions == 1: 1508 ↛ 1514line 1508 didn't jump to line 1514, because the condition on line 1508 was never false
1509 if self.vectorize:
1510 return InfiniteVectorizedSystem(self)
1511 else:
1512 return InfiniteSystem(self)
1513 else:
1514 raise ValueError('Currently, only builders without or with a 1D '
1515 'translational symmetry can be finalized.')
1517 # Protect novice users from confusing error messages if they
1518 # forget to finalize their Builder.
1520 @staticmethod
1521 def _require_system(*args, **kwargs):
1522 """You need a finalized system; Use Builder.finalized() first."""
1523 raise TypeError('You need a finalized system; '
1524 'use Builder.finalized() first.')
1526 hamiltonian = hamiltonian_submatrix = modes = selfenergy = \
1527 inter_cell_hopping = cell_hamiltonian = precalculated = \
1528 _require_system
1531def _add_peierls_phase(syst, peierls_parameter):
1533 @memoize
1534 def wrap_hopping(hop):
1535 def f(*args):
1536 a, b, *args, phi = args
1537 return hop(a, b, *args) * phi(a, b)
1539 params = ('_site0', '_site1')
1540 if callable(hop): 1540 ↛ 1542line 1540 didn't jump to line 1542, because the condition on line 1540 was never false
1541 params += get_parameters(hop)[2:] # cut off both site parameters
1542 params += (peierls_parameter,)
1544 f.__signature__ = inspect.Signature(
1545 inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
1546 for name in params)
1548 return f
1550 const_hoppings = {}
1552 def const_hopping(a, b, phi):
1553 return const_hoppings[a, b] * phi(a, b)
1555 # fix the signature
1556 params = ('_site0', '_site1', peierls_parameter)
1557 const_hopping.__signature__ = inspect.Signature(
1558 inspect.Parameter(name, inspect.Parameter.POSITIONAL_ONLY)
1559 for name in params)
1561 ret = Builder(
1562 syst.symmetry,
1563 # No time reversal, as magnetic fields break it
1564 conservation_law=syst.conservation_law,
1565 particle_hole=syst.particle_hole,
1566 chiral=syst.chiral,
1567 )
1568 for i, lead in enumerate(syst.leads):
1569 ret.leads.append(
1570 BuilderLead(
1571 _add_peierls_phase(
1572 lead.builder,
1573 peierls_parameter + '_lead{}'.format(i),
1574 ),
1575 lead.interface,
1576 lead.padding,
1577 )
1578 )
1580 for site, val in syst.site_value_pairs():
1581 ret[site] = val
1583 for hop, val in syst.hopping_value_pairs():
1584 if callable(val):
1585 ret[hop] = wrap_hopping(val)
1586 else:
1587 const_hoppings[hop] = val
1588 ret[hop] = const_hopping
1590 return ret
1593def add_peierls_phase(syst, peierls_parameter='phi', fix_gauge=True):
1594 """Add a Peierls phase parameter to a Builder.
1596 Parameters
1597 ----------
1598 syst : kwant.Builder
1599 The system to which to add a Peierls phase.
1600 peierls_parameter : str
1601 The name of the Peierls phase parameter to value functions.
1602 fix_gauge : bool (default: True)
1603 If True, fix the gauge of the system and return it also
1605 Returns
1606 -------
1607 syst : kwant.builder.Finitesystem or kwant.builder.InfiniteSystem
1608 A system where all hoppings are given by value functions
1609 that have an additional parameter 'phi' that should be
1610 the Peierls phase returned by `kwant.physics.magnetic_gauge`.
1611 Any leads have similar parameters named 'phi_lead0',
1612 'phi_lead1' etc.
1613 gauge : callable
1614 Only returned if 'fix_gauge' is True. Called with magnetic
1615 field(s), and returns a parameter dictionary that can be
1616 passed to the system as 'params' (see the example below).
1618 Examples
1619 --------
1620 >>> import numpy as np
1621 >>> import kwant
1622 >>>
1623 >>> syst = make_system()
1624 >>> lead = make_lead()
1625 >>> syst.attach_lead(lead)
1626 >>> syst.attach_lead(lead.reversed())
1627 >>>
1628 >>> fsyst, phase = kwant.builder.add_peierls_phase(syst)
1629 >>>
1630 >>> def B_syst(pos):
1631 >>> return np.exp(-np.sum(pos * pos))
1632 >>>
1633 >>> kwant.smatrix(fsyst, parameters=dict(t=-1, **phase(B_syst, 0, 0)))
1634 """
1636 def wrap_gauge(gauge):
1637 phase_names = [peierls_parameter]
1638 phase_names += [peierls_parameter + '_lead{}'.format(i)
1639 for i in range(len(syst.leads))]
1641 doc = textwrap.dedent("""\
1642 Return the Peierls phase for a magnetic field configuration.
1644 Parameters
1645 ----------
1646 syst_field : scalar, vector or callable
1647 The magnetic field to apply to the scattering region.
1648 If callable, takes a position and returns the
1649 magnetic field at that position. Can be a scalar if
1650 the system is 1D or 2D, otherwise must be a vector.
1651 Magnetic field is expressed in units :math:`φ₀ / l²`,
1652 where :math:`φ₀` is the magnetic flux quantum and
1653 :math:`l` is the unit of length.
1654 *lead_fields : scalar, vector or callable
1655 The magnetic fields to apply to each of the leads, in
1656 the same format as 'syst_field'. In addition, if a callable
1657 is provided, then the magnetic field must have the symmetry
1658 of the associated lead.
1659 tol : float, default: 1E-8
1660 The tolerance to which to calculate the flux through each
1661 hopping loop in the system.
1662 average : bool, default: False
1663 If True, estimate the magnetic flux through each hopping loop
1664 in the system by evaluating the magnetic field at a single
1665 position inside the loop and multiplying it by the area of the
1666 loop. If False, then ``scipy.integrate.quad`` is used to integrate
1667 the magnetic field. This parameter is only used when 'syst_field'
1668 or 'lead_fields' are callable.
1670 Returns
1671 -------
1672 phases : dict of callables
1673 To be passed directly as the parameters of a system wrapped with
1674 'kwant.builder.add_peierls_phase'.
1675 """)
1677 @wraps(gauge)
1678 def f(*args, **kwargs):
1679 phases = gauge(*args, **kwargs)
1680 # When there are no leads 'gauge' returns a single callable
1681 if not syst.leads:
1682 phases = (phases,)
1683 return dict(zip(phase_names, phases))
1685 f.__doc__ = doc
1687 return f
1689 if syst.vectorize: 1689 ↛ 1690line 1689 didn't jump to line 1690, because the condition on line 1689 was never true
1690 raise TypeError("'add_peierls_phase' does not work with "
1691 "vectorized Builders")
1693 ret = _add_peierls_phase(syst, peierls_parameter).finalized()
1695 if fix_gauge: 1695 ↛ 1698line 1695 didn't jump to line 1698, because the condition on line 1695 was never false
1696 return ret, wrap_gauge(magnetic_gauge(ret))
1697 else:
1698 return ret
1701################ Finalized systems
1703def _raise_user_error(exc, func):
1704 msg = ('Error occurred in user-supplied value function "{0}".\n'
1705 'See the upper part of the above backtrace for more information.')
1706 raise UserCodeError(msg.format(func.__name__)) from exc
1709def _translate_cons_law(cons_law):
1710 """Translate a conservation law from builder format to something that can
1711 be used to initialize operator.Density.
1712 """
1713 if callable(cons_law):
1714 @wraps(cons_law)
1715 def vals(site, *args, **kwargs):
1716 if site.family.norbs == 1:
1717 return cons_law(site, *args, **kwargs)
1718 return np.diag(np.linalg.eigvalsh(cons_law(site, *args, **kwargs)))
1720 @wraps(cons_law)
1721 def vecs(site, *args, **kwargs):
1722 if site.family.norbs == 1:
1723 return 1
1724 return np.linalg.eigh(cons_law(site, *args, **kwargs))[1]
1726 elif isinstance(cons_law, collections.abc.Mapping):
1727 vals = {family: (value if family.norbs == 1 else
1728 ta.array(np.diag(np.linalg.eigvalsh(value))))
1729 for family, value in cons_law.items()}
1730 vecs = {family: (1 if family.norbs == 1 else
1731 ta.array(np.linalg.eigh(value)[1]))
1732 for family, value in cons_law.items()}
1734 else:
1735 try:
1736 vals, vecs = np.linalg.eigh(cons_law)
1737 vals = np.diag(vals)
1738 except np.linalg.LinAlgError as e:
1739 if '0-dimensional' not in e.args[0]:
1740 raise e # skip coverage
1741 vals, vecs = cons_law, 1
1743 return vals, vecs
1746class _FinalizedBuilderMixin:
1747 """Common functionality for all finalized builders"""
1749 def _init_discrete_symmetries(self, builder):
1750 def operator(op):
1751 return Density(self, op, check_hermiticity=False)
1753 if builder.conservation_law is None:
1754 self._cons_law = None
1755 else:
1756 a = _translate_cons_law(builder.conservation_law)
1757 self._cons_law = tuple(map(operator, a))
1758 self._symmetries = tuple(None if op is None else operator(op)
1759 for op in [builder.time_reversal,
1760 builder.particle_hole,
1761 builder.chiral])
1763 def hamiltonian(self, i, j, *args, params=None):
1764 if args and params:
1765 raise TypeError("'args' and 'params' are mutually exclusive.")
1766 if i == j:
1767 value, param_names = self.onsites[i]
1768 if param_names is not None: # 'value' is callable
1769 site = self.symmetry.to_fd(self.sites[i])
1770 if params:
1771 # See body of _value_params_pair_cache().
1772 if isinstance(param_names, Exception):
1773 raise param_names
1774 args = map(params.__getitem__, param_names)
1775 try:
1776 value = value(site, *args)
1777 except Exception as exc:
1778 if isinstance(exc, KeyError) and params:
1779 missing = [p for p in param_names if p not in params]
1780 if missing: 1780 ↛ 1784line 1780 didn't jump to line 1784, because the condition on line 1780 was never false
1781 msg = ('System is missing required arguments: ',
1782 ', '.join(map('"{}"'.format, missing)))
1783 raise TypeError(''.join(msg))
1784 _raise_user_error(exc, value)
1785 else:
1786 edge_id = self.graph.first_edge_id(i, j)
1787 value, param_names = self.hoppings[edge_id]
1788 conj = value is Other
1789 if conj:
1790 i, j = j, i
1791 edge_id = self.graph.first_edge_id(i, j)
1792 value, param_names = self.hoppings[edge_id]
1793 if param_names is not None: # 'value' is callable
1794 sites = self.sites
1795 site_i, site_j = self.symmetry.to_fd(sites[i], sites[j])
1796 if params:
1797 # See body of _value_params_pair_cache().
1798 if isinstance(param_names, Exception): 1798 ↛ 1799line 1798 didn't jump to line 1799, because the condition on line 1798 was never true
1799 raise param_names
1800 args = map(params.__getitem__, param_names)
1801 try:
1802 value = value(site_i, site_j, *args)
1803 except Exception as exc:
1804 if isinstance(exc, KeyError) and params: 1804 ↛ 1805line 1804 didn't jump to line 1805, because the condition on line 1804 was never true
1805 missing = [p for p in param_names if p not in params]
1806 if missing:
1807 msg = ('System is missing required arguments: ',
1808 ', '.join(map('"{}"'.format, missing)))
1809 raise TypeError(''.join(msg))
1810 _raise_user_error(exc, value)
1811 if conj:
1812 value = herm_conj(value)
1813 return value
1815 @deprecate_args
1816 def discrete_symmetry(self, args=(), *, params=None):
1817 if self._cons_law is not None:
1818 eigvals, eigvecs = self._cons_law
1819 eigvals = eigvals.tocoo(args, params=params)
1820 if not np.allclose(eigvals.data, np.round(eigvals.data)): 1820 ↛ 1821line 1820 didn't jump to line 1821, because the condition on line 1820 was never true
1821 raise ValueError("Conservation law must have integer"
1822 " eigenvalues.")
1823 eigvals = np.round(eigvals).tocsr()
1824 # Avoid appearance of zero eigenvalues
1825 eigvals = eigvals + 0.5 * sparse.identity(eigvals.shape[0])
1826 eigvals.eliminate_zeros()
1827 eigvecs = eigvecs.tocoo(args, params=params)
1828 projectors = [eigvecs.dot(eigvals == val)
1829 for val in sorted(np.unique(eigvals.data))]
1830 else:
1831 projectors = None
1833 def evaluate(op):
1834 return None if op is None else op.tocoo(args, params=params)
1836 return DiscreteSymmetry(projectors, *(evaluate(symm) for symm in
1837 self._symmetries))
1839 def pos(self, i):
1840 return self.sites[i].pos
1843class _VectorizedFinalizedBuilderMixin(_FinalizedBuilderMixin):
1844 """Common functionality for all vectorized finalized builders
1846 Attributes
1847 ----------
1848 _term_values : sequence
1849 Each item is either an array of values (if 'param_names is None')
1850 or a vectorized value function (in which case 'param_names' is a
1851 sequence of strings: the extra parameters to the value function).
1852 _onsite_term_by_site_id : sequence of int
1853 Maps site (integer) to the term that evaluates the associated
1854 onsite matrix element.
1855 _hopping_term_by_edge_id : sequence of int
1856 Maps edge id in the graph to the term that evaluates the associated
1857 hopping matrix element. If negative, then the associated hopping is
1858 the Hermitian conjugate of another hopping; if the number stored is
1859 't' (< 0) then the associated hopping is stored in term '-t - 1'.
1860 """
1861 def hamiltonian(self, i, j, *args, params=None):
1862 warnings.warn(
1863 "Querying individual matrix elements with 'hamiltonian' is "
1864 "deprecated, and now takes O(log N) time where N is the number "
1865 "of matrix elements per hamiltonian term. Query many matrix "
1866 "elements at once using 'hamiltonian_term'.",
1867 KwantDeprecationWarning
1868 )
1869 site_offsets, _, _ = self.site_ranges.transpose()
1870 if i == j:
1871 which_term = self._onsite_term_by_site_id[i]
1872 (w, _), (off, _) = self.subgraphs[self.terms[which_term].subgraph]
1873 i_off = i - site_offsets[w]
1874 selector = np.searchsorted(off, i_off)
1875 onsite = self.hamiltonian_term(
1876 which_term, [selector], args, params=params)
1877 return onsite[0]
1878 else:
1879 edge_id = self.graph.first_edge_id(i, j)
1880 # Map sites in previous cell to fundamental domain; vectorized
1881 # infinite systems only store sites in the FD. We already know
1882 # that this hopping is between unit cells because this is encoded
1883 # in the edge_id and term id.
1884 # Using 'is_infinite' is a bit of a hack, but 'hamiltonian' is
1885 # deprecated anyway, and refactoring everything to avoid this check
1886 # is not worth it.
1887 if system.is_infinite(self): 1887 ↛ 1888line 1887 didn't jump to line 1888, because the condition on line 1887 was never true
1888 i, j = i % self.cell_size, j % self.cell_size
1889 which_term = self._hopping_term_by_edge_id[edge_id]
1890 herm_conj = which_term < 0
1891 if herm_conj:
1892 which_term = -which_term - 1
1893 term = self.terms[which_term]
1894 # To search for hopping (i, j) in hopping_terms, we need
1895 # to treat hopping_terms as a record array of integer pairs
1896 dtype = np.dtype([('f0', int), ('f1', int)])
1897 # For hermitian conjugate terms search through the
1898 # *other* term's hoppings because those are sorted.
1900 (to_w, from_w), hoppings = self.subgraphs[term.subgraph]
1901 hoppings = hoppings.transpose() # to get array-of-pairs
1902 if herm_conj:
1903 hop = (j - site_offsets[to_w], i - site_offsets[from_w])
1904 else:
1905 hop = (i - site_offsets[to_w], j - site_offsets[from_w])
1906 hop = np.array(hop, dtype=dtype)
1907 hoppings = hoppings.view(dtype).reshape(-1)
1908 selector = np.recarray.searchsorted(hoppings, hop)
1909 h = self.hamiltonian_term(
1910 which_term, [selector], args, params=params)
1911 h = h[0]
1912 if herm_conj:
1913 h = h.conjugate().transpose()
1914 return h
1916 def hamiltonian_term(self, index, selector=slice(None),
1917 args=(), params=None):
1918 if args and params:
1919 raise TypeError("'args' and 'params' are mutually exclusive.")
1920 if index < 0: 1920 ↛ 1921line 1920 didn't jump to line 1921, because the condition on line 1920 was never true
1921 raise ValueError("term indices must be non-negative")
1923 term = self.terms[index]
1924 val = self._term_values[index]
1926 if not callable(val):
1927 return val[selector]
1929 # Construct site arrays to pass to the vectorized value function.
1930 subgraph = self.subgraphs[self.terms[index].subgraph]
1931 (to_which, from_which), (to_off, from_off) = subgraph
1932 is_onsite = to_off is from_off
1933 to_off = to_off[selector]
1934 from_off = from_off[selector]
1935 assert len(to_off) == len(from_off)
1937 to_family = self.site_arrays[to_which].family
1938 to_tags = self.site_arrays[to_which].tags
1939 to_site_array = SiteArray(to_family, to_tags[to_off])
1940 site_arrays = (to_site_array,)
1941 if not is_onsite:
1942 from_family = self.site_arrays[from_which].family
1943 from_tags = self.site_arrays[from_which].tags
1944 from_site_array = self.symmetry.act(
1945 term.symmetry_element,
1946 SiteArray(from_family, from_tags[from_off])
1947 )
1948 site_arrays = (to_site_array, from_site_array)
1950 # Construct args from params
1951 if params:
1952 # There was a problem extracting parameter names from the value
1953 # function (probably an illegal signature) and we are using
1954 # keyword parameters.
1955 if self._term_errors[index] is not None:
1956 raise self._term_errors[index]
1957 try:
1958 args = [params[p] for p in term.parameters]
1959 except KeyError:
1960 missing = [p for p in term.parameters if p not in params]
1961 msg = ('System is missing required arguments: ',
1962 ', '.join(map('"{}"'.format, missing)))
1963 raise TypeError(''.join(msg))
1965 try:
1966 ham = val(*site_arrays, *args)
1967 except Exception as exc:
1968 _raise_user_error(exc, val)
1970 expected_shape = (
1971 len(to_site_array),
1972 to_family.norbs,
1973 from_family.norbs if not is_onsite else to_family.norbs,
1974 )
1975 ham = system._normalize_matrix_blocks(ham, expected_shape,
1976 calling_function=val)
1978 return ham
1981# The same (value, parameters) pair will be used for many sites/hoppings,
1982# so we cache it to avoid wasting extra memory.
1983def _value_params_pair_cache(nstrip):
1984 def get(value):
1985 entry = cache.get(id(value))
1986 if entry is None:
1987 if isinstance(value, _Substituted):
1988 entry = value.func, value.params[nstrip:]
1989 elif callable(value):
1990 try:
1991 param_names = get_parameters(value)
1992 except ValueError as ex:
1993 # The parameter names are determined and stored in advance
1994 # for future use. This has failed, but it will only turn
1995 # into a problem if user code ever uses the 'params'
1996 # mechanism. To maintain backwards compatibility, we catch
1997 # and store the exception so that it can be raised whenever
1998 # appropriate.
1999 entry = value, ex
2000 else:
2001 entry = value, param_names[nstrip:]
2002 else:
2003 # None means: value is not callable. (That's faster to check.)
2004 entry = value, None
2005 cache[id(value)] = entry
2006 return entry
2007 assert nstrip in [1, 2]
2008 cache = {}
2009 return get
2012def _make_graph(H, id_by_site):
2013 g = graph.Graph()
2014 g.num_nodes = len(id_by_site) # Some sites could not appear in any edge.
2015 for tail, hvhv in H.items():
2016 for head in islice(hvhv, 2, None, 2):
2017 if tail == head: 2017 ↛ 2018line 2017 didn't jump to line 2018, because the condition on line 2017 was never true
2018 continue
2019 g.add_edge(id_by_site[tail], id_by_site[head])
2020 return g.compressed()
2023def _finalize_leads(leads, id_by_site):
2024 #### Connect leads.
2025 finalized_leads = []
2026 lead_interfaces = []
2027 lead_paddings = []
2028 for lead_nr, lead in enumerate(leads):
2029 try:
2030 with warnings.catch_warnings(record=True) as ws:
2031 warnings.simplefilter("always")
2032 # The following line is the whole "payload" of the entire
2033 # try-block.
2034 finalized_leads.append(lead.finalized())
2035 except ValueError as e:
2036 # Re-raise the exception with an additional message.
2037 msg = 'Problem finalizing lead {0}:'.format(lead_nr)
2038 e.args = (' '.join((msg,) + e.args),)
2039 raise
2040 else:
2041 for w in ws: 2041 ↛ 2044line 2041 didn't jump to line 2044, because the loop on line 2041 never started
2042 # Re-raise any warnings with an additional message and the
2043 # proper stacklevel.
2044 w = w.message
2045 msg = 'When finalizing lead {0}:'.format(lead_nr)
2046 warnings.warn(w.__class__(' '.join((msg,) + w.args)),
2047 stacklevel=3)
2048 try:
2049 interface = [id_by_site[isite] for isite in lead.interface]
2050 except KeyError as e:
2051 msg = ("Lead {0} is attached to a site that does not "
2052 "belong to the scattering region:\n {1}")
2053 raise ValueError(msg.format(lead_nr, e.args[0]))
2055 lead_interfaces.append(np.array(interface))
2057 padding = getattr(lead, 'padding', [])
2058 # Some padding sites might have been removed after the lead was
2059 # attached. Unlike in the case of the interface, this is not a
2060 # problem.
2061 finalized_padding = [
2062 id_by_site[isite] for isite in padding if isite in id_by_site
2063 ]
2065 lead_paddings.append(np.array(finalized_padding))
2067 return finalized_leads, lead_interfaces, lead_paddings
2070class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
2071 """Finalized `Builder` with leads.
2073 Usable as input for the solvers in `kwant.solvers`.
2075 Attributes
2076 ----------
2077 sites : sequence
2078 ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
2079 to the integer-labeled site ``i`` of the low-level system. The sites
2080 are ordered first by their family and then by their tag.
2081 id_by_site : mapping
2082 The inverse of ``sites``; maps high-level `~kwant.system.Site`
2083 instances to their integer label.
2084 Satisfies ``id_by_site[sites[i]] == i``.
2085 """
2087 def __init__(self, builder):
2088 assert builder.symmetry.num_directions == 0
2090 #### Make translation tables.
2091 sites = tuple(sorted(builder.H))
2092 id_by_site = {}
2093 for site_id, site in enumerate(sites):
2094 id_by_site[site] = site_id
2096 graph = _make_graph(builder.H, id_by_site)
2098 finalized_leads, lead_interfaces, lead_paddings = _finalize_leads(
2099 builder.leads, id_by_site)
2101 # Because many onsites/hoppings share the same (value, parameter)
2102 # pairs, we keep them in a cache so that we only store a given pair
2103 # in memory *once*. This is a similar idea to interning strings.
2104 cache = _value_params_pair_cache(1)
2105 onsites = [cache(builder.H[site][1]) for site in sites]
2106 cache = _value_params_pair_cache(2)
2107 hoppings = [cache(builder._get_edge(sites[tail], sites[head]))
2108 for tail, head in graph]
2110 # Compute the union of the parameters of onsites and hoppings. Here,
2111 # 'onsites' and 'hoppings' are pairs whose second element is one of
2112 # three things:
2113 #
2114 # * a tuple of parameter names when the matrix element is a function,
2115 # * 'None' when it is a constant,
2116 # * an exception when the parameter names could not have been
2117 # determined (See body of _value_params_pair_cache()).
2118 parameters = []
2119 for _, names in chain(onsites, hoppings):
2120 if isinstance(names, Exception):
2121 parameters = None
2122 break
2123 if names:
2124 parameters.extend(names)
2125 else:
2126 parameters = frozenset(parameters)
2128 self.graph = graph
2129 self.sites = sites
2130 self.site_ranges = _site_ranges(sites)
2131 self.id_by_site = id_by_site
2132 self.hoppings = hoppings
2133 self.onsites = onsites
2134 self.parameters = parameters
2135 self.symmetry = builder.symmetry
2136 self.leads = finalized_leads
2137 self.lead_interfaces = lead_interfaces
2138 self.lead_paddings = lead_paddings
2139 self._init_discrete_symmetries(builder)
2142def _make_site_arrays(sites):
2143 tags_by_family = {}
2144 for family, tag in sites: # Sites are tuples
2145 tags_by_family.setdefault(family, []).append(tag)
2146 site_arrays = []
2147 for family, tags in sorted(tags_by_family.items()):
2148 site_arrays.append(SiteArray(family, sorted(tags)))
2149 return site_arrays
2152# Wrapper for accessing sites by their sequential number from a
2153# list of SiteArrays.
2154class _Sites:
2156 def __init__(self, site_arrays):
2157 self._site_arrays = site_arrays
2158 lengths = [0] + [len(arr.tags) for arr in site_arrays]
2159 self._start_idxs = np.cumsum(lengths)[:-1]
2161 def __len__(self):
2162 return sum(len(arr.tags) for arr in self._site_arrays)
2164 def __getitem__(self, idx):
2165 if not isinstance(idx, numbers.Integral): 2165 ↛ 2166line 2165 didn't jump to line 2166, because the condition on line 2165 was never true
2166 raise TypeError("Only individual sites may be indexed")
2167 if idx < 0 or idx >= len(self): 2167 ↛ 2168line 2167 didn't jump to line 2168, because the condition on line 2167 was never true
2168 raise IndexError(idx)
2170 which = bisect.bisect(self._start_idxs, idx) - 1
2171 arr = self._site_arrays[which]
2172 start = self._start_idxs[which]
2173 fam = arr.family
2174 tag = arr.tags[idx - start]
2175 return Site(fam, tag)
2177 def __iter__(self):
2178 for arr in self._site_arrays:
2179 for tag in arr.tags:
2180 yield Site(arr.family, tag)
2183# Wrapper for accessing the sequential number of a site, given
2184# Site object, from a list of SiteArrays.
2185class _IdBySite:
2187 def __init__(self, site_arrays):
2188 self._site_arrays = site_arrays
2189 lengths = [0] + [len(arr.tags) for arr in site_arrays]
2190 self._start_idxs = np.cumsum(lengths)[:-1]
2192 def __len__(self):
2193 return sum(len(arr.tags) for arr in self._site_arrays)
2195 def __getitem__(self, site):
2196 # treat tags as 1D sequence by defining custom dtype
2197 tag_dtype = np.asarray(site.tag).dtype
2198 dtype = np.dtype([('f{}'.format(i), tag_dtype)
2199 for i in range(len(site.tag))])
2200 site_tag = np.asarray(site.tag).view(dtype)[0]
2201 # This slightly convoluted logic is necessary because there
2202 # may be >1 site_array with the same family for InfiniteSystems.
2203 for start, arr in zip(self._start_idxs, self._site_arrays): 2203 ↛ 2212line 2203 didn't jump to line 2212, because the loop on line 2203 didn't complete
2204 if arr.family != site.family: 2204 ↛ 2205line 2204 didn't jump to line 2205, because the condition on line 2204 was never true
2205 continue
2206 tags = arr.tags.view(dtype).reshape(-1)
2207 idx_in_fam = np.recarray.searchsorted(tags, site_tag)
2208 if idx_in_fam >= len(tags) or tags[idx_in_fam] != site_tag: 2208 ↛ 2209line 2208 didn't jump to line 2209, because the condition on line 2208 was never true
2209 continue
2210 return start + idx_in_fam
2212 raise KeyError(site)
2215def _sort_term(term, value):
2216 term = np.asarray(term)
2218 if not callable(value):
2219 # Ensure that values still correspond to the correct
2220 # sites in 'term' once the latter has been sorted.
2221 value = value[term.argsort()]
2223 term.sort()
2225 return term, value
2228def _sort_hopping_term(term, value):
2229 term = np.asarray(term)
2230 # We want to sort the hopping terms in the same way
2231 # as tuples (i, j), so we need to use a record array.
2232 dtype = np.dtype([('f0', term.dtype), ('f1', term.dtype)])
2233 term_prime = term.view(dtype=dtype).reshape(-1)
2234 # _normalize_term will sort 'term' in-place via 'term_prime'
2235 _, value = _sort_term(term_prime, value)
2237 return term, value
2240def _make_onsite_terms(builder, sites, site_arrays, term_offset):
2241 # Construct onsite terms.
2242 #
2243 # onsite_subgraphs
2244 # Will contain tuples of the form described in
2245 # kwant.system.VectorizedSystem.subgraphs, but during construction
2246 # contains lists of the sites associated with each onsite term.
2247 #
2248 # onsite_term_values
2249 # Contains a callable or array of constant values for each term.
2250 #
2251 # onsite_terms
2252 # Constructed after the main loop. Contains a kwant.system.Term
2253 # tuple for each onsite term.
2254 #
2255 # _onsite_term_by_site_id
2256 # Maps the site ID to the index of the term that the site is
2257 # a part of.
2258 # lists the index of the
2259 # Hamiltonian term associated with each site/hopping. For
2260 # Hermitian conjugate hoppings "-term - 1" is stored instead.
2262 site_offsets = np.cumsum([0] + [len(sa) for sa in site_arrays])
2264 onsite_subgraphs = []
2265 onsite_term_values = []
2266 onsite_term_parameters = []
2267 # We just use the cache to handle non/callables and exceptions
2268 cache = _value_params_pair_cache(1)
2269 # maps onsite key -> onsite term number
2270 onsite_to_term_nr = collections.OrderedDict()
2271 _onsite_term_by_site_id = []
2272 for site_id, site in enumerate(sites):
2273 val = builder.H[builder.symmetry.to_fd(site)][1]
2274 const_val = not callable(val)
2275 which_site_array = bisect.bisect(site_offsets, site_id) - 1
2276 # "key" uniquely identifies an onsite term.
2277 if const_val:
2278 key = (None, which_site_array)
2279 else:
2280 key = (id(val), which_site_array)
2281 if key not in onsite_to_term_nr:
2282 # Start a new term
2283 onsite_to_term_nr[key] = len(onsite_subgraphs)
2284 onsite_subgraphs.append([])
2285 if const_val:
2286 onsite_term_values.append([])
2287 onsite_term_parameters.append(None)
2288 else:
2289 val, params = cache(val)
2290 onsite_term_values.append(val)
2291 onsite_term_parameters.append(params)
2292 onsite_subgraphs[onsite_to_term_nr[key]].append(site_id)
2293 _onsite_term_by_site_id.append(onsite_to_term_nr[key])
2294 if const_val:
2295 vals = onsite_term_values[onsite_to_term_nr[key]]
2296 vals.append(val)
2297 # Normalize any constant values and check that the shapes are consistent.
2298 onsite_term_values = [
2299 val if callable(val)
2300 else
2301 system._normalize_matrix_blocks(
2302 val,
2303 (
2304 len(val),
2305 site_arrays[sa].family.norbs,
2306 site_arrays[sa].family.norbs,
2307 ),
2308 )
2309 for (_, sa), val in
2310 zip(onsite_to_term_nr.keys(), onsite_term_values)
2311 ]
2312 # Sort the sites in each term, and also sort the values
2313 # in the same way if they are a constant (as opposed to a callable).
2314 onsite_subgraphs, onsite_term_values = zip(*(
2315 _sort_term(term, value)
2316 for term, value in
2317 zip(onsite_subgraphs, onsite_term_values)))
2318 # Store site array index and site offsets rather than sites themselves
2319 tmp = []
2320 for (_, which), s in zip(onsite_to_term_nr, onsite_subgraphs):
2321 s = s - site_offsets[which]
2322 s.flags.writeable = False
2323 tmp.append(((which, which), (s, s)))
2324 onsite_subgraphs = tmp
2325 # onsite_term_errors[i] contains an exception if the corresponding
2326 # term has a value function with an illegal signature. We only raise
2327 # the exception if we actually try to evaluate the offending term
2328 # (to maintain backwards compatibility).
2329 onsite_term_errors = [
2330 err if isinstance(err, Exception) else None
2331 for err in onsite_term_parameters
2332 ]
2334 # We store sites in the fundamental domain, so the symmetry element
2335 # is the same for all onsite terms; it's just the identity.
2336 identity = ta.array([0] * builder.symmetry.num_directions)
2338 onsite_terms = [
2339 system.Term(
2340 subgraph=term_offset + i,
2341 symmetry_element=identity,
2342 hermitian=False,
2343 parameters=(
2344 params if not isinstance(params, Exception) else None
2345 ),
2346 )
2347 for i, params in enumerate(onsite_term_parameters)
2348 ]
2349 _onsite_term_by_site_id = np.array(_onsite_term_by_site_id)
2351 return (onsite_subgraphs, onsite_terms, onsite_term_values,
2352 onsite_term_errors, _onsite_term_by_site_id)
2355def _make_hopping_terms(builder, graph, sites, site_arrays, cell_size, term_offset):
2356 # Construct the hopping terms
2357 #
2358 # The logic is the same as for the onsite terms, with the following
2359 # differences.
2360 #
2361 # _hopping_term_by_edge_id
2362 # Maps hopping edge IDs to the index of the term that the hopping
2363 # is a part of. For Hermitian conjugate hoppings "-term_index -1"
2364 # is stored instead.
2365 #
2366 # NOTE: 'graph' contains site indices >= cell_size, however 'sites'
2367 # and 'site_arrays' only contain information about the sites
2368 # in the fundamental domain.
2369 sym = builder.symmetry
2371 site_offsets = np.cumsum([0] + [len(sa) for sa in site_arrays])
2373 hopping_subgraphs = []
2374 hopping_term_values = []
2375 hopping_term_parameters = []
2376 # We just use the cache to handle non/callables and exceptions
2377 cache = _value_params_pair_cache(2)
2378 # maps hopping key -> hopping term number
2379 hopping_to_term_nr = collections.OrderedDict()
2380 _hopping_term_by_edge_id = []
2381 for tail, head in graph:
2382 # map 'tail' and 'head' so that they are both in the FD, and
2383 # assign the symmetry element to apply to 'sites[head]'
2384 # to make the actual hopping. Here we assume that the symmetry
2385 # may only be NoSymmetry or 1D TranslationalSymmetry.
2386 if isinstance(sym, NoSymmetry):
2387 tail_site = sites[tail]
2388 head_site = sites[head]
2389 sym_el = ta.array([])
2390 else:
2391 if tail < cell_size and head < cell_size:
2392 sym_el = ta.array([0])
2393 elif head >= cell_size:
2394 assert tail < cell_size
2395 head = head - cell_size
2396 sym_el = ta.array([-1])
2397 else:
2398 tail = tail - cell_size
2399 sym_el = ta.array([1])
2400 tail_site = sites[tail]
2401 head_site = sym.act(sym_el, sites[head])
2403 val = builder._get_edge(tail_site, head_site)
2404 herm_conj = val is Other
2405 if herm_conj:
2406 val = builder._get_edge(*builder.symmetry.to_fd(head_site, tail_site))
2407 const_val = not callable(val)
2409 tail_site_array = bisect.bisect(site_offsets, tail) - 1
2410 head_site_array = bisect.bisect(site_offsets, head) - 1
2411 # "key" uniquely identifies a hopping term.
2412 if const_val:
2413 key = (None, sym_el, tail_site_array, head_site_array)
2414 else:
2415 key = (id(val), sym_el, tail_site_array, head_site_array)
2416 if herm_conj:
2417 key = (key[0], -sym_el, head_site_array, tail_site_array)
2419 if key not in hopping_to_term_nr:
2420 # start a new term
2421 hopping_to_term_nr[key] = len(hopping_subgraphs)
2422 hopping_subgraphs.append([])
2423 if const_val:
2424 hopping_term_values.append([])
2425 hopping_term_parameters.append(None)
2426 else:
2427 val, params = cache(val)
2428 hopping_term_values.append(val)
2429 hopping_term_parameters.append(params)
2431 if herm_conj:
2432 # Hopping terms come after all onsite terms, so we need
2433 # to offset the term ID by that many
2434 term_id = -term_offset - hopping_to_term_nr[key] - 1
2435 _hopping_term_by_edge_id.append(term_id)
2436 else:
2437 # Hopping terms come after all onsite terms, so we need
2438 # to offset the term ID by that many
2439 term_id = term_offset + hopping_to_term_nr[key]
2440 _hopping_term_by_edge_id.append(term_id)
2441 hopping_subgraphs[hopping_to_term_nr[key]].append((tail, head))
2442 if const_val:
2443 vals = hopping_term_values[hopping_to_term_nr[key]]
2444 vals.append(val)
2446 if hopping_subgraphs:
2447 # Normalize any constant values and check that the shapes are consistent.
2448 hopping_term_values = [
2449 val if callable(val)
2450 else
2451 system._normalize_matrix_blocks(
2452 val,
2453 (
2454 len(val),
2455 site_arrays[to_sa].family.norbs,
2456 site_arrays[from_sa].family.norbs,
2457 ),
2458 )
2459 for (_, _, to_sa, from_sa), val in
2460 zip(hopping_to_term_nr.keys(), hopping_term_values)
2461 ]
2462 # Sort the hoppings in each term, and also sort the values
2463 # in the same way if they are a constant (as opposed to a callable).
2464 hopping_subgraphs, hopping_term_values = zip(*(
2465 _sort_hopping_term(term, value)
2466 for term, value in
2467 zip(hopping_subgraphs, hopping_term_values)))
2468 # Store site array index and site offsets rather than sites themselves
2469 tmp = []
2470 for (_, _, tail_which, head_which), h in zip(hopping_to_term_nr,
2471 hopping_subgraphs):
2472 start = (site_offsets[tail_which], site_offsets[head_which])
2473 # Transpose to get a pair of arrays rather than array of pairs
2474 # We use the fact that the underlying array is stored in
2475 # array-of-pairs order to search through it in 'hamiltonian'.
2476 pairs = (h - start).transpose()
2477 pairs.flags.writeable = False
2478 tmp.append(((tail_which, head_which), pairs))
2479 hopping_subgraphs = tmp
2480 # hopping_term_errors[i] contains an exception if the corresponding
2481 # term has a value function with an illegal signature. We only raise
2482 # the exception if this becomes a problem (to maintain backwards
2483 # compatibility)
2484 hopping_term_errors = [
2485 err if isinstance(err, Exception) else None
2486 for err in hopping_term_parameters
2487 ]
2488 term_symmetry_elements = [
2489 sym_el for (_, sym_el, *_) in hopping_to_term_nr.keys()
2490 ]
2491 hopping_terms = [
2492 system.Term(
2493 subgraph=term_offset + i,
2494 symmetry_element=sym_el,
2495 hermitian=True, # Builders are always Hermitian
2496 parameters=(
2497 params if not isinstance(params, Exception) else None
2498 ),
2499 )
2500 for i, (sym_el, params) in
2501 enumerate(zip(term_symmetry_elements, hopping_term_parameters))
2502 ]
2503 _hopping_term_by_edge_id = np.array(_hopping_term_by_edge_id)
2505 return (hopping_subgraphs, hopping_terms, hopping_term_values,
2506 hopping_term_errors, _hopping_term_by_edge_id)
2509class FiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.FiniteVectorizedSystem):
2510 """Finalized `Builder` with leads.
2512 Usable as input for the solvers in `kwant.solvers`.
2514 Attributes
2515 ----------
2516 site_arrays : sequence of `~kwant.system.SiteArray`
2517 The sites of the system.
2518 sites : sequence
2519 ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
2520 to the integer-labeled site ``i`` of the low-level system. The sites
2521 are ordered first by their family and then by their tag.
2522 id_by_site : mapping
2523 The inverse of ``sites``; maps high-level `~kwant.system.Site`
2524 instances to their integer label.
2525 Satisfies ``id_by_site[sites[i]] == i``.
2526 """
2528 def __init__(self, builder):
2529 assert builder.symmetry.num_directions == 0
2530 assert builder.vectorize
2532 sites = sorted(builder.H)
2533 id_by_site = {site: site_id for site_id, site in enumerate(sites)}
2535 if not all(s.family.norbs for s in sites):
2536 raise ValueError("All sites must belong to families with norbs "
2537 "specified for vectorized Builders. Specify "
2538 "norbs when creating site families.")
2540 graph = _make_graph(builder.H, id_by_site)
2542 finalized_leads, lead_interfaces, lead_paddings = _finalize_leads(
2543 builder.leads, id_by_site)
2545 del id_by_site # cleanup due to large size
2547 site_arrays = _make_site_arrays(builder.H)
2549 (onsite_subgraphs, onsite_terms, onsite_term_values,
2550 onsite_term_errors, _onsite_term_by_site_id) = _make_onsite_terms(
2551 builder, sites, site_arrays, term_offset=0)
2553 (hopping_subgraphs, hopping_terms, hopping_term_values,
2554 hopping_term_errors, _hopping_term_by_edge_id) = _make_hopping_terms(
2555 builder, graph, sites, site_arrays, len(sites),
2556 term_offset=len(onsite_terms))
2558 # Construct the combined onsite/hopping term datastructures
2559 subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
2560 terms = tuple(onsite_terms) + tuple(hopping_terms)
2561 term_values = tuple(onsite_term_values) + tuple(hopping_term_values)
2562 term_errors = tuple(onsite_term_errors) + tuple(hopping_term_errors)
2564 # Construct system parameters
2565 parameters = set()
2566 for term in terms:
2567 if term.parameters is not None:
2568 parameters.update(term.parameters)
2569 parameters = frozenset(parameters)
2571 self.site_arrays = site_arrays
2572 self.sites = _Sites(self.site_arrays)
2573 self.id_by_site = _IdBySite(self.site_arrays)
2574 self.graph = graph
2575 self.subgraphs = subgraphs
2576 self.terms = terms
2577 self._term_values = term_values
2578 self._term_errors = term_errors
2579 self._hopping_term_by_edge_id = _hopping_term_by_edge_id
2580 self._onsite_term_by_site_id = _onsite_term_by_site_id
2581 self.parameters = parameters
2582 self.symmetry = builder.symmetry
2583 self.leads = finalized_leads
2584 self.lead_interfaces = lead_interfaces
2585 self.lead_paddings = lead_paddings
2586 self._init_discrete_symmetries(builder)
2589def _make_lead_sites(builder, interface_order):
2590 #### For each site of the fundamental domain, determine whether it has
2591 #### neighbors in the previous domain or not.
2592 sym = builder.symmetry
2593 lsites_with = [] # Fund. domain sites with neighbors in prev. dom
2594 lsites_without = [] # Remaining sites of the fundamental domain
2595 for tail in builder.H: # Loop over all sites of the fund. domain.
2596 for head in builder._out_neighbors(tail):
2597 fd = sym.which(head)[0]
2598 if fd == 1:
2599 # Tail belongs to fund. domain, head to the next domain.
2600 lsites_with.append(tail)
2601 break
2602 else:
2603 # Tail is a fund. domain site not connected to prev. domain.
2604 lsites_without.append(tail)
2606 if not lsites_with:
2607 warnings.warn('Infinite system with disconnected cells.',
2608 RuntimeWarning, stacklevel=3)
2610 ### Create list of sites and a lookup table
2611 minus_one = ta.array((-1,))
2612 plus_one = ta.array((1,))
2613 if interface_order is None:
2614 # interface must be sorted
2615 interface = [sym.act(minus_one, s) for s in lsites_with]
2616 interface.sort()
2617 else:
2618 lsites_with_set = set(lsites_with)
2619 lsites_with = []
2620 interface = []
2621 if interface_order:
2622 shift = ta.array((-sym.which(interface_order[0])[0] - 1,))
2623 for shifted_iface_site in interface_order:
2624 # Shift the interface domain before the fundamental domain.
2625 # That's the right place for the interface of a lead to be, but
2626 # the sites of interface_order might live in a different
2627 # domain.
2628 iface_site = sym.act(shift, shifted_iface_site)
2629 lsite = sym.act(plus_one, iface_site)
2631 try:
2632 lsites_with_set.remove(lsite)
2633 except KeyError:
2634 if (-sym.which(shifted_iface_site)[0] - 1,) != shift:
2635 raise ValueError(
2636 'The sites in interface_order do not all '
2637 'belong to the same lead cell.')
2638 else:
2639 raise ValueError('A site in interface_order is not an '
2640 'interface site:\n' + str(iface_site))
2641 interface.append(iface_site)
2642 lsites_with.append(lsite)
2643 if lsites_with_set:
2644 raise ValueError(
2645 'interface_order did not contain all interface sites.')
2646 # `interface_order` *must* be sorted, hence `interface` should also
2647 if interface != sorted(interface):
2648 raise ValueError('Interface sites must be sorted.')
2649 del lsites_with_set
2651 return sorted(lsites_with), sorted(lsites_without), interface
2654def _make_lead_graph(builder, sites, id_by_site, cell_size):
2655 sym = builder.symmetry
2656 g = graph.Graph()
2657 g.num_nodes = len(sites) # Some sites could not appear in any edge.
2658 for tail_id, tail in enumerate(sites[:cell_size]):
2659 for head in builder._out_neighbors(tail):
2660 head_id = id_by_site.get(head)
2661 if head_id is None:
2662 # Head belongs neither to the fundamental domain nor to the
2663 # previous domain. Check that it belongs to the next
2664 # domain and ignore it otherwise as an edge corresponding
2665 # to this one has been added already or will be added.
2666 fd = sym.which(head)[0]
2667 if fd != 1: 2667 ↛ 2671line 2667 didn't jump to line 2671, because the condition on line 2667 was never false
2668 msg = ('Further-than-nearest-neighbor cells '
2669 'are connected by hopping\n{0}.')
2670 raise ValueError(msg.format((tail, head)))
2671 continue
2672 if head_id >= cell_size:
2673 # Head belongs to previous domain. The edge added here
2674 # correspond to one left out just above.
2675 g.add_edge(head_id, tail_id)
2676 g.add_edge(tail_id, head_id)
2678 return g.compressed()
2681class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
2682 """Finalized infinite system, extracted from a `Builder`.
2684 Attributes
2685 ----------
2686 sites : sequence
2687 ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
2688 to the integer-labeled site ``i`` of the low-level system.
2689 id_by_site : mapping
2690 The inverse of ``sites``; maps high-level `~kwant.system.Site`
2691 instances to their integer label.
2692 Satisfies ``id_by_site[sites[i]] == i``.
2694 Notes
2695 -----
2696 In infinite systems ``sites`` consists of 3 parts: sites in the fundamental
2697 domain (FD) with hoppings to neighboring cells, sites in the FD with no
2698 hoppings to neighboring cells, and sites in FD+1 attached to the FD by
2699 hoppings. Each of these three subsequences is individually sorted.
2700 """
2702 def __init__(self, builder, interface_order=None):
2703 """
2704 Finalize a builder instance which has to have exactly a single
2705 symmetry direction.
2707 If interface_order is not set, the order of the interface sites in the
2708 finalized system will be arbitrary. If interface_order is set to a
2709 sequence of interface sites, this order will be kept.
2710 """
2711 sym = builder.symmetry
2712 assert sym.num_directions == 1
2714 lsites_with, lsites_without, interface = _make_lead_sites(
2715 builder, interface_order)
2716 cell_size = len(lsites_with) + len(lsites_without)
2718 # we previously sorted the interface, so don't sort it again
2719 sites = lsites_with + lsites_without + interface
2720 del lsites_with
2721 del lsites_without
2722 del interface
2723 id_by_site = {}
2724 for site_id, site in enumerate(sites):
2725 id_by_site[site] = site_id
2727 graph = _make_lead_graph(builder, sites, id_by_site, cell_size)
2729 # In the following, because many onsites/hoppings share the same
2730 # (value, parameter) pairs, we keep them in 'cache' so that we only
2731 # store a given pair in memory *once*. This is like interning strings.
2733 #### Extract onsites
2734 cache = _value_params_pair_cache(1)
2735 onsites = [cache(builder.H[tail][1]) for tail in sites[:cell_size]]
2737 #### Extract hoppings.
2738 cache = _value_params_pair_cache(2)
2739 hoppings = []
2740 for tail_id, head_id in graph:
2741 tail = sites[tail_id]
2742 head = sites[head_id]
2743 if tail_id >= cell_size:
2744 # The tail belongs to the previous domain. Find the
2745 # corresponding hopping with the tail in the fund. domain.
2746 tail, head = sym.to_fd(tail, head)
2747 hoppings.append(cache(builder._get_edge(tail, head)))
2749 # Compute the union of the parameters of onsites and hoppings. Here,
2750 # 'onsites' and 'hoppings' are pairs whose second element is one of
2751 # three things:
2752 #
2753 # * a tuple of parameter names when the matrix element is a function,
2754 # * 'None' when it is a constant,
2755 # * an exception when the parameter names could not have been
2756 # determined (See body of _value_params_pair_cache()).
2757 parameters = []
2758 for _, names in chain(onsites, hoppings):
2759 if isinstance(names, Exception): 2759 ↛ 2760line 2759 didn't jump to line 2760, because the condition on line 2759 was never true
2760 parameters = None
2761 break
2762 if names:
2763 parameters.extend(names)
2764 else:
2765 parameters = frozenset(parameters)
2767 self.graph = graph
2768 self.sites = sites
2769 self.site_ranges = _site_ranges(sites)
2770 self.id_by_site = id_by_site
2771 self.hoppings = hoppings
2772 self.onsites = onsites
2773 self.parameters = parameters
2774 self.symmetry = builder.symmetry
2775 self.cell_size = cell_size
2776 self._init_discrete_symmetries(builder)
2778 def hamiltonian(self, i, j, *args, params=None):
2779 cs = self.cell_size
2780 if i == j >= cs:
2781 i -= cs
2782 j -= cs
2783 return super().hamiltonian(i, j, *args, params=params)
2786class InfiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.InfiniteVectorizedSystem):
2787 """Finalized infinite system, extracted from a `Builder`.
2789 Attributes
2790 ----------
2791 site_arrays : sequence of `~kwant.system.SiteArray`
2792 The sites of the system.
2793 sites : sequence
2794 ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
2795 to the integer-labeled site ``i`` of the low-level system.
2796 id_by_site : mapping
2797 The inverse of ``sites``; maps high-level `~kwant.system.Site`
2798 instances to their integer label.
2799 Satisfies ``id_by_site[sites[i]] == i``.
2801 Notes
2802 -----
2803 In infinite systems ``sites`` and ``site_arrays`` consists of 3 parts:
2804 sites in the fundamental domain (FD) with hoppings to neighboring cells,
2805 sites in the FD with no hoppings to neighboring cells, and sites in FD+1
2806 attached to the FD by hoppings. Each of these three subsequences is
2807 individually sorted.
2808 """
2810 def __init__(self, builder, interface_order=None):
2811 """
2812 Finalize a builder instance which has to have exactly a single
2813 symmetry direction.
2815 If interface_order is not set, the order of the interface sites in the
2816 finalized system will be arbitrary. If interface_order is set to a
2817 sequence of interface sites, this order will be kept.
2818 """
2819 sym = builder.symmetry
2820 assert sym.num_directions == 1
2821 assert builder.vectorize
2823 lsites_with, lsites_without, interface = _make_lead_sites(
2824 builder, interface_order)
2825 cell_size = len(lsites_with) + len(lsites_without)
2828 fd_sites = lsites_with + lsites_without
2829 sites = fd_sites + interface
2830 id_by_site = {}
2831 for site_id, site in enumerate(sites):
2832 id_by_site[site] = site_id
2833 # these are needed for constructing hoppings
2834 lsites_with = frozenset(lsites_with)
2835 lsites_without = frozenset(lsites_without)
2836 interface = frozenset(interface)
2838 if not all(s.family.norbs for s in sites):
2839 raise ValueError("All sites must belong to families with norbs "
2840 "specified for vectorized Builders. Specify "
2841 "norbs when creating site families.")
2843 graph = _make_lead_graph(builder, sites, id_by_site, cell_size)
2845 del id_by_site # cleanup due to large size
2847 # In order to conform to the kwant.system.InfiniteVectorizedSystem
2848 # interface we need to put the interface sites (which have hoppings
2849 # to the next unit cell) before non-interface sites.
2850 # Note that we do not *store* the interface sites of the previous
2851 # unit cell, but we still need to *handle* site indices >= cell_size
2852 # because those are used e.g. in the graph.
2853 site_arrays = (
2854 _make_site_arrays(lsites_with)
2855 + _make_site_arrays(lsites_without)
2856 )
2858 (onsite_subgraphs, onsite_terms, onsite_term_values,
2859 onsite_term_errors, _onsite_term_by_site_id) = _make_onsite_terms(
2860 builder, fd_sites, site_arrays, term_offset=0)
2862 (hopping_subgraphs, hopping_terms, hopping_term_values,
2863 hopping_term_errors, _hopping_term_by_edge_id) = _make_hopping_terms(
2864 builder, graph, fd_sites, site_arrays, cell_size,
2865 term_offset=len(onsite_terms))
2867 # Construct the combined onsite/hopping term datastructures
2868 subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
2869 terms = tuple(onsite_terms) + tuple(hopping_terms)
2870 term_values = tuple(onsite_term_values) + tuple(hopping_term_values)
2871 term_errors = tuple(onsite_term_errors) + tuple(hopping_term_errors)
2873 # Construct system parameters
2874 parameters = set()
2875 for term in terms:
2876 if term.parameters is not None:
2877 parameters.update(term.parameters)
2878 parameters = frozenset(parameters)
2880 self.site_arrays = site_arrays
2881 # 'sites' and 'id_by_site' have to be backwards compatible with the
2882 # unvectorized interface, so we have to be able to index the interface
2883 # sites in the previous unit cell also
2884 _extended_site_arrays = self.site_arrays + _make_site_arrays(interface)
2885 self.sites = _Sites(_extended_site_arrays)
2886 self.id_by_site = _IdBySite(_extended_site_arrays)
2887 self.graph = graph
2888 self.subgraphs = subgraphs
2889 self.terms = terms
2890 self._term_values = term_values
2891 self._term_errors = term_errors
2892 self._hopping_term_by_edge_id = _hopping_term_by_edge_id
2893 self._onsite_term_by_site_id = _onsite_term_by_site_id
2894 self.parameters = parameters
2895 self.symmetry = builder.symmetry
2896 self.cell_size = cell_size
2897 self._init_discrete_symmetries(builder)
2899 def hamiltonian(self, i, j, *args, params=None):
2900 cs = self.cell_size
2901 if i == j >= cs: 2901 ↛ 2902line 2901 didn't jump to line 2902, because the condition on line 2901 was never true
2902 i -= cs
2903 j -= cs
2904 return super().hamiltonian(i, j, *args, params=params)
2907def is_finite_system(syst):
2908 return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
2911def is_infinite_system(syst):
2912 return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
2915def is_system(syst):
2916 return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem,
2917 InfiniteSystem, InfiniteVectorizedSystem))