builder.py 51 KB
Newer Older
1
# Copyright 2011-2015 Kwant authors.
2
#
Christoph Groth's avatar
Christoph Groth committed
3
4
# This file is part of Kwant.  It is subject to the license terms in the file
# LICENSE.rst found in the top-level directory of this distribution and at
5
# http://kwant-project.org/license.  A list of Kwant authors can be found in
Christoph Groth's avatar
Christoph Groth committed
6
# the file AUTHORS.rst at the top-level directory of this distribution and at
7
8
# http://kwant-project.org/authors.

9
__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
10
           'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead']
11

Anton Akhmerov's avatar
Anton Akhmerov committed
12
import abc
13
import warnings
14
import operator
Joseph Weston's avatar
Joseph Weston committed
15
from itertools import islice, chain
16
import tinyarray as ta
17
import numpy as np
18
from . import system, graph, KwantDeprecationWarning
19
from ._common import ensure_isinstance
20

21

Christoph Groth's avatar
Christoph Groth committed
22

23
################ Sites and site families
24

25
class Site(tuple):
26
    """A site, member of a `SiteFamily`.
27
28
29
30

    Sites are the vertices of the graph which describes the tight binding
    system in a `Builder`.

31
    A site is uniquely identified by its family and its tag.
32
33
34

    Parameters
    ----------
35
    family : an instance of `SiteFamily`
Christoph Groth's avatar
Christoph Groth committed
36
        The 'type' of the site.
37
    tag : a hashable python object
38
        The unique identifier of the site within the site family, typically a
Christoph Groth's avatar
Christoph Groth committed
39
        vector of integers.
40
41
42
43

    Raises
    ------
    ValueError
44
        If `tag` is not a proper tag for `family`.
45
46
47

    Notes
    -----
48
    For convenience, ``family(*tag)`` can be used instead of ``Site(family,
49
50
51
52
    tag)`` to create a site.

    The parameters of the constructor (see above) are stored as instance
    variables under the same names.  Given a site ``site``, common things to
53
    query are thus ``site.family``, ``site.tag``, and ``site.pos``.
54
    """
55
    __slots__ = ()
56

57
58
59
60
    family = property(operator.itemgetter(0),
                      doc="The site family to which the site belongs.")
    tag = property(operator.itemgetter(1), doc="The tag of the site.")

61

62
    def __new__(cls, family, tag, _i_know_what_i_do=False):
Christoph Groth's avatar
Christoph Groth committed
63
        if _i_know_what_i_do:
64
            return tuple.__new__(cls, (family, tag))
65
        try:
66
            tag = family.normalize_tag(tag)
67
        except (TypeError, ValueError) as e:
68
            msg = 'Tag {0} is not allowed for site family {1}: {2}'
69
            raise type(e)(msg.format(repr(tag), repr(family), e.args[0]))
70
        return tuple.__new__(cls, (family, tag))
71
72

    def __repr__(self):
73
        return 'Site({0}, {1})'.format(repr(self.family), repr(self.tag))
74

75
    def __str__(self):
76
        sf = self.family
77
        return '<Site {0} of {1}>'.format(self.tag, sf.name if sf.name else sf)
78

79
80
81
    @property
    def pos(self):
        """Real space position of the site."""
82
        return self.family.pos(self.tag)
83
84


85
class SiteFamily(metaclass=abc.ABCMeta):
86
87
88
89
90
    """Abstract base class for site families.

    Site families are the 'type' of `Site` objects.  Within a family, individual
    sites are uniquely identified by tags.  Valid tags must be hashable Python
    objects, further details are up to the family.
91

92
93
94
95
96
97
98
    Site families must be immutable and fully defined by their initial
    arguments.  They must inherit from this abstract base class and call its
    __init__ function providing it with two arguments: a canonical
    representation and a name.  The canonical representation will be returned as
    the objects representation and must uniquely identify the site family
    instance.  The name is a string used to distinguish otherwise identical site
    families.  It may be empty.
99

100
    All site families must define the method `normalize_tag` which brings a tag
101
    to the standard format for this site family.
102

103
104
105
106
    Site families that are intended for use with plotting should also provide a
    method `pos(tag)`, which returns a vector with real-space coordinates of the
    site belonging to this family with a given tag.

107
108
    """

109
110
111
112
113
    def __init__(self, canonical_repr, name):
        self.canonical_repr = canonical_repr
        self.hash = hash(canonical_repr)
        self.name = name

114
    def __repr__(self):
115
116
117
        return self.canonical_repr

    def __str__(self):
