test_builder.py 58.7 KB
Newer Older
1
# Copyright 2011-2019 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
import warnings
Anton Akhmerov's avatar
Anton Akhmerov committed
10
import pickle
11
import itertools as it
12
13
14
15
16
import functools as ft
from random import Random

import numpy as np
import tinyarray as ta
17
import pytest
Christoph Groth's avatar
Christoph Groth committed
18
from pytest import raises, warns
19
from numpy.testing import assert_almost_equal
20

21
import kwant
22
from kwant import builder, system
23
from kwant._common import ensure_rng, KwantDeprecationWarning
24
25


26
27
28
def test_bad_keys():

    def setitem(key):
29
        syst[key] = None
30

31
    fam = builder.SimpleSiteFamily(norbs=1)
32
    syst = builder.Builder()
33
34
35

    failures = [
        # Invalid single keys
36
        ([syst.__contains__, syst.__getitem__, setitem, syst.__delitem__],
37
38
39
40
41
42
43
44
45
46
         [(TypeError, [123,
                       (0, 1),
                       (fam(0), 123),
                       (123, (fam(0)))]),
          (IndexError, [(fam(0),),
                        (fam(0), fam(1), fam(2))]),
          (ValueError, [(fam(0), fam(0)),
                        (fam(2), fam(2))])]),

        # Hoppings that contain sites that do not belong to the system
47
        ([syst.__getitem__, setitem, syst.__delitem__],
48
49
50
51
         [(KeyError, [(fam(0), fam(3)), (fam(2), fam(1)),
                      (fam(2), fam(3))])]),

        # Sequences containing a bad key.
52
        ([setitem, syst.__delitem__],
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
         [(TypeError, [[fam(0), fam(1), 123],
                       [fam(0), (fam(1),)],
                       [fam(0), (fam(1), fam(2))],
                       [(fam(0), fam(1)), (0, 1)],
                       [(fam(0), fam(1)), (fam(0), 123)],
                       [(fam(0), fam(1)), (123, fam(0))],
                       [(fam(0), fam(1)), fam(2)]]),
          (IndexError, [[(fam(0), fam(1)), (fam(2),)]]),
          (ValueError, [[(fam(0), fam(1)), (fam(2), fam(2))],
                        [(fam(0), fam(0)), (fam(1), fam(0))]]),
          (KeyError, [[(fam(0), fam(1)), (fam(0), fam(3))],
                      [(fam(0), fam(1)), (fam(2), fam(1))],
                      [(fam(1), fam(2)), (fam(0), fam(1))]])]),

        # Sites that do not belong to the system, also as part of a
        # sequence
69
        ([syst.__delitem__],
70
71
72
73
74
         [(KeyError, [fam(123),
                      [fam(0), fam(123)],
                      [fam(123), fam(1)]])]),

        # Various things that are not sites present in the system.
75
        ([syst.degree, lambda site: list(syst.neighbors(site))],
76
77
78
79
80
81
82
83
84
85
86
87
88
         [(TypeError, [123,
                       [0, 1, 2],
                       (0, 1),
                       (fam(0), fam(1)),
                       [fam(0), fam(1)],
                       [fam(1), fam(2)],
                       [fam(3), fam(0)]]),
          (KeyError, [fam(123)])])]

    for funcs, errors in failures:
        for error, keys in errors:
            for key in keys:
                for func in funcs:
89
90
                    syst[[fam(0), fam(1)]] = None
                    syst[fam(0), fam(1)] = None
91
                    try:
Anton Akhmerov's avatar
Anton Akhmerov committed
92
                        raises(error, func, key)
93
                    except AssertionError:
Joseph Weston's avatar
Joseph Weston committed
94
                        print(func, error, key)
95
96
97
                        raise


98
def test_site_families():
99
    syst = builder.Builder()
100
101
102
    fam = builder.SimpleSiteFamily(norbs=1)
    ofam = builder.SimpleSiteFamily(norbs=1)
    yafam = builder.SimpleSiteFamily('another_name', norbs=1)
103

104
    syst[fam(0)] = 7
105
    assert syst[fam(0)] == 7
106

107
    assert len(set([fam, ofam, fam('a'), ofam('a'), yafam])) == 3
108
    syst[fam(1)] = 123
109
110
    assert syst[fam(1)] == 123
    assert syst[ofam(1)] == 123
Anton Akhmerov's avatar
Anton Akhmerov committed
111
    raises(KeyError, syst.__getitem__, yafam(1))
112

113
    # test site families compare equal/not-equal
114
    assert fam == ofam
Anton Akhmerov's avatar
Anton Akhmerov committed
115
116
117
    assert fam != yafam
    assert fam != None
    assert fam != 'a'
118

119
120
121
    # test site families sorting
    fam1 = builder.SimpleSiteFamily(norbs=1)
    fam2 = builder.SimpleSiteFamily(norbs=2)
Anton Akhmerov's avatar
Anton Akhmerov committed
122
    assert fam1 < fam2  # string '1' is lexicographically less than '2'
123
124


125
126
127
class VerySimpleSymmetry(builder.Symmetry):
    def __init__(self, period):
        self.period = period
128

129
130
131
    @property
    def num_directions(self):
        return 1
132

133
    def has_subgroup(self, other):
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
        if isinstance(other, builder.NoSymmetry):
            return True
        elif isinstance(other, VerySimpleSymmetry):
            return not other.period % self.period
        else:
            return False

    def subgroup(self, *generators):
        generators = ta.array(generators)
        assert generators.shape == (1, 1)
        if generators.dtype != int:
            raise ValueError('Generators must be sequences of integers.')
        g = generators[0, 0]
        return VerySimpleSymmetry(g * self.period)

