# Copyright 2011-2019 Kwant authors. # # This file is part of Kwant. It is subject to the license terms in the file # LICENSE.rst found in the top-level directory of this distribution and at # http://kwant-project.org/license. A list of Kwant authors can be found in # the file AUTHORS.rst at the top-level directory of this distribution and at # http://kwant-project.org/authors. import warnings import pickle import itertools as it import functools as ft from random import Random import numpy as np import tinyarray as ta import pytest from pytest import raises, warns from numpy.testing import assert_almost_equal import kwant from kwant import builder, system from kwant._common import ensure_rng, KwantDeprecationWarning @ft.total_ordering class SimpleSiteFamily(system.SiteFamily): """A site family used as an example and for testing. A family of sites tagged by any python objects where object satisfied condition ``object == eval(repr(object))``. It exists to provide a basic site family that can be used for testing the 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 `SimpleSiteFamily` when `kwant.lattice.Monatomic` would also work. """ def __init__(self, name=None, norbs=None): canonical_repr = '{0}({1}, {2})'.format(self.__class__, repr(name), repr(norbs)) super().__init__(canonical_repr, name, norbs) 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 def test_bad_keys(): def setitem(key): syst[key] = None fam = SimpleSiteFamily(norbs=1) syst = builder.Builder() failures = [ # Invalid single keys ([syst.__contains__, syst.__getitem__, setitem, syst.__delitem__], [(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 ([syst.__getitem__, setitem, syst.__delitem__], [(KeyError, [(fam(0), fam(3)), (fam(2), fam(1)), (fam(2), fam(3))])]), # Sequences containing a bad key. ([setitem, syst.__delitem__], [(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 ([syst.__delitem__], [(KeyError, [fam(123), [fam(0), fam(123)], [fam(123), fam(1)]])]), # Various things that are not sites present in the system. ([syst.degree, lambda site: list(syst.neighbors(site))], [(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: syst[[fam(0), fam(1)]] = None syst[fam(0), fam(1)] = None try: raises(error, func, key) except AssertionError: print(func, error, key) raise def test_site_families(): syst = builder.Builder() fam = SimpleSiteFamily(norbs=1) ofam = SimpleSiteFamily(norbs=1) yafam = SimpleSiteFamily('another_name', norbs=1) syst[fam(0)] = 7 assert syst[fam(0)] == 7 assert len(set([fam, ofam, fam('a'), ofam('a'), yafam])) == 3 syst[fam(1)] = 123 assert syst[fam(1)] == 123 assert syst[ofam(1)] == 123 raises(KeyError, syst.__getitem__, yafam(1)) # test site families compare equal/not-equal assert fam == ofam assert fam != yafam assert fam != None assert fam != 'a' # test site families sorting fam1 = SimpleSiteFamily(norbs=1) fam2 = SimpleSiteFamily(norbs=2) assert fam1 < fam2 # string '1' is lexicographically less than '2' class VerySimpleSymmetry(builder.Symmetry): def __init__(self, period): self.period = period @property def num_directions(self): return 1 def has_subgroup(self, other): 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) def which(self, site): return ta.array((site.tag[0] // self.period,), int) def act(self, element, a, b=None): shifted = lambda site, delta: site.family(*ta.add(site.tag, delta)) delta = (self.period * element[0],) + (len(a.tag) - 1) * (0,) if b is None: return shifted(a, delta) else: return shifted(a, delta), shifted(b, delta) # 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, unknown_hoppings, sym=None): fam = SimpleSiteFamily(norbs=1) syst = builder.Builder(sym) t, V = 1.0j, 0.0 syst[sites] = V for site in sites: syst[site] = V syst[hoppings] = t for hopping in hoppings: syst[hopping] = t for hopping in unknown_hoppings: raises(KeyError, syst.__setitem__, hopping, t) 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 for site in sites: assert site in syst assert syst[site] == V for hop in hoppings: rev_hop = hop[1], hop[0] assert hop in syst assert rev_hop in syst assert syst[hop] == t assert syst[rev_hop] == t.conjugate() assert syst.degree(sites[0]) == 2 assert (sorted(s for s in syst.neighbors(sites[0])) == sorted([sites[1], sites[-1]])) del syst[hoppings] assert list(syst.hoppings()) == [] syst[hoppings] = t del syst[sites[0]] 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])) 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())) def test_construction_and_indexing(): # Without symmetry fam = SimpleSiteFamily(norbs=1) 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))] unknown_hoppings = [(fam(0, 1), fam(7, 8)), (fam(12, 14), fam(0, 1))] check_construction_and_indexing(sites, sites, hoppings, hoppings, unknown_hoppings) # With symmetry 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))] unknown_hoppings = [(fam(0, 0), fam(0, 3)), (fam(0, 4), fam(0, 0)), (fam(0, 0), fam(2, 3)), (fam(2, 4), fam(0, 0)), (fam(4, 2), fam(6, 3)), (fam(6, 4), fam(4, 2))] sym = VerySimpleSymmetry(2) check_construction_and_indexing(sites, sites_fd, hoppings, hoppings_fd, unknown_hoppings, sym) def test_hermitian_conjugation(): def f(i, j, arg): i, j = i.tag, j.tag if j[0] == i[0] + 1: return arg * ta.array([[1, 2j], [3 + 1j, 4j]]) else: raise ValueError syst = builder.Builder() fam = SimpleSiteFamily(norbs=1) syst[fam(0)] = syst[fam(1)] = ta.identity(2) syst[fam(0), fam(1)] = f assert syst[fam(0), fam(1)] is f assert isinstance(syst[fam(1), fam(0)], builder.HermConjOfFunc) assert (syst[fam(1), fam(0)](fam(1), fam(0), 2) == syst[fam(0), fam(1)](fam(0), fam(1), 2).conjugate().transpose()) 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 def test_value_equality_and_identity(): m = ta.array([[1, 2], [3j, 4j]]) syst = builder.Builder() fam = SimpleSiteFamily(norbs=1) syst[fam(0)] = m syst[fam(1)] = m assert syst[fam(1)] is m syst[fam(0), fam(1)] = m assert syst[fam(1), fam(0)] == m.transpose().conjugate() assert syst[fam(0), fam(1)] is m syst[fam(1), fam(0)] = m assert syst[fam(0), fam(1)] == m.transpose().conjugate() assert syst[fam(1), fam(0)] is m def random_onsite_hamiltonian(rng): return 2 * rng.random() - 1 def random_hopping_integral(rng): return complex(2 * rng.random() - 1, 2 * rng.random() - 1) def check_onsite(fsyst, sites, subset=False, check_values=True): vectorized = system.is_vectorized(fsyst) if vectorized: site_offsets, _, _ = fsyst.site_ranges.transpose() freq = {} for node in range(fsyst.graph.num_nodes): site = fsyst.sites[node].tag freq[site] = freq.get(site, 0) + 1 if check_values and site in sites: if vectorized: term = fsyst._onsite_term_by_site_id[node] value = fsyst._term_values[term] if callable(value): assert value is sites[site] else: (w, _), (off, _) = fsyst.subgraphs[fsyst.terms[term].subgraph] node_off = node - site_offsets[w] selector = np.searchsorted(off, node_off) assert np.allclose(value[selector], sites[site]) else: assert fsyst.onsites[node][0] is sites[site] if not subset: # Check that all sites of `fsyst` are in `sites`. for site in freq.keys(): assert site in sites # Check that all sites of `sites` are in `fsyst`. for site in sites: assert freq[site] == 1 def check_hoppings(fsyst, hops): vectorized = system.is_vectorized(fsyst) if vectorized: site_offsets, _, _ = fsyst.site_ranges.transpose() assert fsyst.graph.num_edges == 2 * len(hops) for edge_id, edge in enumerate(fsyst.graph): i, j = edge tail = fsyst.sites[i].tag head = fsyst.sites[j].tag if vectorized: term = fsyst._hopping_term_by_edge_id[edge_id] if term < 0: # Hermitian conjugate assert (head, tail) in hops else: value = fsyst._term_values[term] assert (tail, head) in hops if callable(value): assert value is hops[tail, head] else: dtype = np.dtype([('f0', int), ('f1', int)]) subgraph = fsyst.terms[term].subgraph (to_w, from_w), hoppings = fsyst.subgraphs[subgraph] hop = (i - site_offsets[to_w], j - site_offsets[from_w]) hop = np.array(hop, dtype=dtype) hoppings = hoppings.transpose().view(dtype).reshape(-1) selector = np.recarray.searchsorted(hoppings, hop) assert np.allclose(value[selector], hops[tail, head]) else: value = fsyst.hoppings[edge_id][0] if value is builder.Other: assert (head, tail) in hops else: assert (tail, head) in hops assert value is hops[tail, head] def check_id_by_site(fsyst): for i, site in enumerate(fsyst.sites): assert fsyst.id_by_site[site] == i @pytest.mark.parametrize("vectorize", [False, True]) def test_finalization(vectorize): """Test the finalization of finite and infinite systems. In order to exactly verify the finalization, low-level features of the 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: a, b = rng.sample(list(sites), 2) 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() for i in range(n_hops): 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. syst = builder.Builder(vectorize=vectorize) fam = kwant.lattice.general(ta.identity(2), norbs=1) for site, value in sr_sites.items(): syst[fam(*site)] = value for hop, value in sr_hops.items(): syst[fam(*hop[0]), fam(*hop[1])] = value fsyst = syst.finalized() check_id_by_site(fsyst) check_onsite(fsyst, sr_sites) check_hoppings(fsyst, sr_hops) # check that sites are sorted assert tuple(fsyst.sites) == tuple(sorted(fam(*site) for site in sr_sites)) # Build lead from blueprint and test it. lead = builder.Builder(kwant.TranslationalSymmetry((size, 0))) for site, value in lead_sites.items(): shift = rng.randrange(-5, 6) * size site = site[0] + shift, site[1] lead[fam(*site)] = value with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always") lead.finalized() # Trigger the warning. assert len(w) == 1 assert issubclass(w[0].category, RuntimeWarning) assert "disconnected" in str(w[0].message) for (a, b), value in lead_hops.items(): shift = rng.randrange(-5, 6) * size a = a[0] + shift, a[1] b = b[0] + shift, b[1] lead[fam(*a), fam(*b)] = value flead = lead.finalized() all_sites = list(lead_sites) all_sites.extend((x - size, y) for (x, y) in neighbors) check_id_by_site(fsyst) check_onsite(flead, all_sites, check_values=False) check_onsite(flead, lead_sites, subset=True) check_hoppings(flead, lead_hops) # Attach lead to system with empty interface. syst.leads.append(builder.BuilderLead(lead, ())) raises(ValueError, syst.finalized) # Attach lead with improper interface. syst.leads[-1] = builder.BuilderLead( lead, 2 * tuple(system.Site(fam, n) for n in neighbors)) raises(ValueError, syst.finalized) # Attach lead properly. syst.leads[-1] = builder.BuilderLead( lead, (system.Site(fam, n) for n in neighbors)) fsyst = syst.finalized() assert len(fsyst.lead_interfaces) == 1 assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] == neighbors) # test that we cannot finalize a system with a badly sorted interface order raises(ValueError, builder.InfiniteSystem, lead, [system.Site(fam, n) for n in reversed(neighbors)]) # site ordering independent of whether interface was specified flead_order = builder.InfiniteSystem(lead, [system.Site(fam, n) for n in neighbors]) assert tuple(flead.sites) == tuple(flead_order.sites) syst.leads[-1] = builder.BuilderLead( lead, (system.Site(fam, n) for n in neighbors)) fsyst = syst.finalized() assert len(fsyst.lead_interfaces) == 1 assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] == neighbors) # Add a hopping to the lead which couples two next-nearest cells and check # whether this leads to an error. a = rng.choice(lead_sites_list) b = rng.choice(possible_neighbors) b = b[0] + 2 * size, b[1] lead[fam(*a), fam(*b)] = random_hopping_integral(rng) raises(ValueError, lead.finalized) def _make_system_from_sites(sites, vectorize): syst = builder.Builder(vectorize=vectorize) for s in sites: norbs = s.family.norbs if norbs: syst[s] = np.eye(s.family.norbs) else: syst[s] = None return syst.finalized() @pytest.mark.parametrize("vectorize", [False, True]) def test_site_ranges(vectorize): lat1a = kwant.lattice.chain(norbs=1, name='a') lat1b = kwant.lattice.chain(norbs=1, name='b') lat2 = kwant.lattice.chain(norbs=2) # simple case -- single site family for lat in (lat1a, lat2): sites = list(map(lat, range(10))) syst = _make_system_from_sites(sites, vectorize) ranges = syst.site_ranges expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)] assert np.array_equal(ranges, expected) # pair of site families sites = it.chain(map(lat1a, range(4)), map(lat1b, range(6))) syst = _make_system_from_sites(sites, vectorize) expected = [(0, 1, 0), (4, 1, 4), (10, 0, 10)] assert np.array_equal(expected, syst.site_ranges) sites = it.chain(map(lat2, range(4)), map(lat1a, range(6)), map(lat1b, range(4))) syst = _make_system_from_sites(sites, vectorize) expected = [(0, 2, 0), (4, 1, 4*2), (10, 1, 4*2+6), (14, 0, 4*2+10)] assert np.array_equal(expected, syst.site_ranges) # test with an actual builder for lat in (lat1a, lat2): sites = list(map(lat, range(10))) syst = kwant.Builder(vectorize=vectorize) syst[sites] = np.eye(lat.norbs) ranges = syst.finalized().site_ranges expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)] assert np.array_equal(ranges, expected) if not vectorize: # vectorized systems *must* have all norbs # poison system with a single site with no norbs defined. # Also catch the deprecation warning. with warnings.catch_warnings(): warnings.simplefilter("ignore") syst[kwant.lattice.chain(norbs=None)(0)] = 1 ranges = syst.finalized().site_ranges assert ranges is None @pytest.mark.parametrize("vectorize", [False, True]) def test_hamiltonian_evaluation(vectorize): def f_onsite(site): return site.tag[0] def f_hopping(a, b): a, b = a.tag, b.tag return complex(a[0] + b[0], a[1] - b[1]) def f_onsite_vectorized(sites): return sites.tags[:, 0] def f_hopping_vectorized(a, b): a, b = a.tags, b.tags return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1]) tags = [(0, 0), (1, 1), (2, 2), (3, 3)] edges = [(0, 1), (0, 2), (0, 3), (1, 2)] syst = builder.Builder(vectorize=vectorize) fam = SimpleSiteFamily(norbs=1) sites = [fam(*tag) for tag in tags] hoppings = [(sites[i], sites[j]) for i, j in edges] if vectorize: syst[sites] = f_onsite_vectorized syst[hoppings] = f_hopping_vectorized else: syst[sites] = f_onsite syst[hoppings] = f_hopping fsyst = syst.finalized() assert fsyst.graph.num_nodes == len(tags) assert fsyst.graph.num_edges == 2 * len(edges) for i in range(len(tags)): site = fsyst.sites[i] assert site in sites with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) assert fsyst.hamiltonian(i, i) == f_onsite(site) for t, h in fsyst.graph: tsite = fsyst.sites[t] hsite = fsyst.sites[h] with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) assert fsyst.hamiltonian(t, h) == f_hopping(tsite, hsite) # test when user-function raises errors def onsite_raises(site): 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 with raises(kwant.UserCodeError) as exc_info: with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) fsyst.hamiltonian(a, a) msg = 'Error occurred in user-supplied value function "onsite_raises"' assert msg in str(exc_info.value) for hop in [(a, b), (b, a)]: with raises(kwant.UserCodeError) as exc_info: with warnings.catch_warnings(): # Ignore deprecation warnings for 'hamiltonian' warnings.filterwarnings("ignore", category=KwantDeprecationWarning) fsyst.hamiltonian(*hop) msg = ('Error occurred in user-supplied ' 'value function "hopping_raises"') assert msg in str(exc_info.value) # 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() hop = tuple(map(fsyst.id_by_site.__getitem__, new_hop)) test_raising(fsyst, hop) # test with infinite system inf_syst = kwant.Builder(VerySimpleSymmetry(2)) for k, v in it.chain(syst.site_value_pairs(), syst.hopping_value_pairs()): inf_syst[k] = v inf_fsyst = inf_syst.finalized() hop = tuple(map(inf_fsyst.id_by_site.__getitem__, new_hop)) test_raising(inf_fsyst, hop) def test_vectorized_hamiltonian_evaluation(): def onsite(site): return site.tag[0] def vectorized_onsite(sites): return sites.tags[:, 0] def hopping(to_site, from_site): a, b = to_site.tag, from_site.tag return a[0] + b[0] + 1j * (a[1] - b[1]) def vectorized_hopping(to_sites, from_sites): a, b = to_sites.tags, from_sites.tags return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1]) tags = [(0, 0), (1, 1), (2, 2), (3, 3)] edges = [(0, 1), (0, 2), (0, 3), (1, 2)] fam = SimpleSiteFamily(norbs=1) sites = [fam(*tag) for tag in tags] hops = [(fam(*tags[i]), fam(*tags[j])) for (i, j) in edges] syst_simple = builder.Builder(vectorize=False) syst_simple[sites] = onsite syst_simple[hops] = hopping fsyst_simple = syst_simple.finalized() syst_vectorized = builder.Builder(vectorize=True) syst_vectorized[sites] = vectorized_onsite syst_vectorized[hops] = vectorized_hopping fsyst_vectorized = syst_vectorized.finalized() assert fsyst_vectorized.graph.num_nodes == len(tags) assert fsyst_vectorized.graph.num_edges == 2 * len(edges) assert len(fsyst_vectorized.site_arrays) == 1 assert fsyst_vectorized.site_arrays[0] == system.SiteArray(fam, tags) assert np.allclose( fsyst_simple.hamiltonian_submatrix(), fsyst_vectorized.hamiltonian_submatrix(), ) for i in range(len(tags)): site = fsyst_vectorized.sites[i] assert site in sites with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) assert ( fsyst_vectorized.hamiltonian(i, i) == fsyst_simple.hamiltonian(i, i)) for t, h in fsyst_vectorized.graph: with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) assert ( fsyst_vectorized.hamiltonian(t, h) == fsyst_simple.hamiltonian(t, h)) # Test infinite system, including hoppings that go both ways into # the next cell lat = kwant.lattice.square(norbs=1) syst_vectorized = builder.Builder(kwant.TranslationalSymmetry((-1, 0)), vectorize=True) syst_vectorized[lat(0, 0)] = 4 syst_vectorized[lat(0, 1)] = 5 syst_vectorized[lat(0, 2)] = vectorized_onsite syst_vectorized[(lat(1, 0), lat(0, 0))] = 1j syst_vectorized[(lat(2, 1), lat(1, 1))] = vectorized_hopping fsyst_vectorized = syst_vectorized.finalized() syst_simple = builder.Builder(kwant.TranslationalSymmetry((-1, 0)), vectorize=False) syst_simple[lat(0, 0)] = 4 syst_simple[lat(0, 1)] = 5 syst_simple[lat(0, 2)] = onsite syst_simple[(lat(1, 0), lat(0, 0))] = 1j syst_simple[(lat(2, 1), lat(1, 1))] = hopping fsyst_simple = syst_simple.finalized() assert np.allclose( fsyst_vectorized.hamiltonian_submatrix(), fsyst_simple.hamiltonian_submatrix(), ) assert np.allclose( fsyst_vectorized.cell_hamiltonian(), fsyst_simple.cell_hamiltonian(), ) assert np.allclose( fsyst_vectorized.inter_cell_hopping(), fsyst_simple.inter_cell_hopping(), ) def test_vectorized_value_normalization(): # Here we test whether all legal shapes for values for vectorized Builders # are accepted: # + single scalars are interpreted as 1x1 matrices, and broadcast into an # (N, 1, 1) array # + single (n, m) matrices are broadcast into an (N, n, m) array # + a shape (N,) array is interpreted as an (N, 1, 1) array sz = np.array([[1, 0], [0, -1]]) scalars = [ 2, lambda s: 2, lambda s: [2] * len(s) ] matrices = [ sz, lambda s: sz, lambda s: [sz] * len(s) ] inter_lattice_matrices =[ [[1, 0]], lambda a, b: [[1, 0]], lambda a, b: [[[1, 0]]] * len(a) ] lat_a = kwant.lattice.chain(norbs=1) lat_b = kwant.lattice.chain(norbs=2) hams = [] for s, m, v in it.product(scalars, matrices, inter_lattice_matrices): syst = builder.Builder(vectorize=True) syst[map(lat_a, range(10))] = s syst[map(lat_b, range(10))] = m syst[((lat_a(i), lat_b(i)) for i in range(10))] = v h = syst.finalized().hamiltonian_submatrix() hams.append(h) assert np.all(hams == hams[0]) @pytest.mark.parametrize("sym", [ builder.NoSymmetry(), kwant.TranslationalSymmetry([-1]), ]) def test_vectorized_requires_norbs(sym): # Catch deprecation warning for lack of norbs with warnings.catch_warnings(): warnings.simplefilter("ignore") lat = kwant.lattice.chain() syst = builder.Builder(sym, vectorize=True) syst[lat(0)] = syst[lat(1)] = 2 syst[lat(1), lat(0)] = -1 raises(ValueError, syst.finalized) def test_dangling(): def make_system(): # 1 # / \ # 3-0---2-4-5 6-7 8 syst = builder.Builder() fam = SimpleSiteFamily(norbs=1) 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() assert (sorted(site.tag for site in syst0.dangling()) == sorted([(3,), (5,), (6,), (7,), (8,)])) syst0.eradicate_dangling() syst1 = make_system() while True: dangling = list(syst1.dangling()) if not dangling: break del syst1[dangling] 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())) def test_builder_with_symmetry(): g = kwant.lattice.general(ta.identity(3), norbs=1) sym = kwant.TranslationalSymmetry((0, 0, 3), (0, 2, 0)) syst = builder.Builder(sym) t, V = 1.0j, 0.0 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))] syst[(a for a, b in hoppings)] = V syst[hoppings] = t # TODO: Once Tinyarray supports "<" the conversion to tuple can be removed. assert (sorted(tuple(site.tag) for site in syst.sites()) == sorted(set(tuple(hop[0].tag) for hop in hoppings_fd))) for sites in hoppings_fd: for site in sites: assert site in syst assert syst[site] == V # TODO: Once Tinyarray supports "<" the conversion to tuple can be removed. 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)) for hop in hoppings_fd: rhop = hop[1], hop[0] assert hop in syst assert rhop in syst assert syst[hop] == t assert syst[rhop] == t.conjugate() del syst[g(0, 6, -4), g(0, 11, -9)] assert (g(0, 1, 0), g(0, -4, 5)) not in syst del syst[g(0, 3, -3)] assert (list((a.tag, b.tag) for a, b in syst.hoppings()) == [((0, 0, 2), (5, 0, -3))]) def test_fill(): g = kwant.lattice.square(norbs=1) sym_x = kwant.TranslationalSymmetry((-1, 0)) sym_xy = kwant.TranslationalSymmetry((-1, 0), (0, 1)) template_1d = builder.Builder(sym_x) template_1d[g(0, 0)] = None template_1d[g.neighbors()] = None def line_200(site): return -100 <= site.pos[0] < 100 ## Test that copying a builder by "fill" preserves everything. for sym, func in [(kwant.TranslationalSymmetry(*np.diag([3, 4, 5])), lambda pos: True), (builder.NoSymmetry(), lambda pos: ta.dot(pos, pos) < 17)]: cubic = kwant.lattice.general(ta.identity(3), norbs=1) # 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())) ## Test for warning when "start" is outside the filling shape. target = builder.Builder() for start in [(-101, 0), (101, 0)]: with warns(RuntimeWarning): target.fill(template_1d, line_200, start) ## Test filling of infinite builder. for n in [1, 2, 4]: sym_n = kwant.TranslationalSymmetry((n, 0)) 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()) ## test max_sites target = builder.Builder() for max_sites in (-1, 0): with raises(ValueError): target.fill(template_1d, lambda site: True, g(0, 0), max_sites=max_sites) assert len(list(target.sites())) == 0 target = builder.Builder() with raises(RuntimeError): target.fill(template_1d, line_200, g(0, 0) , max_sites=10) ## test filling target = builder.Builder() added_sites = target.fill(template_1d, line_200, g(0, 0)) assert len(added_sites) == 200 # raise warning if target already contains all starting sites with warns(RuntimeWarning): target.fill(template_1d, line_200, g(0, 0)) ## test multiplying unit cell size in 1D n_cells = 10 sym_nx = kwant.TranslationalSymmetry(*(sym_x.periods * n_cells)) target = builder.Builder(sym_nx) target.fill(template_1d, lambda site: True, g(0, 0)) should_be_syst = builder.Builder(sym_nx) should_be_syst[(g(i, 0) for i in range(n_cells))] = None should_be_syst[g.neighbors()] = None 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) template_2d[g(0, 0)] = None template_2d[g.neighbors()] = None template_2d[builder.HoppingKind((2, 2), g)] = None nm_cells = (3, 5) sym_nmxy = kwant.TranslationalSymmetry(*(sym_xy.periods * nm_cells)) target = builder.Builder(sym_nmxy) target.fill(template_2d, lambda site: True, g(0, 0)) should_be_syst = builder.Builder(sym_nmxy) 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 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() target.fill(template_2d, square_shape, g(0, 0)) should_be_syst = builder.Builder() 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 assert sorted(target.sites()) == sorted(should_be_syst.sites()) assert sorted(target.hoppings()) == sorted(should_be_syst.hoppings()) ## test that 'fill' respects the symmetry of the target builder lat = kwant.lattice.chain(a=1, norbs=1) 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 # Refuses to fill the target because target already contains the starting # site. with warns(RuntimeWarning): target.fill(template, lambda x: True, lat(0)) # should only add a single site (and hopping) new_sites = target.fill(template, lambda x: True, lat(1)) assert target[lat(0)] is None # should not be overwritten by template assert target[lat(-1)] == template[lat(0)] assert len(new_sites) == 1 assert to_target_fd(new_sites[0]) == to_target_fd(lat(-1)) # Test for warning with an empty template template = builder.Builder(kwant.TranslationalSymmetry((-1,))) target = builder.Builder() with warns(RuntimeWarning): target.fill(template, lambda x: True, lat(0)) # Test for warning when one of the starting sites is outside the template lat = kwant.lattice.square(norbs=1) template = builder.Builder(kwant.TranslationalSymmetry((-1, 0))) template[lat(0, 0)] = None template[lat.neighbors()] = None target = builder.Builder() with warns(RuntimeWarning): target.fill(template, lambda x: -1 < x.pos[0] < 1, [lat(0, 0), lat(0, 1)]) def test_fill_sticky(): """Test that adjacent regions are properly interconnected when filled separately. """ # Generate model template. lat = kwant.lattice.kagome(norbs=1) template = kwant.Builder(kwant.TranslationalSymmetry( lat.vec((1, 0)), lat.vec((0, 1)))) for i, sl in enumerate(lat.sublattices): template[sl(0, 0)] = i for i in range(1, 3): for j, hop in enumerate(template.expand(lat.neighbors(i))): template[hop] = j * 1j def disk(site): pos = site.pos return ta.dot(pos, pos) < 13 def halfplane(site): return ta.dot(site.pos - (-1, 1), (-0.9, 0.63)) > 0 # Fill in one go. syst0 = kwant.Builder() syst0.fill(template, disk, (0, 0)) # Fill in two stages. syst1 = kwant.Builder() syst1.fill(template, lambda s: disk(s) and halfplane(s), (-2, 1)) syst1.fill(template, lambda s: disk(s) and not halfplane(s), (0, 0)) # Verify that both results are identical. assert set(syst0.site_value_pairs()) == set(syst1.site_value_pairs()) assert (set(syst0.hopping_value_pairs()) == set(syst1.hopping_value_pairs())) def test_attach_lead(): fam = SimpleSiteFamily(norbs=1) fam_noncommensurate = SimpleSiteFamily(name='other', norbs=1) syst = builder.Builder() syst[fam(1)] = 0 lead = builder.Builder(VerySimpleSymmetry(-2)) raises(ValueError, syst.attach_lead, lead) lead[fam(0)] = 1 raises(ValueError, syst.attach_lead, lead) lead[fam(1)] = 1 syst.attach_lead(lead) raises(ValueError, syst.attach_lead, lead, fam(5)) syst = builder.Builder() # The tag of the site that is added in the following line is an empty tuple. # This simulates a site family that is not commensurate with the symmetry of # the lead. Such sites may be present in the system, as long as there are # other sites that will interrupt the lead. syst[fam_noncommensurate()] = 2 syst[fam(1)] = 0 syst[fam(0)] = 1 lead[fam(0), fam(1)] = lead[fam(0), fam(2)] = 1 syst.attach_lead(lead) assert len(list(syst.sites())) == 4 assert set(syst.leads[0].interface) == set([fam(-1), fam(0)]) syst[fam(-10)] = syst[fam(-11)] = 0 syst.attach_lead(lead) assert set(syst.leads[1].interface) == set([fam(-10), fam(-11)]) assert len(list(syst.sites())) == 6 syst.attach_lead(lead, fam(-5)) assert set(syst.leads[0].interface) == set([fam(-1), fam(0)]) # add some further-than-nearest-neighbor hoppings hop_range = 3 lead = builder.Builder( VerySimpleSymmetry(1), conservation_law=np.eye(1), time_reversal=np.eye(1), particle_hole=np.eye(1), chiral=np.eye(1)) lead[fam(0)] = 1 for i in range(1, hop_range + 1): lead[fam(0), fam(i)] = 1 syst.attach_lead(lead) expanded_lead = syst.leads[-1].builder assert expanded_lead.symmetry.period == hop_range assert len(list(expanded_lead.sites())) == hop_range assert expanded_lead.conservation_law is lead.conservation_law assert expanded_lead.time_reversal is lead.time_reversal assert expanded_lead.particle_hole is lead.particle_hole assert expanded_lead.chiral is lead.chiral # check that we can actually finalize the system syst.finalized() def test_attach_lead_incomplete_unit_cell(): lat = kwant.lattice.chain(norbs=1) syst = kwant.Builder() lead = kwant.Builder(kwant.TranslationalSymmetry((2,))) syst[lat(1)] = lead[lat(0)] = lead[lat(1)] = 0 lead[lat.neighbors()] = 0 assert(len(syst.attach_lead(lead)) == 0) def test_neighbors_not_in_single_domain(): sr = builder.Builder() lead = builder.Builder(VerySimpleSymmetry(-1)) fam = SimpleSiteFamily(norbs=1) sr[(fam(x, y) for x in range(3) for y in range(3) if x >= y)] = 0 sr[builder.HoppingKind((1, 0), fam)] = 1 sr[builder.HoppingKind((0, 1), fam)] = 1 lead[(fam(0, y) for y in range(3))] = 0 lead[((fam(0, y), fam(1, y)) for y in range(3))] = 1 lead[((fam(0, y), fam(0, y + 1)) for y in range(2))] = 1 sr.leads.append(builder.BuilderLead(lead, [fam(i, i) for i in range(3)])) raises(ValueError, sr.finalized) def inside_disc(center, rr): def shape(site): d = site.pos - center dd = ta.dot(d, d) return dd <= rr return shape def test_closest(): rng = ensure_rng(10) for sym_dim in range(1, 4): for space_dim in range(sym_dim, 4): lat = kwant.lattice.general(ta.identity(space_dim), norbs=1) # Choose random periods. while True: periods = rng.randint(-10, 11, (sym_dim, space_dim)) if np.linalg.det(np.dot(periods, periods.T)) > 0.1: # Periods are reasonably linearly independent. break syst = builder.Builder(kwant.TranslationalSymmetry(*periods)) for tag in rng.randint(-30, 31, (4, space_dim)): # Add site and connect it to the others. old_sites = list(syst.sites()) new_site = lat(*tag) syst[new_site] = None syst[((new_site, os) for os in old_sites)] = None # Test consistency with fill(). for point in 200 * rng.random_sample((10, space_dim)) - 100: closest = syst.closest(point) dist = closest.pos - point dist = ta.dot(dist, dist) syst2 = builder.Builder() syst2.fill(syst, inside_disc(point, 2 * dist), closest) assert syst2.closest(point) == closest for site in syst2.sites(): dd = site.pos - point dd = ta.dot(dd, dd) assert dd >= 0.999999 * dist def test_update(): lat = SimpleSiteFamily(norbs=1) syst = builder.Builder() syst[[lat(0,), lat(1,)]] = 1 syst[lat(0,), lat(1,)] = 1 other_syst = builder.Builder() other_syst[[lat(1,), lat(2,)]] = 2 other_syst[lat(1,), lat(2,)] = 1 lead0 = builder.Builder(VerySimpleSymmetry(-1)) lead0[lat(0,)] = 1 lead0[(lat(0,), lat(1,))] = 1 lead0 = builder.BuilderLead(lead0, [lat(0,)]) syst.leads.append(lead0) lead1 = builder.Builder(VerySimpleSymmetry(1)) lead1[lat(2,)] = 1 lead1[(lat(2,), lat(1,))] = 1 lead1 = builder.BuilderLead(lead1, [lat(2,)]) other_syst.leads.append(lead1) syst.update(other_syst) assert syst.leads == [lead0, lead1] expected = sorted([((0,), 1), ((1,), 2), ((2,), 2)]) assert sorted(((s.tag, v) for s, v in syst.site_value_pairs())) == expected expected = sorted([((0,), (1,), 1), ((1,), (2,), 1)]) assert (sorted(((a.tag, b.tag, v) for (a, b), v in syst.hopping_value_pairs())) == expected) # y=0: y=1: # # hhhh hhhh # gggh gggh # hhgh hhgh # ghgh hhgh # def test_HoppingKind(): g = kwant.lattice.general(ta.identity(3), name='some_lattice', norbs=1) h = kwant.lattice.general(ta.identity(3), name='another_lattice', norbs=1) sym = kwant.TranslationalSymmetry((0, 2, 0)) syst = builder.Builder(sym) syst[((h if max(x, y, z) % 2 else g)(x, y, z) for x in range(4) for y in range(2) for z in range(4))] = None for delta, ga, gb, n in [((1, 0, 0), g, h, 4), ((1, 0, 0), h, g, 7), ((0, 1, 0), g, h, 1), ((0, 4, 0), h, h, 21), ((0, 0, 1), g, h, 4)]: ph = list(builder.HoppingKind(delta, ga, gb)(syst)) assert len(ph) == n ph = set(ph) assert len(ph) == n ph2 = list(( sym.to_fd(b, a) for a, b in builder.HoppingKind(ta.negative(delta), gb, ga)(syst))) assert len(ph2) == n ph2 = set(ph2) assert ph2 == ph for a, b in ph: assert a.family == ga assert b.family == gb assert sym.to_fd(a) == a assert a.tag - b.tag == delta # test hashability and equality hk = builder.HoppingKind((1, 0, 0), g) hk2 = builder.HoppingKind((1, 0, 0), g) hk3 = builder.HoppingKind((1, 0, 0), g, h) assert hk == hk2 assert hash(hk) == hash(hk2) assert hk != hk3 assert hash(hk) != hash(hk3) assert len({hk: 0, hk2:1, hk3: 2}) == 2 def test_invalid_HoppingKind(): g = kwant.lattice.general(ta.identity(3), norbs=1) h = kwant.lattice.general(np.identity(3)[:-1], norbs=1) # 2D lattice in 3D delta = (1, 0, 0) # families have incompatible tags with raises(ValueError): builder.HoppingKind(delta, g, h) # delta is incompatible with tags with raises(ValueError): builder.HoppingKind(delta, h) def test_ModesLead_and_SelfEnergyLead(): lat = SimpleSiteFamily(norbs=1) hoppings = [builder.HoppingKind((1, 0), lat), builder.HoppingKind((0, 1), lat)] rng = Random(123) L = 5 t = 1 energies = [0.9, 1.7] syst = builder.Builder() for x in range(L): for y in range(L): syst[lat(x, y)] = 4 * t + rng.random() - 0.5 syst[hoppings] = -t # Attach a lead from the left. lead = builder.Builder(VerySimpleSymmetry(-1)) for y in range(L): lead[lat(0, y)] = 4 * t lead[hoppings] = -t syst.attach_lead(lead) # Make the right lead and attach it. lead = builder.Builder(VerySimpleSymmetry(1)) for y in range(L): lead[lat(0, y)] = 4 * t lead[hoppings] = -t syst.attach_lead(lead) fsyst = syst.finalized() ts = [kwant.smatrix(fsyst, e).transmission(1, 0) for e in energies] # Replace lead with it's finalized copy. lead = fsyst.leads[1] interface = [lat(L-1, lead.sites[i].tag[1]) for i in range(L)] # Re-attach right lead as ModesLead. syst.leads[1] = builder.ModesLead(lead.modes, interface, lead.parameters) fsyst = syst.finalized() ts2 = [kwant.smatrix(fsyst, e).transmission(1, 0) for e in energies] assert_almost_equal(ts2, ts) # Re-attach right lead as ModesLead with old-style modes API # that does not take a 'params' keyword parameter. syst.leads[1] = builder.ModesLead( lambda energy, args: lead.modes(energy, args), interface, lead.parameters) fsyst = syst.finalized() ts2 = [kwant.smatrix(fsyst, e).transmission(1, 0) for e in energies] assert_almost_equal(ts2, ts) # Re-attach right lead as SelfEnergyLead. syst.leads[1] = builder.SelfEnergyLead(lead.selfenergy, interface, lead.parameters) fsyst = syst.finalized() ts2 = [kwant.greens_function(fsyst, e).transmission(1, 0) for e in energies] assert_almost_equal(ts2, ts) # Re-attach right lead as SelfEnergyLead with old-style selfenergy API # that does not take a 'params' keyword parameter. syst.leads[1] = builder.SelfEnergyLead( lambda energy, args: lead.selfenergy(energy, args), interface, lead.parameters) fsyst = syst.finalized() ts2 = [kwant.greens_function(fsyst, e).transmission(1, 0) for e in energies] assert_almost_equal(ts2, ts) # Append a virtual (=zero self energy) lead. This should have no effect. # Also verifies that the selfenergy callback function can return exotic # arraylikes. syst.leads.append(builder.SelfEnergyLead( lambda *args: list(ta.zeros((L, L))), interface, lead.parameters)) fsyst = syst.finalized() ts2 = [kwant.greens_function(fsyst, e).transmission(1, 0) for e in energies] assert_almost_equal(ts2, ts) def test_site_pickle(): site = kwant.lattice.square(norbs=1)(0, 0) assert pickle.loads(pickle.dumps(site)) == site def test_discrete_symmetries(): lat = SimpleSiteFamily(name='ccc', norbs=2) lat2 = SimpleSiteFamily(name='bla', norbs=1) lat3 = SimpleSiteFamily(name='dd', norbs=4) cons_law = {lat: np.diag([1, 2]), lat2: 0} syst = builder.Builder(conservation_law=cons_law, time_reversal=(lambda site, p: np.exp(1j*p) * np.identity(site.family.norbs))) syst[lat(1)] = np.identity(2) syst[lat2(1)] = 1 params=dict(p=0) sym = syst.finalized().discrete_symmetry(params=params) for proj, should_be in zip(sym.projectors, np.identity(3)): assert np.allclose(proj.toarray(), should_be.reshape((3, 1))) assert np.allclose(sym.time_reversal.toarray(), np.identity(3)) syst.conservation_law = lambda site, p: cons_law[site.family] sym = syst.finalized().discrete_symmetry(params=params) for proj, should_be in zip(sym.projectors, np.identity(3)): assert np.allclose(proj.toarray(), should_be.reshape((-1, 1))) syst = builder.Builder(conservation_law=np.diag([-1, 1])) syst[lat(1)] = np.identity(2) sym = syst.finalized().discrete_symmetry() for proj, should_be in zip(sym.projectors, np.identity(2)): assert np.allclose(proj.toarray(), should_be.reshape((-1, 1))) syst = builder.Builder(conservation_law=1) syst[lat2(1)] = 0 sym = syst.finalized().discrete_symmetry() [proj] = sym.projectors assert np.allclose(proj.toarray(), [[1]]) syst = kwant.Builder(conservation_law=np.diag([-1, 1, -1, 1])) syst[lat3(0)] = np.eye(4) sym = syst.finalized().discrete_symmetry() p1 = np.zeros((4, 2)) p1[0, 0] = p1[2, 1] = 1 assert np.allclose(sym.projectors[0].toarray(), p1) p2 = np.zeros((4, 2)) p2[1, 0] = p2[3, 1] = 1 assert np.allclose(sym.projectors[1].toarray(), p2) # test parameter passing to conservation_law syst = builder.Builder(conservation_law=lambda site, b: b) syst[lat2(1)] = 0 sym = syst.finalized().discrete_symmetry(params=dict(a=None, b=1)) [proj] = sym.projectors assert np.allclose(proj.toarray(), [[1]]) # We need to keep testing 'args', but we don't want to see # all the deprecation warnings in the test logs @pytest.mark.filterwarnings("ignore:.*'args' parameter") @pytest.mark.parametrize("vectorize", [False, True]) def test_argument_passing(vectorize): chain = kwant.lattice.chain(norbs=1) # Test for passing parameters to hamiltonian matrix elements def onsite(site, p1, p2): if vectorize: return (p1 + p2) * np.ones(len(site)) else: return p1 + p2 def hopping(site1, site2, p1, p2): if vectorize: return (p1 - p2) * np.ones(len(site1)) else: return p1 - p2 def gen_fill_syst(onsite, hopping, syst): syst[(chain(i) for i in range(3))] = onsite syst[chain.neighbors()] = hopping return syst.finalized() fill_syst = ft.partial(gen_fill_syst, onsite, hopping) syst = fill_syst(kwant.Builder(vectorize=vectorize)) inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,)), vectorize=vectorize)) tests = ( syst.hamiltonian_submatrix, inf_syst.cell_hamiltonian, inf_syst.inter_cell_hopping, inf_syst.selfenergy, lambda *args, **kw: inf_syst.modes(*args, **kw)[0].wave_functions, ) for test in tests: np.testing.assert_array_equal( test(args=(2, 1)), test(params=dict(p1=2, p2=1))) # test that mixing 'args' and 'params' raises TypeError with raises(TypeError): with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1)) with raises(TypeError): with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) inf_syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1)) # Missing parameters raises TypeError with raises(TypeError): with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) syst.hamiltonian(0, 0, params=dict(p1=2)) with raises(TypeError): with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=KwantDeprecationWarning) syst.hamiltonian_submatrix(params=dict(p1=2)) # test that passing parameters without default values works, and that # passing parameters with default values fails fill_syst = ft.partial(gen_fill_syst, onsite, hopping) syst = fill_syst(kwant.Builder(vectorize=vectorize)) inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,)), vectorize=vectorize)) tests = ( syst.hamiltonian_submatrix, inf_syst.cell_hamiltonian, inf_syst.inter_cell_hopping, inf_syst.selfenergy, lambda *args, **kw: inf_syst.modes(*args, **kw)[0].wave_functions, ) for test in tests: np.testing.assert_array_equal( test(args=(1, 2)), test(params=dict(p1=1, p2=2))) # Some common, some different args for value functions def onsite2(site, a, b): if vectorize: return site.positions()[:, 0] + a + b else: return site.pos[0] + a + b def hopping2(site1, site2, a, c, b): if vectorize: return (a + b + c) * np.ones(len(site1)) else: return a + b + c syst = kwant.Builder(vectorize=vectorize) syst[(chain(i) for i in range(3))] = onsite2 syst[((chain(i), chain(i + 1)) for i in range(2))] = hopping2 fsyst = syst.finalized() def expected_hamiltonian(a, b, c): return [[a + b, a + b + c, 0], [a + b + c, 1 + a + b, a + b + c], [0, a+ b + c, 2 + a + b]] params = dict(a=1, b=2, c=3) np.testing.assert_array_equal( fsyst.hamiltonian_submatrix(params=params), expected_hamiltonian(**params) ) @pytest.mark.parametrize("vectorize", [False, True]) def test_invalid_value_functions(vectorize): invalid_value_functions = [ lambda _, *args: 1, # uses *args lambda _, **kwargs: 1, # uses **kwargs lambda _, d=10: 1, # has default arguments lambda _, *, d: 1 # has keyword-only parameters ] lat = kwant.lattice.chain(norbs=1) for f in invalid_value_functions: syst = builder.Builder(vectorize=vectorize) syst[lat(0)] = f syst = syst.finalized() with pytest.raises(ValueError): syst.hamiltonian_submatrix(params=dict(d=1)) def test_parameter_substitution(): subs = builder._substitute_params def f(x, y): return (('x', x), ('y', y)) # 'f' already has a parameter 'y' assert raises(ValueError, subs, f, dict(x='y')) # Swap argument names. g = subs(f, dict(x='y', y='x')) assert g(1, 2) == f(1, 2) # Swap again. h = subs(g, dict(x='y', y='x')) assert h(1, 2) == f(1, 2) # don't nest wrappers inside each other assert h.func is f # Try different names. g = subs(f, dict(x='a', y='b')) assert g(1, 2) == f(1, 2) # Can substitutions be used in sets/dicts? g = subs(f, dict(x='a')) h = subs(f, dict(x='a')) assert len(set([f, g, h])) == 2 def test_subs(): # Simple case def onsite(site, a, b): salt = str(a) + str(b) return kwant.digest.uniform(site.tag, salt=salt) def hopping(sitea, siteb, b, c): salt = str(b) + str(c) return kwant.digest.uniform(ta.array((sitea.tag, siteb.tag)), salt=salt) lat = kwant.lattice.chain(norbs=1) def make_system(sym=kwant.builder.NoSymmetry(), n=3): syst = kwant.Builder(sym) syst[(lat(i) for i in range(n))] = onsite syst[lat.neighbors()] = hopping return syst def hamiltonian(syst, **kwargs): return syst.finalized().hamiltonian_submatrix(params=kwargs) syst = make_system() # substituting a paramter that doesn't exist produces a warning warns(RuntimeWarning, syst.substituted, fakeparam='yes') # name clash in value functions raises(ValueError, syst.substituted, b='a') raises(ValueError, syst.substituted, b='c') raises(ValueError, syst.substituted, a='site') raises(ValueError, syst.substituted, c='sitea') # cannot call 'substituted' on systems with attached leads, because # it is not clear whether the substitutions should propagate # into the leads too. syst = make_system() lead = make_system(kwant.TranslationalSymmetry((-1,)), n=1) syst.attach_lead(lead) raises(ValueError, syst.substituted, a='d') # test basic substitutions syst = make_system() assert syst.finalized().parameters == {'a', 'b', 'c'} expected = hamiltonian(syst, a=1, b=2, c=3) # 1 level of substitutions sub_syst = syst.substituted(a='d', b='e') assert sub_syst.finalized().parameters == {'d', 'e', 'c'} assert np.allclose(hamiltonian(sub_syst, d=1, e=2, c=3), expected) # 2 levels of substitution sub_sub_syst = sub_syst.substituted(d='g', c='h') assert np.allclose(hamiltonian(sub_sub_syst, g=1, e=2, h=3), expected) assert sub_sub_syst.finalized().parameters == {'g', 'e', 'h'} # very confusing but technically valid. 'a' does not appear in 'hopping', # so the signature of 'onsite' is valid. sub_syst = syst.substituted(a='sitea') assert sub_syst.finalized().parameters == {'sitea', 'b', 'c'} assert np.allclose(hamiltonian(sub_syst, sitea=1, b=2, c=3), expected) # Check that this also works for infinite systems, as their finalization # follows a different code path. lead = make_system(kwant.TranslationalSymmetry((-1,)), n=1) lead = lead.substituted(a='lead_a', b='lead_b', c='lead_c') lead = lead.finalized() assert lead.parameters == {'lead_a', 'lead_b', 'lead_c'} def test_attach_stores_padding(): lat = kwant.lattice.chain(norbs=1) syst = kwant.Builder() syst[lat(0)] = 0 lead = kwant.Builder(kwant.TranslationalSymmetry(lat.prim_vecs[0])) lead[lat(0)] = 0 lead[lat(1), lat(0)] = 0 added_sites = syst.attach_lead(lead, add_cells=2) assert syst.leads[0].padding == sorted(added_sites) def test_finalization_preserves_padding(): lat = kwant.lattice.chain(norbs=1) syst = kwant.Builder() for i in range(10): syst[lat(i)] = 0 lead = kwant.Builder(kwant.TranslationalSymmetry(lat.prim_vecs[0])) lead[lat(0)] = 0 lead[lat(0), lat(1)] = 0 # We use a low level way to provide a lead to directly check that the # padding is preserved. We also check that the sites that do not belong to # the system don't end up in the padding of the finalized system. padding = [lat(0), lat(3), lat(5), lat(11)] syst.leads.append(kwant.builder.BuilderLead(lead, [lat(0)], padding)) syst = syst.finalized() # The order is guaranteed because the paddings are sorted. assert [syst.sites[i] for i in syst.lead_paddings[0]] == padding[:-1] def test_add_peierls_phase(): lat = kwant.lattice.square(norbs=1) syst = kwant.Builder() syst[(lat(i, j) for i in range(5) for j in range(5))] = 4 syst[lat.neighbors()] = lambda a, b, t: -t lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0))) lead[(lat(0, j) for j in range(5))] = 4 lead[lat.neighbors()] = -1 # System without leads syst_without_lead, phase_without_lead = builder.add_peierls_phase(syst) assert isinstance(syst_without_lead, builder.FiniteSystem) params_without_lead = phase_without_lead(1) assert all(p in params_without_lead for p in ('phi',)) # System with leads syst.attach_lead(lead) syst.attach_lead(lead.reversed()) syst, phase = builder.add_peierls_phase(syst) assert isinstance(syst, builder.FiniteSystem) params = phase(1, 0, 0) assert all(p in params for p in ('phi', 'phi_lead0', 'phi_lead1')) kwant.smatrix(syst, energy=0.1, params=dict(t=-1, **params))