118
119
120
121
122
        if self.name:
            msg = '<{0} site family {1}>'
        else:
            msg = '<unnamed {0} site family>'
        return msg.format(self.__class__.__name__, self.name)
123

124
    def __hash__(self):
125
        return self.hash
126
127
128
129
130
131
132

    def __eq__(self, other):
        try:
            return self.canonical_repr == other.canonical_repr
        except AttributeError:
            return False

133
    def __ne__(self, other):
134
135
136
        try:
            return self.canonical_repr != other.canonical_repr
        except AttributeError:
137
            return True
138

139
    @abc.abstractmethod
140
141
142
143
144
    def normalize_tag(self, tag):
        """Return a normalized version of the tag.

        Raises TypeError or ValueError if the tag is not acceptable.
        """
145
146
147
148
149
150
        pass

    def __call__(self, *tag):
        """
        A convenience function.

151
        This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).
152
153
154
        """
        # Catch a likely and difficult to find mistake.
        if tag and isinstance(tag[0], tuple):
155
156
            raise ValueError('Use site_family(1, 2) instead of '
                             'site_family((1, 2))!')
157
158
159
        return Site(self, tag)


160
161
class SimpleSiteFamily(SiteFamily):
    """A site family used as an example and for testing.
162

163
    A family of sites tagged by any python objects where object satisfied
164
165
    condition ``object == eval(repr(object))``.

166
    It exists to provide a basic site family that can be used for testing the
167
168
169
170
    builder module without other dependencies.  It can be also used to tag
    sites with non-numeric objects like strings should this every be useful.

    Due to its low storage efficiency for numbers it is not recommended to use
171
    `SimpleSiteFamily` when `kwant.lattice.Monatomic` would also work.
172
    """
173
174

    def __init__(self, name=None):
175
        canonical_repr = '{0}({1})'.format(self.__class__, repr(name))
Anton Akhmerov's avatar
Anton Akhmerov committed
176
        super().__init__(canonical_repr, name)
177

178
179
180
181
182
183
184
185
186
    def normalize_tag(self, tag):
        tag = tuple(tag)
        try:
            if eval(repr(tag)) != tag:
                raise RuntimeError()
        except:
            raise TypeError('It must be possible to recreate the tag from '
                            'its representation.')
        return tag
187

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220

def validate_hopping(hopping):
    """Verify that the argument is a valid hopping."""

    # This check is essential to maintaining the requirement that hoppings must
    # be tuples.  Without it, list "hoppings" would work in some cases
    # (e.g. with Builder.__contains__).
    if not isinstance(hopping, tuple):
        raise TypeError("{0} object is not a valid key.".format(
            type(hopping).__name__))

    # The following check is not strictly necessary (there would be an error
    # anyway), but the error message would be confusing.
    try:
        a, b = hopping
    except:
        raise IndexError("Only length 2 tuples (=hoppings) are valid keys.")

    # This check is essential for Builder.__contains__ - without it a builder
    # would simply "not contain" such invalid hoppings.  In other cases it
    # provides a nicer error message.
    for site in hopping:
        if not isinstance(site, Site):
            raise TypeError("Hopping elements must be Site objects, not {0}."
                            .format(type(site).__name__))

    # This again is an essential check.  Without it, builders would accept loop
    # hoppings.
    if a == b:
        raise ValueError("A hopping connects the following site to itself:\n"
                         "{0}".format(a))


Christoph Groth's avatar
Christoph Groth committed
221

222
223
################ Symmetries

224
class Symmetry(metaclass=abc.ABCMeta):
225
226
    """Abstract base class for spatial symmetries.

227
228
    Many physical systems possess a discrete spatial symmetry, which results in
    special properties of these systems.  This class is the standard way to
229
    describe discrete spatial symmetries in Kwant.  An instance of this class
230
231
232
    can be passed to a `Builder` instance at its creation.  The most important
    kind of symmetry is translational symmetry, used to define scattering
    leads.
233
234
235

    Each symmetry has a fundamental domain -- a set of sites and hoppings,
    generating all the possible sites and hoppings upon action of symmetry
236
237
238
239
240
241
242
243
244
    group elements.  A class derived from `Symmetry` has to implement mapping
    of any site or hopping into the fundamental domain, applying a symmetry
    group element to a site or a hopping, and a method `which` to determine the
    group element bringing some site from the fundamental domain to the
    requested one.  Additionally, it has to have a property `num_directions`
    returning the number of independent symmetry group generators (number of
    elementary periods for translational symmetry).

    A ``ValueError`` must be raised by the symmetry class whenever a symmetry
245
    is used together with sites whose site family is not compatible with it.  A
246
247
    typical example of this is when the vector defining a translational
    symmetry is not a lattice vector.
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
    """

    @abc.abstractproperty
    def num_directions(self):
        """Number of elementary periods of the symmetry."""
        pass

    @abc.abstractmethod
    def which(self, site):
        """Calculate the domain of the site.

        Return the group element whose action on a certain site from the
        fundamental domain will result in the given `site`.
        """
        pass

    @abc.abstractmethod
    def act(self, element, a, b=None):
        """Act with a symmetry group element on a site or hopping."""
        pass

    def to_fd(self, a, b=None):
        """Map a site or hopping to the fundamental domain.

        If `b` is None, return a site equivalent to `a` within the fundamental
        domain.  Otherwise, return a hopping equivalent to `(a, b)` but where
        the first element belongs to the fundamental domain.

        This default implementation works but may be not efficient.
277

278
        """