149
150
151
152
    def which(self, site):
        return ta.array((site.tag[0] // self.period,), int)

    def act(self, element, a, b=None):
153
        shifted = lambda site, delta: site.family(*ta.add(site.tag, delta))
154
155
        delta = (self.period * element[0],) + (len(a.tag) - 1) * (0,)
        if b is None:
Anton Akhmerov's avatar
Anton Akhmerov committed
156
            return shifted(a, delta)
157
        else:
Anton Akhmerov's avatar
Anton Akhmerov committed
158
            return shifted(a, delta), shifted(b, delta)
159
160
161
162
163


# The hoppings have to form a ring.  Some other implicit assumptions are also
# made.
def check_construction_and_indexing(sites, sites_fd, hoppings, hoppings_fd,
164
                                    unknown_hoppings, sym=None):
165
    fam = builder.SimpleSiteFamily(norbs=1)
166
    syst = builder.Builder(sym)
167
    t, V = 1.0j, 0.0
168
    syst[sites] = V
169
    for site in sites:
170
171
        syst[site] = V
    syst[hoppings] = t
172
    for hopping in hoppings:
173
        syst[hopping] = t
174

175
    for hopping in unknown_hoppings:
Anton Akhmerov's avatar
Anton Akhmerov committed
176
        raises(KeyError, syst.__setitem__, hopping, t)
177

178
179
180
    assert (fam(5), fam(123)) not in syst
    assert (sites[0], fam(5, 123)) not in syst
    assert (fam(7, 8), sites[0]) not in syst
181
    for site in sites:
182
        assert site in syst
Anton Akhmerov's avatar
Anton Akhmerov committed
183
        assert syst[site] == V
184
185
    for hop in hoppings:
        rev_hop = hop[1], hop[0]
186
187
        assert hop in syst
        assert rev_hop in syst
188
189
        assert syst[hop] == t
        assert syst[rev_hop] == t.conjugate()
190

191
192
193
    assert syst.degree(sites[0]) == 2
    assert (sorted(s for s in syst.neighbors(sites[0])) ==
            sorted([sites[1], sites[-1]]))
194

195
    del syst[hoppings]
196
    assert list(syst.hoppings()) == []
197
    syst[hoppings] = t
198

199
    del syst[sites[0]]
200
201
202
    assert sorted(tuple(s) for s in syst.sites()) == sorted(sites_fd[1:])
    assert (sorted((a, b) for a, b in syst.hoppings()) ==
            sorted(hoppings_fd[1:-1]))
203

204
205
206
207
208
209
210
    assert (sorted((tuple(site.tag), value)
                        for site, value in syst.site_value_pairs()) ==
            sorted((tuple(site.tag), syst[site]) for site in syst.sites()))
    assert (sorted((tuple(a.tag), tuple(b.tag), value)
                   for (a, b), value in syst.hopping_value_pairs()) ==
            sorted((tuple(a.tag), tuple(b.tag), syst[a, b])
                   for a, b in syst.hoppings()))
211

Anton Akhmerov's avatar
Anton Akhmerov committed
212

213
214
def test_construction_and_indexing():
    # Without symmetry
215
    fam = builder.SimpleSiteFamily(norbs=1)
216
217
218
219
    sites = [fam(0, 0), fam(0, 1), fam(1, 0)]
    hoppings = [(fam(0, 0), fam(0, 1)),
                (fam(0, 1), fam(1, 0)),
                (fam(1, 0), fam(0, 0))]
220
    unknown_hoppings = [(fam(0, 1), fam(7, 8)),
221
                        (fam(12, 14), fam(0, 1))]
222
    check_construction_and_indexing(sites, sites, hoppings, hoppings,
223
                                    unknown_hoppings)
224
225

    # With symmetry
226
227
228
229
230
231
232
233
234
235
    sites = [fam(0, 0), fam(1, 1), fam(2, 1), fam(4, 2)]
    sites_fd = [fam(0, 0), fam(1, 1), fam(0, 1), fam(0, 2)]
    hoppings = [(fam(0, 0), fam(1, 1)),
                (fam(1, 1), fam(2, 1)),
                (fam(2, 1), fam(4, 2)),
                (fam(4, 2), fam(0, 0))]
    hoppings_fd = [(fam(0, 0), fam(1, 1)),
                   (fam(1, 1), fam(2, 1)),
                   (fam(0, 1), fam(2, 2)),
                   (fam(0, 2), fam(-4, 0))]
236
    unknown_hoppings = [(fam(0, 0), fam(0, 3)), (fam(0, 4), fam(0, 0)),
237
238
                        (fam(0, 0), fam(2, 3)), (fam(2, 4), fam(0, 0)),
                        (fam(4, 2), fam(6, 3)), (fam(6, 4), fam(4, 2))]
239
240
    sym = VerySimpleSymmetry(2)
    check_construction_and_indexing(sites, sites_fd, hoppings, hoppings_fd,
241
                                    unknown_hoppings, sym)
242
243
244


def test_hermitian_conjugation():
245
    def f(i, j, arg):
246
        i, j = i.tag, j.tag
247
        if j[0] == i[0] + 1:
248
            return arg * ta.array([[1, 2j], [3 + 1j, 4j]])
249
250
251
        else:
            raise ValueError

252
    syst = builder.Builder()
253
    fam = builder.SimpleSiteFamily(norbs=1)
254
    syst[fam(0)] = syst[fam(1)] = ta.identity(2)
255

256
257
258
    syst[fam(0), fam(1)] = f
    assert syst[fam(0), fam(1)] is f
    assert isinstance(syst[fam(1), fam(0)], builder.HermConjOfFunc)
259
260
    assert (syst[fam(1), fam(0)](fam(1), fam(0), 2) ==
            syst[fam(0), fam(1)](fam(0), fam(1), 2).conjugate().transpose())
261
262
263
    syst[fam(0), fam(1)] = syst[fam(1), fam(0)]
    assert isinstance(syst[fam(0), fam(1)], builder.HermConjOfFunc)
    assert syst[fam(1), fam(0)] is f
264
265
266


def test_value_equality_and_identity():
267
    m = ta.array([[1, 2], [3j, 4j]])
268
    syst = builder.Builder()
269
    fam = builder.SimpleSiteFamily(norbs=1)
270

271
272
273
    syst[fam(0)] = m
    syst[fam(1)] = m
    assert syst[fam(1)] is m
274

275
    syst[fam(0), fam(1)] = m
276
    assert syst[fam(1), fam(0)] == m.transpose().conjugate()
277
    assert syst[fam(0), fam(1)] is m
278

279
    syst[fam(1), fam(0)] = m
280
    assert syst[fam(0), fam(1)] == m.transpose().conjugate()
281
    assert syst[fam(1), fam(0)] is m
282
283
284
285
286


def random_onsite_hamiltonian(rng):
    return 2 * rng.random() - 1

Anton Akhmerov's avatar
Anton Akhmerov committed
287

288
289
290
def random_hopping_integral(rng):
    return complex(2 * rng.random() - 1, 2 * rng.random() - 1)

Anton Akhmerov's avatar
Anton Akhmerov committed
291

292
def check_onsite(fsyst, sites, subset=False, check_values=True):
293
294
295
296
297
    vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))

    if vectorized:
        site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])