279
        return self.act(-self.which(a), a, b)
280
281
282
283
284
285
286
287
288
289
290
291

    def in_fd(self, site):
        """Tell whether `site` lies within the fundamental domain."""
        for d in self.which(site):
            if d != 0:
                return False
        return True


class NoSymmetry(Symmetry):
    """A symmetry with a trivial symmetry group."""

Daniel Jaschke's avatar
Daniel Jaschke committed
292
293
294
295
296
297
    def __eq__(self, other):
        return isinstance(other, NoSymmetry)

    def __ne__(self, other):
        return not self.__eq__(other)

298
299
300
301
302
303
304
    def __repr__(self):
        return 'NoSymmetry()'

    @property
    def num_directions(self):
        return 0

305
306
    _empty_array = ta.array((), int)

307
    def which(self, site):
308
        return self._empty_array
309
310
311
312

    def act(self, element, a, b=None):
        if element:
            raise ValueError('`element` must be empty for NoSymmetry.')
313
        return a if b is None else (a, b)
314
315
316
317
318
319
320

    def to_fd(self, a, b=None):
        return a if b is None else (a, b)

    def in_fd(self, site):
        return True

Christoph Groth's avatar
Christoph Groth committed
321

322
323
################ Hopping kinds

324
class HoppingKind:
325
326
    """A pattern for matching hoppings.

327
328
    A hopping ``(a, b)`` matches precisely when the site family of ``a`` equals
    `family_a` and that of ``b`` equals `family_b` and ``(a.tag - b.tag)`` is
329
    equal to `delta`.  In other words, the matching hoppings have the form:
330
    ``(family_a(x + delta), family_b(x))``
331
332
333
334
335

    Parameters
    ----------
    delta : Sequence of integers
        The sequence is interpreted as a vector with integer elements.
336
337
338
    family_a : `~kwant.builder.SiteFamily`
    family_b : `~kwant.builder.SiteFamily` or ``None`` (default)
        The default value means: use the same family as `family_a`.
339
340
341

    Notes
    -----
Christoph Groth's avatar
Christoph Groth committed
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
    A ``HoppingKind`` is a callable object: When called with a
    `~kwant.builder.Builder` as sole argument, an instance of this class will
    return an iterator over all possible matching hoppings whose sites are
    already present in the system.  The hoppings do *not* have to be already
    present in the system.  For example::

        kind = kwant.builder.HoppingKind((1, 0), lat)
        sys[kind(sys)] = 1

    Because a `~kwant.builder.Builder` can be indexed with functions or
    iterables of functions, ``HoppingKind`` instances (or any non-tuple
    iterables of them, e.g. a list) can be used directly as "wildcards" when
    setting or deleting hoppings::

        kinds = [kwant.builder.HoppingKind(v, lat) for v in [(1, 0), (0, 1)]]
        sys[kinds] = 1
358
    """
359
    __slots__ = ('delta', 'family_a', 'family_b')
360

361
    def __init__(self, delta, family_a, family_b=None):
362
        self.delta = ta.array(delta, int)
363
        ensure_isinstance(family_a, SiteFamily)
364
        self.family_a = family_a
365
366
367
368
369
        if family_b is None:
            self.family_b = family_a
        else:
            ensure_isinstance(family_b, SiteFamily)
            self.family_b = family_b
370

371
    def __call__(self, builder):
372
        delta = self.delta
373
374
        family_a = self.family_a
        family_b = self.family_b
375
376
377
378
        H = builder.H
        symtofd = builder.symmetry.to_fd

        for a in H:
379
            if a.family != family_a:
380
                continue
381
            b = Site(family_b, a.tag - delta, True)
382
383
384
385
            if symtofd(b) in H:
                yield a, b

    def __repr__(self):
Christoph Groth's avatar
Christoph Groth committed
386
        return '{0}({1}, {2}{3})'.format(
387
            self.__class__.__name__, repr(tuple(self.delta)),
388
389
            repr(self.family_a),
            ', ' + repr(self.family_b) if self.family_a != self.family_b else '')
390

391
392
    def __str__(self):
        return '{0}({1}, {2}{3})'.format(
393
394
            self.__class__.__name__, tuple(self.delta),
            self.family_a,
395
396
397
            ', ' + str(self.family_b) if self.family_a != self.family_b else '')


Christoph Groth's avatar
Christoph Groth committed
398

399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
################ Support for Hermitian conjugation

def herm_conj(value):
    """
    Calculate the hermitian conjugate of a python object.

    If the object is neither a complex number nor a matrix, the original value
    is returned.  In the context of this module, this is the correct behavior
    for functions.
    """
    if hasattr(value, 'conjugate'):
        value = value.conjugate()
        if hasattr(value, 'transpose'):
            value = value.transpose()
    return value


416
class HermConjOfFunc:
417
418
    """Proxy returning the hermitian conjugate of the original result."""
    __slots__ = ('function')
Anton Akhmerov's avatar
Anton Akhmerov committed
419

420
421
422
    def __init__(self, function):
        self.function = function

423
424
    def __call__(self, i, j, *args):
        return herm_conj(self.function(j, i, *args))
425

Christoph Groth's avatar
Christoph Groth committed
426

427
428
################ Leads

429
class Lead(metaclass=abc.ABCMeta):
430
431
    """Abstract base class for leads that can be attached to a `Builder`.

432
433
434
435
    To attach a lead to a builder, append it to the builder's `~Builder.leads`
    instance variable.  See the documentation of `kwant.builder` for the
    concrete classes of leads derived from this one.

436
437
    Attributes
    ----------
438
    interface : sequence of sites
439

440
441
442
    """

    @abc.abstractmethod
443
    def finalized(self):
444
445
446
447
448
449
450
451
        """Return a finalized version of the lead.

        Returns
        -------
        finalized_lead

        Notes
        -----
Christoph Groth's avatar
Christoph Groth committed
452
453
454
        The finalized lead must be an object that can be used as a lead in a
        `kwant.system.FiniteSystem`.  It could be an instance of
        `kwant.system.InfiniteSystem` for example.
455

Christoph Groth's avatar
Christoph Groth committed
456
457
        The order of sites for the finalized lead must be the one specified in
        `interface`.
458
459
460
461
462
463
464
465
466
467
468
469
        """
        pass


class BuilderLead(Lead):
    """A lead made from a `Builder` with a spatial symmetry.

    Parameters
    ----------
    builder : `Builder`
        The tight-binding system of a lead. It has to possess appropriate
        symmetry, and it may not contain hoppings between further than
Christoph Groth's avatar
Christoph Groth committed
470
        neighboring images of the fundamental domain.
471
    interface : sequence of `Site` instances
472
473
474
475
476
        Sequence of sites in the scattering region to which the lead is
        attached.

    Notes
    -----
477
478
479
480
    The hopping from the scattering region to the lead is assumed to be equal to
    the hopping from a lead unit cell to the next one in the direction of the
    symmetry vector (i.e. the lead is 'leaving' the system and starts with a
    hopping).
481

482
    The given order of interface sites is preserved throughout finalization.
483
484
485
486
487

    Every system has an attribute `leads`, which stores a list of
    `BuilderLead` objects with all the information about the leads that are
    attached.
    """
488
    def __init__(self, builder, interface):
489
        self.builder = builder
490
        self.interface = tuple(interface)
491
492
493
494
495

    def finalized(self):
        """Return a `kwant.system.InfiniteSystem` corresponding to the
        compressed lead.

496
        The order of interface sites is kept during finalization.
497
        """
498
        return self.builder._finalized_infinite(self.interface)
499
500


501
class SelfEnergyLead(Lead):
502
503
504
505
    """A general lead defined by its self energy.

    Parameters
    ----------
506
    selfenergy_func : function
507
        Function which returns the self energy matrix for the interface sites
508
        given the energy and optionally a sequence of extra arguments.
509
    interface : sequence of `Site` instances
510
    """
511
    def __init__(self, selfenergy_func, interface):
512
        self.selfenergy_func = selfenergy_func
513
        self.interface = tuple(interface)
514
515
516
517
518

    def finalized(self):
        """Trivial finalization: the object is returned itself."""
        return self

519
    def selfenergy(self, energy, args=()):
520
        return self.selfenergy_func(energy, args)