298
    freq = {}
299
300
    for node in range(fsyst.graph.num_nodes):
        site = fsyst.sites[node].tag
301
302
        freq[site] = freq.get(site, 0) + 1
        if check_values and site in sites:
303
304
305
306
307
308
309
310
311
312
313
314
            if vectorized:
                term = fsyst._onsite_term_by_site_id[node]
                value = fsyst._term_values[term]
                if callable(value):
                    assert value is sites[site]
                else:
                    (w, _), (off, _) = fsyst.subgraphs[fsyst.terms[term].subgraph]
                    node_off = node - site_offsets[w]
                    selector = np.searchsorted(off, node_off)
                    assert np.allclose(value[selector], sites[site])
            else:
                assert fsyst.onsites[node][0] is sites[site]
315
    if not subset:
316
        # Check that all sites of `fsyst` are in `sites`.
Joseph Weston's avatar
Joseph Weston committed
317
        for site in freq.keys():
318
            assert site in sites
319
    # Check that all sites of `sites` are in `fsyst`.
320
    for site in sites:
321
        assert freq[site] == 1
322

Anton Akhmerov's avatar
Anton Akhmerov committed
323

324
def check_hoppings(fsyst, hops):
325
326
327
328
329
    vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))

    if vectorized:
        site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])

330
    assert fsyst.graph.num_edges == 2 * len(hops)
331
    for edge_id, edge in enumerate(fsyst.graph):
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
        i, j = edge
        tail = fsyst.sites[i].tag
        head = fsyst.sites[j].tag

        if vectorized:
            term = fsyst._hopping_term_by_edge_id[edge_id]
            if term < 0:  # Hermitian conjugate
                assert (head, tail) in hops
            else:
                value = fsyst._term_values[term]
                assert (tail, head) in hops
                if callable(value):
                    assert value is hops[tail, head]
                else:
                    dtype = np.dtype([('f0', int), ('f1', int)])
                    subgraph = fsyst.terms[term].subgraph
                    (to_w, from_w), hoppings = fsyst.subgraphs[subgraph]
                    hop = (i - site_offsets[to_w], j - site_offsets[from_w])
                    hop = np.array(hop, dtype=dtype)
                    hoppings = hoppings.transpose().view(dtype).reshape(-1)
                    selector = np.recarray.searchsorted(hoppings, hop)
                    assert np.allclose(value[selector], hops[tail, head])
354
        else:
355
356
357
358
359
360
            value = fsyst.hoppings[edge_id][0]
            if value is builder.Other:
                assert (head, tail) in hops
            else:
                assert (tail, head) in hops
                assert value is hops[tail, head]
361

362
363
def check_id_by_site(fsyst):
    for i, site in enumerate(fsyst.sites):
364
        assert fsyst.id_by_site[site] == i
365

Anton Akhmerov's avatar
Anton Akhmerov committed
366

367
368
@pytest.mark.parametrize("vectorize", [False, True])
def test_finalization(vectorize):
369
370
371
372
373
374
375
376
377
378
379
380
381
382
    """Test the finalization of finite and infinite systems.

    In order to exactly verify the finalization, low-level features of the
    build module are used directly.  This is not the way one would use a
    finalized system in normal code.
    """
    def set_sites(dest):
        while len(dest) < n_sites:
            site = rng.randrange(size), rng.randrange(size)
            if site not in dest:
                dest[site] = random_onsite_hamiltonian(rng)

    def set_hops(dest, sites):
        while len(dest) < n_hops:
383
            a, b = rng.sample(list(sites), 2)
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
            if (a, b) not in dest and (b, a) not in dest:
                dest[a, b] = random_hopping_integral(rng)

    rng = Random(123)
    size = 20
    n_sites = 120
    n_hops = 500

    # Make scattering region blueprint.
    sr_sites = {}
    set_sites(sr_sites)
    sr_hops = {}
    set_hops(sr_hops, sr_sites)

    # Make lead blueprint.
    possible_neighbors = rng.sample(list(sr_sites), n_sites // 2)
    lead_sites = {}
    for pn in possible_neighbors:
        lead_sites[pn] = random_hopping_integral(rng)
    set_sites(lead_sites)
    lead_hops = {}        # Hoppings within a single lead unit cell
    set_hops(lead_hops, lead_sites)
    lead_sites_list = list(lead_sites)
    neighbors = set()
Joseph Weston's avatar
Joseph Weston committed
408
    for i in range(n_hops):
409
410
411
412
413
414
415
416
417
418
419
420
421
        while True:
            a = rng.choice(lead_sites_list)
            b = rng.choice(possible_neighbors)
            neighbors.add(b)
            b = b[0] - size, b[1]
            if rng.randrange(2):
                a, b = b, a
            if (a, b) not in lead_hops and (b, a) not in lead_hops:
                break
        lead_hops[a, b] = random_hopping_integral(rng)
    neighbors = sorted(neighbors)

    # Build scattering region from blueprint and test it.
422
    syst = builder.Builder(vectorize=vectorize)
423
    fam = kwant.lattice.general(ta.identity(2), norbs=1)
Joseph Weston's avatar
Joseph Weston committed
424
    for site, value in sr_sites.items():
425
        syst[fam(*site)] = value
Joseph Weston's avatar
Joseph Weston committed
426
    for hop, value in sr_hops.items():
427
428
        syst[fam(*hop[0]), fam(*hop[1])] = value
    fsyst = syst.finalized()
429
    check_id_by_site(fsyst)
430
431
    check_onsite(fsyst, sr_sites)
    check_hoppings(fsyst, sr_hops)
432
    # check that sites are sorted
433
    assert tuple(fsyst.sites) == tuple(sorted(fam(*site) for site in sr_sites))
434
435

    # Build lead from blueprint and test it.
436
    lead = builder.Builder(kwant.TranslationalSymmetry((size, 0)))
Joseph Weston's avatar
Joseph Weston committed
437
    for site, value in lead_sites.items():
438
439
        shift = rng.randrange(-5, 6) * size
        site = site[0] + shift, site[1]
440
        lead[fam(*site)] = value
441
442
443
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always")
        lead.finalized()        # Trigger the warning.
444
        assert len(w) == 1
445
446
        assert issubclass(w[0].category, RuntimeWarning)
        assert "disconnected" in str(w[0].message)
Joseph Weston's avatar
Joseph Weston committed
447
    for (a, b), value in lead_hops.items():
448
449
450
        shift = rng.randrange(-5, 6) * size
        a = a[0] + shift, a[1]
        b = b[0] + shift, b[1]
451
        lead[fam(*a), fam(*b)] = value
452
453
454
    flead = lead.finalized()
    all_sites = list(lead_sites)
    all_sites.extend((x - size, y) for (x, y) in neighbors)
455
    check_id_by_site(fsyst)
456
457
458
459
    check_onsite(flead, all_sites, check_values=False)
    check_onsite(flead, lead_sites, subset=True)
    check_hoppings(flead, lead_hops)

460
    # Attach lead to system with empty interface.
461
    syst.leads.append(builder.BuilderLead(lead, ()))
Anton Akhmerov's avatar
Anton Akhmerov committed
462
    raises(ValueError, syst.finalized)
463
464

    # Attach lead with improper interface.
465
    syst.leads[-1] = builder.BuilderLead(
466
        lead, 2 * tuple(system.Site(fam, n) for n in neighbors))
Anton Akhmerov's avatar
Anton Akhmerov committed
467
    raises(ValueError, syst.finalized)
468
469

    # Attach lead properly.
470
    syst.leads[-1] = builder.BuilderLead(
471
        lead, (system.Site(fam, n) for n in neighbors))
472
    fsyst = syst.finalized()
473
474
475
    assert len(fsyst.lead_interfaces) == 1
    assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] ==
            neighbors)
476
477

    # test that we cannot finalize a system with a badly sorted interface order
478
    raises(ValueError, builder.InfiniteSystem, lead,
479
           [system.Site(fam, n) for n in reversed(neighbors)])
480
    # site ordering independent of whether interface was specified
481
    flead_order = builder.InfiniteSystem(lead, [system.Site(fam, n)
482
                                                for n in neighbors])
483
    assert tuple(flead.sites) == tuple(flead_order.sites)
484
485


486
    syst.leads[-1] = builder.BuilderLead(
487
        lead, (system.Site(fam, n) for n in neighbors))
488
    fsyst = syst.finalized()
489
490
491
    assert len(fsyst.lead_interfaces) == 1
    assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] ==
            neighbors)
492

493
    # Add a hopping to the lead which couples two next-nearest cells and check
494
495
496
497
    # whether this leads to an error.
    a = rng.choice(lead_sites_list)
    b = rng.choice(possible_neighbors)
    b = b[0] + 2 * size, b[1]
498
    lead[fam(*a), fam(*b)] = random_hopping_integral(rng)