521
522
523
524
525
526
527
528


class ModesLead(Lead):
    """A general lead defined by its modes wave functions.

    Parameters
    ----------
    modes_func : function
529
530
        Function which returns the modes of the lead as a tuple of
        `~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes`
531
        given the energy and optionally a sequence of extra arguments.
532
533
534
    interface :
        sequence of `Site` instances

535
536
    """
    def __init__(self, modes_func, interface):
537
        self.modes_func = modes_func
538
539
540
541
542
543
544
        self.interface = tuple(interface)

    def finalized(self):
        """Trivial finalization: the object is returned itself."""
        return self

    def modes(self, energy, args=()):
545
        return self.modes_func(energy, args)
546

547
    def selfenergy(self, energy, args=()):
548
        stabilized = self.modes(energy, args)[1]
549
        return stabilized.selfenergy()
Christoph Groth's avatar
Christoph Groth committed
550

551

Christoph Groth's avatar
Christoph Groth committed
552

553
554
################ Builder class

555
# A marker, meaning for hopping (i, j): this value is given by the Hermitian
556
# conjugate the value of the hopping (j, i).  Used by Builder and System.
557
class Other:
Christoph Groth's avatar
Christoph Groth committed
558
    pass
559

560

561
562
563
564
def edges(seq):
    # izip, when given the same iterator twice, turns a sequence into a
    # sequence of pairs.
    seq_iter = iter(seq)
Joseph Weston's avatar
Joseph Weston committed
565
    result = zip(seq_iter, seq_iter)
566
567
    next(result)                # Skip the special loop edge.
    return result
568

569

570
class Builder:
571
572
    """A tight binding system defined on a graph.

573
    This is one of the central types in Kwant.  It is used to construct tight
574
575
576
577
    binding systems in a flexible way.

    The nodes of the graph are `Site` instances.  The edges, i.e. the hoppings,
    are pairs (2-tuples) of sites.  Each node and each edge has a value
Christoph Groth's avatar
Christoph Groth committed
578
    associated with it.  The values associated with nodes are interpreted as
579
580
    on-site Hamiltonians, the ones associated with edges as hopping integrals.

Christoph Groth's avatar
Christoph Groth committed
581
    To make the graph accessible in a way that is natural within the Python
582
    language it is exposed as a *mapping* (much like a built-in Python
Christoph Groth's avatar
Christoph Groth committed
583
    dictionary).  Keys are sites or hoppings.  Values are 2d arrays
584
    (e.g. NumPy or Tinyarray) or numbers (interpreted as 1 by 1 matrices).
585
586
587
588
589
590
591
592

    Parameters
    ----------
    symmetry : `Symmetry` or `None`
        The symmetry of the system.

    Notes
    -----
Christoph Groth's avatar
Christoph Groth committed
593
594
595
596
597
598
    Values can be also functions that receive the site or the hopping (passed
    to the function as two sites) and possibly additional arguments and are
    expected to return a valid value.  This allows to define systems quickly,
    to modify them without reconstructing, and to save memory for many-orbital
    models.

599
600
    In addition to simple keys (single sites and hoppings) more powerful keys
    are possible as well that allow to manipulate multiple sites/hoppings in a
Christoph Groth's avatar
Christoph Groth committed
601
602
603
    single operation.  Such keys are internally expanded into a sequence of
    simple keys by using the method `Builder.expand`.  For example,
    ``sys[general_key] = value`` is equivalent to ::
604

Christoph Groth's avatar
Christoph Groth committed
605
        for simple_key in sys.expand(general_key):
606
            sys[simple_key] = value
Christoph Groth's avatar
Christoph Groth committed
607

608
609
610
611
    Builder instances automatically ensure that every hopping is Hermitian, so
    that if ``builder[a, b]`` has been set, there is no need to set
    ``builder[b, a]``.

Christoph Groth's avatar
Christoph Groth committed
612
613
614
    Builder instances can be made to automatically respect a `Symmetry` that is
    passed to them during creation.  The behavior of builders with a symmetry
    is slightly more sophisticated.  First of all, it is implicitly assumed
615
    throughout Kwant that **every** function assigned as a value to a builder
Christoph Groth's avatar
Christoph Groth committed
616
617
618
    with a symmetry possesses the same symmetry.  Secondly, all keys are mapped
    to the fundamental domain of the symmetry before storing them.  This may
    produce confusing results when neighbors of a site are queried.
619

Christoph Groth's avatar
Christoph Groth committed
620
621
622
    The method `attach_lead` *works* only if the sites affected by them have
    tags which are sequences of integers.  It *makes sense* only when these
    sites live on a regular lattice, like the ones provided by `kwant.lattice`.
623

624
625
626
    Attaching a lead manually (without the use of `~Builder.attach_lead`)
    amounts to creating a `Lead` object and appending it to this list.

Christoph Groth's avatar
Christoph Groth committed
627
628
629
630
    ``builder0 += builder1`` adds all the sites, hoppings, and leads of
    ``builder1`` to ``builder0``.  Sites and hoppings present in both systems
    are overwritten by those in ``builder1``.  The leads of ``builder1`` are
    appended to the leads of the system being extended.
631

632
633
634
635
636
637
    .. warning::

        If functions are used to set values in a builder with a symmetry, then
        they must satisfy the same symmetry.  There is (currently) no check and
        wrong results will be the consequence of a misbehaving function.

638
639
640
641
642
    Attributes
    ----------
    leads : list of `Lead` instances
        The leads that are attached to the system.

643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
    Examples
    --------
    Define a site.

    >>> builder[site] = value

    Print the value of a site.

    >>> print builder[site]

    Define a hopping.

    >>> builder[site1, site2] = value

    Delete a site.

    >>> del builder[site3]

661
662
663
664
665
    Detach the last lead.  (This does not remove the sites that were added to
    the scattering region by `~Builder.attach_lead`.)

    >>> del builder.leads[-1]

666
    """
667

668
669
670
    def __init__(self, symmetry=None):
        if symmetry is None:
            symmetry = NoSymmetry()
671
672
        else:
            ensure_isinstance(symmetry, Symmetry)
673
674
        self.symmetry = symmetry
        self.leads = []
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
        self.H = {}

    #### Note on H ####
    #
    # This dictionary stores a directed graph optimized for efficient querying
    # and modification.  The nodes are instances of `Site`.
    #
    # Each edge, specified by a ``(tail, head)`` pair of nodes, holds an object
    # as a value.  Likewise, each tail which occurs in the graph also holds a
    # value.  (Nodes which only occur as heads are not required to have
    # values.)
    #
    # For a given `tail` site, H[tail] is a list alternately storing
    # heads and values.  (The heads occupy even locations followed by the
    # values at odd locations.)  Each pair of entries thus describes a single
    # directed edge of the graph.
    #
    # The first pair of entries in each list is special: it always
    # corresponds to a loop edge.  (The head is equal to the tail.)  This
    # special edge has two purposes: It is used to store the value
    # associated with the tail node itself, and it is necessary for the
    # method getkey_tail which helps to conserve memory by storing equal
    # node label only once.

    def _get_edge(self, tail, head):
        for h, value in edges(self.H[tail]):
            if h == head:
                return value
703
704
705
706
707
708
709
710

        # (tail, head) is not present in the system, but tail is.
        if head in self.H:
            raise KeyError((tail, head))
        else:
            # If already head is missing, we only report this.  This way the
            # behavior is symmetric with regard to tail and head.
            raise KeyError(head)
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725

    def _set_edge(self, tail, head, value):
        hvhv = self.H[tail]
        heads = hvhv[2::2]
        if head in heads:
            i = 2 + 2 * heads.index(head)
            hvhv[i] = head
            hvhv[i + 1] = value
        else:
            hvhv.append(head)
            hvhv.append(value)

    def _del_edge(self, tail, head):
        hvhv = self.H[tail]
        heads = hvhv[2::2]
726
727
728
729
730
731
732
733
734
735
736
737

        try:
            i = 2 + 2 * heads.index(head)
        except ValueError:
            # (tail, head) is not present in the system, but tail is.
            if head in self.H:
                raise KeyError((tail, head))
            else:
                # If already head is missing, we only report this.  This way the
                # behavior is symmetric with regard to tail and head.
                raise KeyError(head)

738
739
740
        del hvhv[i : i + 2]

    def _out_neighbors(self, tail):
741
        hvhv = self.H[tail]
742
743
744
        return islice(hvhv, 2, None, 2)

    def _out_degree(self, tail):
745
        hvhv = self.H[tail]
746
        return len(hvhv) // 2 - 1
747

748
    # TODO: write a test for this method.
749
750
751
752
753
754
755
756
757
758
759
760
    def reversed(self):
        """Return a shallow copy of the builder with the symmetry reversed.

        This method can be used to attach the same infinite system as lead from
        two opposite sides.  It requires a builder to which an infinite
        symmetry is associated.
        """
        result = object.__new__(Builder)
        result.symmetry = self.symmetry.reversed()
        if self.leads:
            raise ValueError('System to be reversed may not have leads.')
        result.leads = []
761
        result.H = self.H