Anton Akhmerov's avatar
Anton Akhmerov committed
499
    raises(ValueError, lead.finalized)
500
501


502
503
504
505
506
507
508
509
510
511
512
513
514
def _make_system_from_sites(sites, vectorize):
    syst = builder.Builder(vectorize=vectorize)
    for s in sites:
        norbs = s.family.norbs
        if norbs:
            syst[s] = np.eye(s.family.norbs)
        else:
            syst[s] = None
    return syst.finalized()


@pytest.mark.parametrize("vectorize", [False, True])
def test_site_ranges(vectorize):
515
516
517
518
519
520
521
    lat1a = kwant.lattice.chain(norbs=1, name='a')
    lat1b = kwant.lattice.chain(norbs=1, name='b')
    lat2 = kwant.lattice.chain(norbs=2)

    # simple case -- single site family
    for lat in (lat1a, lat2):
        sites = list(map(lat, range(10)))
522
523
        syst = _make_system_from_sites(sites, vectorize)
        ranges = syst.site_ranges
524
        expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)]
525
        assert np.array_equal(ranges, expected)
526
527

    # pair of site families
528
529
530
531
    sites = it.chain(map(lat1a, range(4)), map(lat1b, range(6)))
    syst = _make_system_from_sites(sites, vectorize)
    expected = [(0, 1, 0), (4, 1, 4), (10, 0, 10)]
    assert np.array_equal(expected, syst.site_ranges)
532
533
    sites = it.chain(map(lat2, range(4)), map(lat1a, range(6)),
                     map(lat1b, range(4)))
534
    syst = _make_system_from_sites(sites, vectorize)
535
    expected = [(0, 2, 0), (4, 1, 4*2), (10, 1, 4*2+6), (14, 0, 4*2+10)]
536
    assert np.array_equal(expected, syst.site_ranges)
537
538
539
540

    # test with an actual builder
    for lat in (lat1a, lat2):
        sites = list(map(lat, range(10)))
541
        syst = kwant.Builder(vectorize=vectorize)
542
543
544
        syst[sites] = np.eye(lat.norbs)
        ranges = syst.finalized().site_ranges
        expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)]
545
546
547
548
549
550
551
552
553
554
555
556
557
        assert np.array_equal(ranges, expected)
        if not vectorize:  # vectorized systems *must* have all norbs
            # poison system with a single site with no norbs defined.
            # Also catch the deprecation warning.
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                syst[kwant.lattice.chain(norbs=None)(0)] = 1
            ranges = syst.finalized().site_ranges
            assert ranges is None


@pytest.mark.parametrize("vectorize", [False, True])
def test_hamiltonian_evaluation(vectorize):
558
559
560
561
562
563
564
    def f_onsite(site):
        return site.tag[0]

    def f_hopping(a, b):
        a, b = a.tag, b.tag
        return complex(a[0] + b[0], a[1] - b[1])

565
566
567
568
569
570
571
    def f_onsite_vectorized(sites):
        return sites.tags[:, 0]

    def f_hopping_vectorized(a, b):
        a, b = a.tags, b.tags
        return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1])

572
573
574
    tags = [(0, 0), (1, 1), (2, 2), (3, 3)]
    edges = [(0, 1), (0, 2), (0, 3), (1, 2)]

575
    syst = builder.Builder(vectorize=vectorize)
576
    fam = builder.SimpleSiteFamily(norbs=1)
577
    sites = [fam(*tag) for tag in tags]
578
579
580
581
582
583
584
    hoppings = [(sites[i], sites[j]) for i, j in edges]
    if vectorize:
        syst[sites] = f_onsite_vectorized
        syst[hoppings] = f_hopping_vectorized
    else:
        syst[sites] = f_onsite
        syst[hoppings] = f_hopping
585
    fsyst = syst.finalized()
586

587
588
    assert fsyst.graph.num_nodes == len(tags)
    assert fsyst.graph.num_edges == 2 * len(edges)
589
590

    for i in range(len(tags)):
591
        site = fsyst.sites[i]
592
        assert site in sites
593
594
595
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
            assert fsyst.hamiltonian(i, i) == f_onsite(site)
596

597
598
599
    for t, h in fsyst.graph:
        tsite = fsyst.sites[t]
        hsite = fsyst.sites[h]
600
601
602
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
            assert fsyst.hamiltonian(t, h) == f_hopping(tsite, hsite)
603

604
605
606
607
608
609
610
611
612
613
    # test when user-function raises errors
    def onsite_raises(site):
        raise ValueError()

    def hopping_raises(a, b):
        raise ValueError('error message')

    def test_raising(fsyst, hop):
        a, b = hop
        # exceptions are converted to kwant.UserCodeError and we add our message
614
        with raises(kwant.UserCodeError) as exc_info:
615
616
617
            with warnings.catch_warnings():
                warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
                fsyst.hamiltonian(a, a)
618
        msg = 'Error occurred in user-supplied value function "onsite_raises"'
619
        assert msg in str(exc_info.value)
620
621

        for hop in [(a, b), (b, a)]:
622
            with raises(kwant.UserCodeError) as exc_info:
623
624
625
626
                with warnings.catch_warnings():
                    # Ignore deprecation warnings for 'hamiltonian'
                    warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
                    fsyst.hamiltonian(*hop)
627
628
            msg = ('Error occurred in user-supplied '
                   'value function "hopping_raises"')
629
            assert msg in str(exc_info.value)