762
763
        return result

Joseph Weston's avatar
Joseph Weston committed
764
    def __bool__(self):
765
        return bool(self.H)
766

767
768
    def expand(self, key):
        """
Christoph Groth's avatar
Christoph Groth committed
769
        Expand a general key into an iterator over simple keys.
770

Christoph Groth's avatar
Christoph Groth committed
771
772
773
774
        Parameters
        ----------
        key : builder key (see notes)
            The key to be expanded
775

Christoph Groth's avatar
Christoph Groth committed
776
777
        Notes
        -----
778
779
780
781
782
783
784
        Keys are (recursively):
            * Simple keys: sites or 2-tuples of sites (=hoppings).
            * Any (non-tuple) iterable of keys, e.g. a list or a generator
              expression.
            * Any function that returns a key when passed a builder as sole
              argument, e.g. a `HoppingKind` instance or the function returned
              by `~kwant.lattice.Polyatomic.shape`.
Christoph Groth's avatar
Christoph Groth committed
785
786
787
788
789

        This method is internally used to expand the keys when getting or
        deleting items of a builder (i.e. ``sys[key] = value`` or ``del
        sys[key]``).

790
791
792
793
794
        """
        itr = iter((key,))
        iter_stack = [None]
        while iter_stack:
            for key in itr:
795
                while callable(key):
Christoph Groth's avatar
Christoph Groth committed
796
                    key = key(self)
797
798
                if isinstance(key, tuple):
                    # Site instances are also tuples.
799
                    yield key
800
                else:
801
                    iter_stack.append(itr)
802
803
804
805
806
                    try:
                        itr = iter(key)
                    except TypeError:
                        raise TypeError("{0} object is not a valid key."
                                        .format(type(key).__name__))
807
808
809
                    break
            else:
                itr = iter_stack.pop()
810

811
812
813
814
    def __getitem__(self, key):
        """Get the value of a single site or hopping."""
        if isinstance(key, Site):
            site = self.symmetry.to_fd(key)
815
            return self.H[site][1]
816
817

        sym = self.symmetry
818
819
820
        validate_hopping(key)
        a, b = sym.to_fd(*key)
        value = self._get_edge(a, b)
821
        if value is Other:
822
823
824
            if not sym.in_fd(b):
                b, a = sym.to_fd(b, a)
                assert not sym.in_fd(a)
825
            value = self._get_edge(b, a)
826
            if callable(value):
827
828
829
830
831
832
833
834
                assert not isinstance(value, HermConjOfFunc)
                value = HermConjOfFunc(value)
            else:
                value = herm_conj(value)
        return value

    def __contains__(self, key):
        """Tell whether the system contains a site or hopping."""
835
        if isinstance(key, Site):
836
837
838
839
840
841
842
            key = self.symmetry.to_fd(key)
            return key in self.H

        validate_hopping(key)
        a, b = self.symmetry.to_fd(*key)
        hvhv = self.H.get(a, ())
        return b in islice(hvhv, 2, None, 2)
843

844
    def _set_site(self, site, value):
845
        """Set a single site."""
846
847
        if not isinstance(site, Site):
            raise TypeError('Expecting a site, got {0} instead.'.format(type(site).__name__))
848
        site = self.symmetry.to_fd(site)
849
850
851
852
853
        hvhv = self.H.setdefault(site, [])
        if hvhv:
            hvhv[1] = value
        else:
            hvhv[:] = [site, value]
854

855
    def _set_hopping(self, hopping, value):
856
857
        """Set a single hopping."""
        sym = self.symmetry
858
859
860
861
862
863
864
865
866
867
868
869
870
        validate_hopping(hopping)
        a, b = sym.to_fd(*hopping)

        if sym.in_fd(b):
            # Make sure that we do not waste space by storing multiple instances
            # of identical sites.
            a2 = a = self.H[a][0]
            b2 = b = self.H[b][0]
        else:
            b2, a2 = sym.to_fd(b, a)
            assert not sym.in_fd(a2)
            if b2 not in self.H:
                raise KeyError(b)
871

872
873
874
875
876
877
878
879
880
        # It's important that we have verified that b2 belongs to the system.
        # Otherwise, we risk ending up with a half-added hopping.
        if isinstance(value, HermConjOfFunc):
            # Avoid nested HermConjOfFunc instances.
            self._set_edge(a, b, Other)            # May raise KeyError(a).
            self._set_edge(b2, a2, value.function) # Must succeed.
        else:
            self._set_edge(a, b, value)            # May raise KeyError(a).
            self._set_edge(b2, a2, Other)          # Must succeed.
881
882

    def __setitem__(self, key, value):
883
884
885
886
        """Set a single site/hopping or a bunch of them."""
        func = None
        for sh in self.expand(key):
            if func is None:
887
888
                func = (self._set_site if isinstance(sh, Site)
                        else self._set_hopping)
889
            func(sh, value)
890

891
    def _del_site(self, site):
892
        """Delete a single site and all associated hoppings."""
893
894
895
896
        if not isinstance(site, Site):
            raise TypeError('Expecting a site, got {0} instead.'.format(
                type(site).__name__))

897
        tfd = self.symmetry.to_fd
898
        site = tfd(site)
899
900
901
902
903
904
905
906
907
908
909

        out_neighbors = self._out_neighbors(site)

        for neighbor in out_neighbors:
            if neighbor in self.H:
                self._del_edge(neighbor, site)
            else:
                assert not self.symmetry.in_fd(neighbor)
                self._del_edge(*tfd(neighbor, site))

        del self.H[site]
910

911
    def _del_hopping(self, hopping):
912
913
        """Delete a single hopping."""
        sym = self.symmetry
914
915
916
        validate_hopping(hopping)
        a, b = sym.to_fd(*hopping)
        self._del_edge(a, b)
917

918
919
920
921
922
923
        if sym.in_fd(b):
            self._del_edge(b, a)
        else:
            b, a = sym.to_fd(b, a)
            assert not sym.in_fd(a)
            self._del_edge(b, a)
924
925

    def __delitem__(self, key):
926
927
928
929
        """Delete a single site/hopping or bunch of them."""
        func = None
        for sh in self.expand(key):
            if func is None:
930
931
                func = (self._del_site if isinstance(sh, Site)
                        else self._del_hopping)
932
            func(sh)
933
934
935

    def eradicate_dangling(self):
        """Keep deleting dangling sites until none are left."""
936
937
        sites = list(site for site in self.H
                     if self._out_degree(site) < 2)
938
        for site in sites:
Anton Akhmerov's avatar
Anton Akhmerov committed
939
940
            if site not in self.H:
                continue
941
            while site:
942
943
944
945
946
947
948
                neighbors = tuple(self._out_neighbors(site))
                if neighbors:
                    assert len(neighbors) == 1
                    neighbor = neighbors[0]
                    self._del_edge(neighbor, site)
                    if self._out_degree(neighbor) > 1:
                        neighbor = False
949
                else:
950
                    neighbor = False
951
                del self.H[site]
952
                site = neighbor
953
954
955

    def __iter__(self):
        """Return an iterator over all sites and hoppings."""
956
        return chain(self.H, self.hoppings())
957
958

    def sites(self):
Anton Akhmerov's avatar
Anton Akhmerov committed
959
960
961
962
963
964
        """Return a read-only set over all sites.

        The sites that are returned belong to the fundamental domain of the
        `Builder` symmetry, and are not necessarily the ones that were set
        initially (but always the equivalent ones).
        """
965
        try:
Joseph Weston's avatar
Joseph Weston committed
966
            return self.H.keys()
967
968
        except AttributeError:
            return frozenset(self.H)
969
970
971

    def site_value_pairs(self):
        """Return an iterator over all (site, value) pairs."""
Joseph Weston's avatar
Joseph Weston committed
972
        for site, hvhv in self.H.items():
973
            yield site, hvhv[1]
974
975

    def hoppings(self):
Anton Akhmerov's avatar
Anton Akhmerov committed
976
977
978
979
980
981
        """Return an iterator over all Builder hoppings.

        The hoppings that are returned belong to the fundamental domain of the
        `Builder` symmetry, and are not necessarily the ones that were set
        initially (but always the equivalent ones).
        """
Joseph Weston's avatar
Joseph Weston committed
982
        for tail, hvhv in self.H.items():
983
            for head, value in edges(hvhv):
984
                if value is Other:
Anton Akhmerov's avatar
Anton Akhmerov committed
985
                    continue
986
                yield tail, head
987
988
989

    def hopping_value_pairs(self):
        """Return an iterator over all (hopping, value) pairs."""
Joseph Weston's avatar
Joseph Weston committed
990
        for tail, hvhv in self.H.items():
991
            for head, value in edges(hvhv):
992
                if value is Other:
Anton Akhmerov's avatar
Anton Akhmerov committed
993
                    continue
994
                yield (tail, head), value
995
996
997

    def dangling(self):
        """Return an iterator over all dangling sites."""
998
999
        for site in self.H:
            if self._out_degree(site) < 2:
1000
                yield site
For faster browsing, not all history is shown. View entire blame