630
631
632
633
634
635

    # test with finite system
    new_hop = (fam(-1, 0), fam(0, 0))
    syst[new_hop[0]] = onsite_raises
    syst[new_hop] = hopping_raises
    fsyst = syst.finalized()
636
    hop = tuple(map(fsyst.id_by_site.__getitem__, new_hop))
637
638
639
640
    test_raising(fsyst, hop)

    # test with infinite system
    inf_syst = kwant.Builder(VerySimpleSymmetry(2))
641
642
    for k, v in it.chain(syst.site_value_pairs(), syst.hopping_value_pairs()):
        inf_syst[k] = v
643
    inf_fsyst = inf_syst.finalized()
644
    hop = tuple(map(inf_fsyst.id_by_site.__getitem__, new_hop))
645
646
    test_raising(inf_fsyst, hop)

647

648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
def test_vectorized_hamiltonian_evaluation():

    def onsite(site):
        return site.tag[0]

    def vectorized_onsite(sites):
        return sites.tags[:, 0]

    def hopping(to_site, from_site):
        a, b = to_site.tag, from_site.tag
        return a[0] + b[0] + 1j * (a[1] - b[1])

    def vectorized_hopping(to_sites, from_sites):
        a, b = to_sites.tags, from_sites.tags
        return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1])

    tags = [(0, 0), (1, 1), (2, 2), (3, 3)]
    edges = [(0, 1), (0, 2), (0, 3), (1, 2)]

    fam = builder.SimpleSiteFamily(norbs=1)
    sites = [fam(*tag) for tag in tags]
    hops = [(fam(*tags[i]), fam(*tags[j])) for (i, j) in edges]

    syst_simple = builder.Builder(vectorize=False)
    syst_simple[sites] = onsite
    syst_simple[hops] = hopping
    fsyst_simple = syst_simple.finalized()

    syst_vectorized = builder.Builder(vectorize=True)
    syst_vectorized[sites] = vectorized_onsite
    syst_vectorized[hops] = vectorized_hopping
    fsyst_vectorized = syst_vectorized.finalized()

    assert fsyst_vectorized.graph.num_nodes == len(tags)
    assert fsyst_vectorized.graph.num_edges == 2 * len(edges)
    assert len(fsyst_vectorized.site_arrays) == 1
    assert fsyst_vectorized.site_arrays[0] == system.SiteArray(fam, tags)

    assert np.allclose(
        fsyst_simple.hamiltonian_submatrix(),
        fsyst_vectorized.hamiltonian_submatrix(),
    )

    for i in range(len(tags)):
        site = fsyst_vectorized.sites[i]
        assert site in sites
694
695
696
697
698
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
            assert (
                fsyst_vectorized.hamiltonian(i, i)
                == fsyst_simple.hamiltonian(i, i))
699
700

    for t, h in fsyst_vectorized.graph:
701
702
703
704
705
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
            assert (
                fsyst_vectorized.hamiltonian(t, h)
                == fsyst_simple.hamiltonian(t, h))
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742

    # Test infinite system, including hoppings that go both ways into
    # the next cell
    lat = kwant.lattice.square(norbs=1)

    syst_vectorized = builder.Builder(kwant.TranslationalSymmetry((-1, 0)),
                                      vectorize=True)
    syst_vectorized[lat(0, 0)] = 4
    syst_vectorized[lat(0, 1)] = 5
    syst_vectorized[lat(0, 2)] = vectorized_onsite
    syst_vectorized[(lat(1, 0), lat(0, 0))] = 1j
    syst_vectorized[(lat(2, 1), lat(1, 1))] = vectorized_hopping
    fsyst_vectorized = syst_vectorized.finalized()

    syst_simple = builder.Builder(kwant.TranslationalSymmetry((-1, 0)),
                                      vectorize=False)
    syst_simple[lat(0, 0)] = 4
    syst_simple[lat(0, 1)] = 5
    syst_simple[lat(0, 2)] = onsite
    syst_simple[(lat(1, 0), lat(0, 0))] = 1j
    syst_simple[(lat(2, 1), lat(1, 1))] = hopping
    fsyst_simple = syst_simple.finalized()

    assert np.allclose(
        fsyst_vectorized.hamiltonian_submatrix(),
        fsyst_simple.hamiltonian_submatrix(),
    )
    assert np.allclose(
        fsyst_vectorized.cell_hamiltonian(),
        fsyst_simple.cell_hamiltonian(),
    )
    assert np.allclose(
        fsyst_vectorized.inter_cell_hopping(),
        fsyst_simple.inter_cell_hopping(),
    )


743
744
745
746
747
@pytest.mark.parametrize("sym", [
    builder.NoSymmetry(),
    kwant.TranslationalSymmetry([-1]),
])
def test_vectorized_requires_norbs(sym):
748
749
750
751

    # Catch deprecation warning for lack of norbs
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
752
        lat = kwant.lattice.chain()
753

754
    syst = builder.Builder(sym, vectorize=True)
755
756
    syst[lat(0)] = syst[lat(1)] = 2
    syst[lat(1), lat(0)] = -1
757
758
759
760

    raises(ValueError, syst.finalized)


761
762
763
764
765
def test_dangling():
    def make_system():
        #        1
        #       / \
        #    3-0---2-4-5  6-7  8
766
        syst = builder.Builder()
767
        fam = builder.SimpleSiteFamily(norbs=1)
768
769
770
771
772
773
774
        syst[(fam(i) for i in range(9))] = None
        syst[[(fam(0), fam(1)), (fam(1), fam(2)), (fam(2), fam(0))]] = None
        syst[[(fam(0), fam(3)), (fam(2), fam(4)), (fam(4), fam(5))]] = None
        syst[fam(6), fam(7)] = None
        return syst

    syst0 = make_system()
775
776
    assert (sorted(site.tag for site in syst0.dangling()) ==
            sorted([(3,), (5,), (6,), (7,), (8,)]))
777
    syst0.eradicate_dangling()
778

779
    syst1 = make_system()
780
    while True:
781
        dangling = list(syst1.dangling())
Anton Akhmerov's avatar
Anton Akhmerov committed
782
783
        if not dangling:
            break
784
        del syst1[dangling]
785

786
787
788
789
    assert (sorted(site.tag for site in syst0.sites()) ==
            sorted([(0,), (1,), (2,)]))
    assert (sorted(site.tag for site in syst0.sites()) ==
            sorted(site.tag for site in syst1.sites()))
790
791
792


def test_builder_with_symmetry():
793
    g = kwant.lattice.general(ta.identity(3), norbs=1)
794
    sym = kwant.TranslationalSymmetry((0, 0, 3), (0, 2, 0))
795
    syst = builder.Builder(sym)
796
797

    t, V = 1.0j, 0.0
798
799
800
801
802
803
804
805
806
807
808
    hoppings = [(g(5, 0, 0), g(0, 5, 0)),
                (g(0, 5, 0), g(0, 0, 5)),
                (g(0, 0, 5), g(5, 0, 0)),
                (g(0, 3, 0), g(0, 0, 5)),
                (g(0, 7, -6), g(5, 6, -6))]
    hoppings_fd = [(g(5, 0, 0), g(0, 5, 0)),
                   (g(0, 1, 0), g(0, -4, 5)),
                   (g(0, 0, 2), g(5, 0, -3)),
                   (g(0, 1, 0), g(0, -2, 5)),
                   (g(0, 1, 0), g(5, 0, 0))]

809
810
    syst[(a for a, b in hoppings)] = V
    syst[hoppings] = t
811

812
    # TODO: Once Tinyarray supports "<" the conversion to tuple can be removed.
813
814
    assert (sorted(tuple(site.tag) for site in syst.sites()) ==
            sorted(set(tuple(hop[0].tag) for hop in hoppings_fd)))
815
816
    for sites in hoppings_fd:
        for site in sites:
817
            assert site in syst
818
            assert syst[site] == V
819

820
    # TODO: Once Tinyarray supports "<" the conversion to tuple can be removed.
821
822
    assert (sorted((tuple(a.tag), tuple(b.tag)) for a, b in syst.hoppings()) ==
            sorted((tuple(a.tag), tuple(b.tag)) for a, b in hoppings_fd))
823
824
    for hop in hoppings_fd:
        rhop = hop[1], hop[0]
825
826
        assert hop in syst
        assert rhop in syst
827
828
        assert syst[hop] == t
        assert syst[rhop] == t.conjugate()
829

830
831
    del syst[g(0, 6, -4), g(0, 11, -9)]
    assert (g(0, 1, 0), g(0, -4, 5)) not in syst
832

833
    del syst[g(0, 3, -3)]
834
835
    assert (list((a.tag, b.tag) for a, b in syst.hoppings()) == [((0, 0, 2),
                                                                  (5, 0, -3))])
836
837


Anton Akhmerov's avatar
Anton Akhmerov committed
838
def test_fill():
839
    g = kwant.lattice.square(norbs=1)
Anton Akhmerov's avatar
Anton Akhmerov committed
840
841
842
843
    sym_x = kwant.TranslationalSymmetry((-1, 0))
    sym_xy = kwant.TranslationalSymmetry((-1, 0), (0, 1))

    template_1d = builder.Builder(sym_x)
844
845
    template_1d[g(0, 0)] = None
    template_1d[g.neighbors()] = None
Anton Akhmerov's avatar
Anton Akhmerov committed
846
847
848
849

    def line_200(site):
        return -100 <= site.pos[0] < 100

850
    ## Test that copying a builder by "fill" preserves everything.
851
852
853
854
    for sym, func in [(kwant.TranslationalSymmetry(*np.diag([3, 4, 5])),
                       lambda pos: True),
                      (builder.NoSymmetry(),
                       lambda pos: ta.dot(pos, pos) < 17)]:
855
        cubic = kwant.lattice.general(ta.identity(3), norbs=1)
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880

        # Make a weird system.
        orig = kwant.Builder(sym)
        sites = cubic.shape(func, (0, 0, 0))
        for i, site in enumerate(orig.expand(sites)):
            if i % 7 == 0:
                continue
            orig[site] = i
        for i, hopp in enumerate(orig.expand(cubic.neighbors(1))):
            if i % 11 == 0:
                continue
            orig[hopp] = i * 1.2345
        for i, hopp in enumerate(orig.expand(cubic.neighbors(2))):
            if i % 13 == 0:
                continue
            orig[hopp] = i * 1j

        # Clone the original using fill.
        clone = kwant.Builder(sym)
        clone.fill(orig, lambda s: True, (0, 0, 0))

        # Verify that both are identical.
        assert set(clone.site_value_pairs()) == set(orig.site_value_pairs())
        assert (set(clone.hopping_value_pairs())
                == set(orig.hopping_value_pairs()))
881

882
    ## Test for warning when "start" is outside the filling shape.
883
    target = builder.Builder()
Christoph Groth's avatar
Christoph Groth committed
884
885
886
    for start in [(-101, 0), (101, 0)]:
        with warns(RuntimeWarning):
            target.fill(template_1d, line_200, start)
887

Christoph Groth's avatar
Christoph Groth committed
888
889
890
    ## Test filling of infinite builder.
    for n in [1, 2, 4]:
        sym_n = kwant.TranslationalSymmetry((n, 0))
891
892
893
894
895
896
897
        for start in [g(0, 0), g(20, 0)]:
            target = builder.Builder(sym_n)
            sites = target.fill(template_1d, lambda s: True, start,
                                max_sites=10)
            assert len(sites) == n
            assert len(list(target.hoppings())) == n
            assert set(sym_n.to_fd(s) for s in sites) == set(target.sites())
Christoph Groth's avatar
Christoph Groth committed
898

Anton Akhmerov's avatar
Anton Akhmerov committed
899
900
901
902
    ## test max_sites
    target = builder.Builder()
    for max_sites in (-1, 0):
        with raises(ValueError):
903
904
            target.fill(template_1d, lambda site: True, g(0, 0),
                        max_sites=max_sites)
905
        assert len(list(target.sites())) == 0
Anton Akhmerov's avatar
Anton Akhmerov committed
906
907
    target = builder.Builder()
    with raises(RuntimeError):
908
        target.fill(template_1d, line_200, g(0, 0) , max_sites=10)
Anton Akhmerov's avatar
Anton Akhmerov committed
909
910
    ## test filling
    target = builder.Builder()
911
    added_sites = target.fill(template_1d, line_200, g(0, 0))
Anton Akhmerov's avatar
Anton Akhmerov committed
912
    assert len(added_sites) == 200
913
    # raise warning if target already contains all starting sites
Christoph Groth's avatar
Christoph Groth committed
914
    with warns(RuntimeWarning):
915
        target.fill(template_1d, line_200, g(0, 0))
Anton Akhmerov's avatar
Anton Akhmerov committed
916
917
918
919
920

    ## test multiplying unit cell size in 1D
    n_cells = 10
    sym_nx = kwant.TranslationalSymmetry(*(sym_x.periods * n_cells))
    target = builder.Builder(sym_nx)
921
    target.fill(template_1d, lambda site: True, g(0, 0))
Anton Akhmerov's avatar
Anton Akhmerov committed
922
923

    should_be_syst = builder.Builder(sym_nx)
924
925
    should_be_syst[(g(i, 0) for i in range(n_cells))] = None
    should_be_syst[g.neighbors()] = None
Anton Akhmerov's avatar
Anton Akhmerov committed
926
927
928
929
930
931

    assert sorted(target.sites()) == sorted(should_be_syst.sites())
    assert sorted(target.hoppings()) == sorted(should_be_syst.hoppings())

    ## test multiplying unit cell size in 2D
    template_2d = builder.Builder(sym_xy)
932
933
934
    template_2d[g(0, 0)] = None
    template_2d[g.neighbors()] = None
    template_2d[builder.HoppingKind((2, 2), g)] = None
Anton Akhmerov's avatar
Anton Akhmerov committed
935
936
937
938

    nm_cells = (3, 5)
    sym_nmxy = kwant.TranslationalSymmetry(*(sym_xy.periods * nm_cells))
    target = builder.Builder(sym_nmxy)
939
    target.fill(template_2d, lambda site: True, g(0, 0))
Anton Akhmerov's avatar
Anton Akhmerov committed
940
941

    should_be_syst = builder.Builder(sym_nmxy)
942
943
944
    should_be_syst[(g(i, j) for i in range(10) for j in range(10))] = None
    should_be_syst[g.neighbors()] = None
    should_be_syst[builder.HoppingKind((2, 2), g)] = None
Anton Akhmerov's avatar
Anton Akhmerov committed
945
946
947
948
949
950
951
952
953
954

    assert sorted(target.sites()) == sorted(should_be_syst.sites())
    assert sorted(target.hoppings()) == sorted(should_be_syst.hoppings())

    ## test filling 0D builder with 2D builder
    def square_shape(site):
        x, y = site.tag
        return 0 <= x < 10 and 0 <= y < 10

    target = builder.Builder()
955
    target.fill(template_2d, square_shape, g(0, 0))
Anton Akhmerov's avatar
Anton Akhmerov committed
956
957

    should_be_syst = builder.Builder()
958
959
960
    should_be_syst[(g(i, j) for i in range(10) for j in range(10))] = None
    should_be_syst[g.neighbors()] = None
    should_be_syst[builder.HoppingKind((2, 2), g)] = None
Anton Akhmerov's avatar
Anton Akhmerov committed
961
962
963
964

    assert sorted(target.sites()) == sorted(should_be_syst.sites())
    assert sorted(target.hoppings()) == sorted(should_be_syst.hoppings())

965
    ## test that 'fill' respects the symmetry of the target builder
966
    lat = kwant.lattice.chain(a=1, norbs=1)
967
968
969
970
971
972
973
    template = builder.Builder(kwant.TranslationalSymmetry((-1,)))
    template[lat(0)] = 2
    template[lat.neighbors()] = -1

    target = builder.Builder(kwant.TranslationalSymmetry((-2,)))
    target[lat(0)] = None
    to_target_fd = target.symmetry.to_fd
974
975
    # Refuses to fill the target because target already contains the starting
    # site.
Christoph Groth's avatar
Christoph Groth committed
976
977
    with warns(RuntimeWarning):
        target.fill(template, lambda x: True, lat(0))