add kwant.qsymm module

Co-authored-by: Daniel Varjas <dvarjas@gmail.com>

kwant/qsymm.py

+# Copyright 2011-2018 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.
+
+
+__all__ = ['builder_to_model', 'model_to_builder', 'find_builder_symmetries']
+
+import itertools as it
+from collections import OrderedDict, defaultdict
+
+import numpy as np
+import tinyarray as ta
+import scipy.linalg as la
+
+try:
+    import sympy
+    import qsymm
+    from qsymm.model import Model, BlochModel, BlochCoeff
+    from qsymm.groups import PointGroupElement, ContinuousGroupGenerator
+    from qsymm.symmetry_finder import bravais_point_group
+    from qsymm.linalg import allclose
+    from qsymm.hamiltonian_generator import hamiltonian_from_family
+except ImportError as error:
+    msg = ("'kwant.qsymm' is not available because one or more of its "
+           "dependencies is not installed.")
+    raise ImportError(msg) from error
+
+from kwant import lattice, builder
+from kwant._common import get_parameters
+
+
+def builder_to_model(syst, momenta=None, real_space=True,
+                     params=None):
+    """Make a qsymm.BlochModel out of a Builder.
+
+    Parameters
+    ----------
+    syst: kwant.Builder
+        May have translational symmetries.
+    momenta: list of strings or None
+        Names of momentum variables. If None, 'k_x', 'k_y', ... is used.
+    real_space: bool (default True)
+        If False, use the unit cell convention for Bloch basis, the
+        exponential has the difference in the unit cell coordinates and
+        k is expressed in the reciprocal lattice basis. This is consistent
+        with kwant.wraparound.
+        If True, the difference in the real space coordinates is used
+        and k is given in an absolute basis.
+        Only the default choice guarantees that qsymm is able to find
+        nonsymmorphic symmetries.
+    params: dict, optional
+        Dictionary of parameter names and their values; used when
+        evaluating the Hamiltonian matrix elements.
+
+    Returns
+    -------
+    qsymm.BlochModel
+        Model representing the tight-binding Hamiltonian.
+
+    Notes
+    -----
+    The sites in the the builder are in lexicographical order, i.e. ordered
+    first by their family and then by their tag. This is the same ordering that
+    is used in finalized kwant systems.
+    """
+    def term_to_model(d, par, matrix):
+        if np.allclose(matrix, 0):
+            result = BlochModel({})
+        else:
+            result = BlochModel({BlochCoeff(d, qsymm.sympify(par)): matrix},
+                                momenta=momenta)
+        return result
+
+    def hopping_to_model(hop, value, proj, params):
+        site1, site2 = hop
+        if real_space:
+            d = proj @ np.array(site2.pos - site1.pos)
+        else:
+            # site in the FD
+            d = np.array(syst.symmetry.which(site2))
+
+        slice1, slice2 = slices[to_fd(site1)], slices[to_fd(site2)]
+        if callable(value):
+            return sum(term_to_model(d, par, set_block(slice1, slice2, val))
+                       for par, val in function_to_terms(hop, value, params))
+        else:
+            matrix = set_block(slice1, slice2, value)
+            return term_to_model(d, '1', matrix)
+
+    def onsite_to_model(site, value, params):
+        d = np.zeros((dim, ))
+        slice1 = slices[to_fd(site)]
+        if callable(value):
+            return sum(term_to_model(d, par, set_block(slice1, slice1, val))
+                       for par, val in function_to_terms(site, value, params))
+        else:
+            return term_to_model(d, '1', set_block(slice1, slice1, value))
+
+    def function_to_terms(site_or_hop, value, fixed_params):
+        assert callable(value)
+        parameters = get_parameters(value)
+        # remove site or site1, site2 parameters
+        if isinstance(site_or_hop, builder.Site):
+            parameters = parameters[1:]
+            site_or_hop = (site_or_hop,)
+        else:
+            parameters = parameters[2:]
+        free_parameters = (par for par in parameters
+                           if par not in fixed_params.keys())
+        # first set all free parameters to 0
+        args = ((fixed_params[par] if par in fixed_params.keys() else 0)
+                for par in parameters)
+        h_0 = value(*site_or_hop, *args)
+        # set one of the free parameters to 1 at a time, the rest 0
+        terms = []
+        for p in free_parameters:
+            args = ((fixed_params[par] if par in fixed_params.keys() else
+                     (1 if par == p else 0)) for par in parameters)
+            terms.append((p, value(*site_or_hop, *args) - h_0))
+        return terms + [('1', h_0)]
+
+    def orbital_slices(syst):
+        orbital_slices = {}
+        start_orb = 0
+
+        for site in sorted(syst.sites()):
+            n = site.family.norbs
+            if n is None:
+                raise ValueError('norbs must be provided for every lattice.')
+            orbital_slices[site] = slice(start_orb, start_orb + n)
+            start_orb += n
+        return orbital_slices, start_orb
+
+    def set_block(slice1, slice2, val):
+        matrix = np.zeros((N, N), dtype=complex)
+        matrix[slice1, slice2] = val
+        return matrix
+
+    if params is None:
+        params = dict()
+
+    periods = np.array(syst.symmetry.periods)
+    dim = len(periods)
+    to_fd = syst.symmetry.to_fd
+    if momenta is None:
+        momenta = ['k_x', 'k_y', 'k_z'][:dim]
+    # If the system is higher dimensional than the number of translation
+    # vectors, we need to project onto the subspace spanned by the
+    # translation vectors.
+    if dim == 0:
+        proj = np.empty((0, len(list(syst.sites())[0].pos)))
+    elif dim < len(list(syst.sites())[0].pos):
+        proj, r = la.qr(np.array(periods).T, mode='economic')
+        sign = np.diag(np.diag(np.sign(r)))
+        proj = sign @ proj.T
+    else:
+        proj = np.eye(dim)
+
+    slices, N = orbital_slices(syst)
+
+    one_way_hoppings = [hopping_to_model(hop, value, proj, params)
+                        for hop, value in syst.hopping_value_pairs()]
+    other_way_hoppings = [term.T().conj() for term in one_way_hoppings]
+    hoppings = one_way_hoppings + other_way_hoppings
+
+    onsites = [onsite_to_model(site, value, params)
+               for site, value in syst.site_value_pairs()]
+
+    result = sum(onsites) + sum(hoppings)
+
+    return result
+
+
+def model_to_builder(model, norbs, lat_vecs, atom_coords, *, coeffs=None):
+    """Make a Builder out of qsymm.Models or qsymm.BlochModels.
+
+    Parameters
+    ----------
+    model: qsymm.model.Model, qsymm.model.BlochModel, or an iterable thereof
+        The Hamiltonian (or terms of the Hamiltonian) to convert to a
+        Builder.
+    norbs: OrderedDict or sequence of pairs
+        Maps sites to the number of orbitals per site in a unit cell.
+    lat_vecs: list of arrays
+        Lattice vectors of the underlying tight binding lattice.
+    atom_coords: list of arrays
+        Positions of the sites (or atoms) within a unit cell.
+        The ordering of the atoms is the same as in norbs.
+    coeffs: list of sympy.Symbol, default None.
+        Constant prefactors for the individual terms in model, if model
+        is a list of multiple objects. If model is a single Model or BlochModel
+        object, this argument is ignored. By default assigns the coefficient
+        c_n to element model[n].
+
+    Returns
+    -------
+    syst: kwant.Builder
+        The unfinalized Kwant system representing the qsymm Model(s).
+
+    Notes
+    -----
+    Onsite terms that are not provided in the input model are set
+    to zero by default.
+    The input Model(s) representing the tight binding Hamiltonian in
+    Bloch form should follow the convention where the difference in the real
+    space atomic positions appear in the Bloch factors.
+    """
+
+    def make_int(R):
+        # If close to an integer array convert to integer tinyarray, else
+        # return None
+        R_int = ta.array(np.round(R), int)
+        if qsymm.linalg.allclose(R, R_int):
+            return R_int
+        else:
+            return None
+
+    def term_onsite(onsites_dict, hopping_dict, hop_mat, atoms,
+                    sublattices, coords_dict):
+        """Find the Kwant onsites and hoppings in a qsymm.BlochModel term
+        that has no lattice translation in the Bloch factor.
+        """
+        for atom1, atom2 in it.product(atoms, atoms):
+            # Subblock within the same sublattice is onsite
+            hop = hop_mat[ranges[atom1], ranges[atom2]]
+            if sublattices[atom1] == sublattices[atom2]:
+                onsites_dict[atom1] += Model({coeff: hop}, momenta=momenta)
+            # Blocks between sublattices are hoppings between sublattices
+            # at the same position.
+            # Only include nonzero hoppings
+            elif not allclose(hop, 0):
+                if not allclose(np.array(coords_dict[atom1]),
+                                np.array(coords_dict[atom2])):
+                    raise ValueError(
+                        "Position of sites not compatible with qsymm model.")
+                lat_basis = np.array(zer)
+                hop = Model({coeff: hop}, momenta=momenta)
+                hop_dir = builder.HoppingKind(-lat_basis, sublattices[atom1],
+                                              sublattices[atom2])
+                hopping_dict[hop_dir] += hop
+        return onsites_dict, hopping_dict
+
+    def term_hopping(hopping_dict, hop_mat, atoms,
+                     sublattices, coords_dict):
+        """Find Kwant hoppings in a qsymm.BlochModel term that has a lattice
+        translation in the Bloch factor.
+        """
+        # Iterate over combinations of atoms, set hoppings between each
+        for atom1, atom2 in it.product(atoms, atoms):
+            # Take the block from atom1 to atom2
+            hop = hop_mat[ranges[atom1], ranges[atom2]]
+            # Only include nonzero hoppings
+            if allclose(hop, 0):
+                continue
+            # Adjust hopping vector to Bloch form basis
+            r_lattice = (
+                r_vec
+                + np.array(coords_dict[atom1])
+                - np.array(coords_dict[atom2])
+            )
+            # Bring vector to basis of lattice vectors
+            lat_basis = np.linalg.solve(np.vstack(lat_vecs).T, r_lattice)
+            lat_basis = make_int(lat_basis)
+            # Should only have hoppings that are integer multiples of
+            # lattice vectors
+            if lat_basis is not None:
+                hop_dir = builder.HoppingKind(-lat_basis,
+                                              sublattices[atom1],
+                                              sublattices[atom2])
+                # Set the hopping as the matrix times the hopping amplitude
+                hopping_dict[hop_dir] += Model({coeff: hop}, momenta=momenta)
+            else:
+                raise RuntimeError('A nonzero hopping not matching a '
+                                   'lattice vector was found.')
+        return hopping_dict
+
+    # Disambiguate single model instances from iterables thereof. Because
+    # Model is itself iterable (subclasses dict) this is a bit cumbersome.
+    if isinstance(model, Model):
+        # BlochModel can't yet handle getting a Blochmodel as input
+        if not isinstance(model, BlochModel):
+            model = BlochModel(model)
+    else:
+        model = BlochModel(hamiltonian_from_family(
+            model, coeffs=coeffs, nsimplify=False, tosympy=False))
+
+
+    # 'momentum' and 'zer' are used in the closures defined above, so don't
+    # move these declarations down.
+    momenta = model.momenta
+    if len(momenta) != len(lat_vecs):
+        raise ValueError("Dimension of the lattice and number of "
+                         "momenta do not match.")
+    zer = [0] * len(momenta)
+
+
+    # Subblocks of the Hamiltonian for different atoms.
+    N = 0
+    if not any([isinstance(norbs, OrderedDict), isinstance(norbs, list),
+                isinstance(norbs, tuple)]):
+        raise ValueError('norbs must be OrderedDict, tuple, or list.')
+    else:
+        norbs = OrderedDict(norbs)
+    ranges = dict()
+    for a, n in norbs.items():
+        ranges[a] = slice(N, N + n)
+        N += n
+
+    # Extract atoms and number of orbitals per atom,
+    # store the position of each atom
+    atoms, orbs = zip(*norbs.items())
+    coords_dict = dict(zip(atoms, atom_coords))
+
+    # Make the kwant lattice
+    lat = lattice.general(lat_vecs, atom_coords, norbs=orbs)
+    # Store sublattices by name
+    sublattices = dict(zip(atoms, lat.sublattices))
+
+    # Keep track of the hoppings and onsites by storing those
+    # which have already been set.
+    hopping_dict = defaultdict(dict)
+    onsites_dict = defaultdict(dict)
+
+    # Iterate over all terms in the model.
+    for key, hop_mat in model.items():
+        # Determine whether this term is an onsite or a hopping, extract
+        # overall symbolic coefficient if any, extract the exponential
+        # part describing the hopping if present.
+        r_vec, coeff = key
+        # Onsite term; modifies onsites_dict and hopping_dict in-place
+        if allclose(r_vec, 0):
+            term_onsite(
+                onsites_dict, hopping_dict, hop_mat,
+                atoms, sublattices, coords_dict)
+        # Hopping term; modifies hopping_dict in-place
+        else:
+            term_hopping(hopping_dict, hop_mat, atoms,
+                         sublattices, coords_dict)
+
+    # If some onsite terms are not set, we set them to zero.
+    for atom in atoms:
+        if atom not in onsites_dict:
+            onsites_dict[atom] = Model(
+                {sympy.numbers.One(): np.zeros((norbs[atom], norbs[atom]))},
+                momenta=momenta)
+
+    # Make the Kwant system, and set all onsites and hoppings.
+
+    sym = lattice.TranslationalSymmetry(*lat_vecs)
+    syst = builder.Builder(sym)
+
+    # Iterate over all onsites and set them
+    for atom, onsite in onsites_dict.items():
+        syst[sublattices[atom](*zer)] = onsite.lambdify(onsite=True)
+
+    # Finally, iterate over all the hoppings and set them
+    for direction, hopping in hopping_dict.items():
+        syst[direction] = hopping.lambdify(hopping=True)
+
+    return syst
+
+
+# This may be useful in the future, so we'll keep it as internal for now,
+# and can make it part of the API in the future if we wish.
+def _get_builder_symmetries(builder):
+    """Extract the declared symmetries of a Kwant builder.
+
+    Parameters
+    ----------
+    builder: kwant.Builder
+
+    Returns
+    -------
+    builder_symmetries: dict
+        Dictionary of the discrete symmetries that the builder has.
+        The symmetries can be particle-hole, time-reversal or chiral,
+        which are returned as qsymm.PointGroupElements, or
+        a conservation law, which is returned as a
+        qsymm.ContinuousGroupGenerators.
+    """
+
+    dim = len(np.array(builder.symmetry.periods))
+    symmetry_names = ['time_reversal', 'particle_hole', 'chiral',
+                      'conservation_law']
+    builder_symmetries = {name: getattr(builder, name)
+                          for name in symmetry_names
+                          if getattr(builder, name) is not None}
+    for name, symmetry in builder_symmetries.items():
+        if name == 'time_reversal':
+            builder_symmetries[name] = PointGroupElement(np.eye(dim),
+                                                         True, False, symmetry)
+        elif name == 'particle_hole':
+            builder_symmetries[name] = PointGroupElement(np.eye(dim),
+                                                         True, True, symmetry)
+        elif name == 'chiral':
+            builder_symmetries[name] = PointGroupElement(np.eye(dim),
+                                                         False, True, symmetry)
+        elif name == 'conservation_law':
+            builder_symmetries[name] = ContinuousGroupGenerator(R=None,
+                                                                U=symmetry)
+        else:
+            raise ValueError("Invalid symmetry name.")
+    return builder_symmetries
+
+
+def find_builder_symmetries(builder, momenta=None, params=None,
+                            spatial_symmetries=True, prettify=True):
+    """Finds the symmetries of a Kwant system using qsymm.
+
+    Parameters
+    ----------
+    builder: kwant.Builder
+    momenta: list of strings or None
+        Names of momentum variables, if None 'k_x', 'k_y', ... is used.
+    params: dict, optional
+        Dictionary of parameter names and their values; used when
+        evaluating the Hamiltonian matrix elements.
+    spatial_symmetries: bool (default True)
+        If True, search for all symmetries.
+        If False, only searches for the symmetries that are declarable in
+        kwant.Builder objects, i.e. time-reversal symmetry, particle-hole
+        symmetry, chiral symmetry, or conservation laws. This can save
+        computation time.
+    prettify : bool (default True)
+        Whether to carry out sparsification of the continuous symmetry
+        generators, in general an arbitrary linear combination of the
+        symmetry generators is returned.
+
+    Returns
+    -------
+    symmetries: list of qsymm.PointGroupElements and/or
+                qsymm.ContinuousGroupElement
+        The symmetries of the Kwant system.
+    """
+
+    if params is None:
+        params = dict()
+
+    ham = builder_to_model(builder, momenta=momenta,
+                           real_space=False, params=params)
+
+    dim = len(np.array(builder.symmetry.periods))
+
+    if spatial_symmetries:
+        candidates = bravais_point_group(builder.symmetry.periods, tr=True,
+                                         ph=True, generators=False,
+                                         verbose=False)
+    else:
+        candidates = [
+            qsymm.PointGroupElement(np.eye(dim), True, False, None),  # T
+            qsymm.PointGroupElement(np.eye(dim), True, True, None),   # P
+            qsymm.PointGroupElement(np.eye(dim), False, True, None)]  # C
+    sg, cg = qsymm.symmetries(ham, candidates, prettify=prettify,
+                              continuous_rotations=False)
+    return list(sg) + list(cg)

kwant/tests/test_qsymm.py

+# Copyright 2011-2018 Kwant authors.
+#
+# This file is part of Kwant.  It is subject to the license terms in the file
+# LICENSE.rst found in the top-level directory of this distribution and at
+# http://kwant-project.org/license.  A list of Kwant authors can be found in
+# the file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from collections import OrderedDict
+
+import numpy as np
+import sympy
+
+import qsymm
+from qsymm.symmetry_finder import symmetries
+from qsymm.hamiltonian_generator import bloch_family, hamiltonian_from_family
+from qsymm.groups import (hexagonal, PointGroupElement, spin_matrices,
+                          spin_rotation, ContinuousGroupGenerator)
+from qsymm.model import Model, e, I, _commutative_momenta
+from qsymm.linalg import allclose
+
+import kwant
+from kwant._common import ensure_rng
+from kwant.lattice import TranslationalSymmetry
+from kwant.qsymm import (builder_to_model, model_to_builder,
+                         _get_builder_symmetries, find_builder_symmetries)
+
+
+def test_honeycomb():
+    lat = kwant.lattice.honeycomb(norbs=1)
+
+    # Test simple honeycomb model with constant terms
+    # Add discrete symmetries to the kwant builder as well, to check that they are
+    # returned as well.
+    syst = kwant.Builder(symmetry=TranslationalSymmetry(*lat.prim_vecs))
+    syst[lat.a(0, 0)] = 1
+    syst[lat.b(0, 0)] = 1
+    syst[lat.neighbors(1)] = -1
+
+    H = builder_to_model(syst)
+    sg, cs = symmetries(H, hexagonal(sympy_R=False), prettify=True)
+    assert len(sg) == 24
+    assert len(cs) == 0
+
+    # Test simple honeycomb model with value functions
+    syst = kwant.Builder(symmetry=TranslationalSymmetry(*lat.prim_vecs))
+    syst[lat.a(0, 0)] = lambda site, ma: ma
+    syst[lat.b(0, 0)] = lambda site, mb: mb
+    syst[lat.neighbors(1)] = lambda site1, site2, t: t
+
+    H = builder_to_model(syst)
+    sg, cs = symmetries(H, hexagonal(sympy_R=False), prettify=True)
+    assert len(sg) == 12
+    assert len(cs) == 0
+
+
+def test_get_builder_symmetries():
+    syst = kwant.Builder(kwant.TranslationalSymmetry((1, 0), (0, 1)),
+                         particle_hole=np.eye(2),
+                         conservation_law=2*np.eye(2))
+    builder_symmetries = _get_builder_symmetries(syst)
+
+    assert len(builder_symmetries) == 2
+    P = builder_symmetries['particle_hole']
+    assert isinstance(P, PointGroupElement)
+    assert allclose(P.U, np.eye(2))
+    assert P.conjugate and P.antisymmetry
+    assert allclose(P.R, np.eye(2))
+
+    cons = builder_symmetries['conservation_law']
+    assert isinstance(cons, ContinuousGroupGenerator)
+    assert allclose(cons.U, 2*np.eye(2))
+    assert cons.R is None
+
+    syst = kwant.Builder()
+    builder_symmetries = _get_builder_symmetries(syst)
+    assert len(builder_symmetries) == 0
+
+
+def test_higher_dim():
+    # Test 0D finite system
+    lat = kwant.lattice.cubic(norbs=1)
+    syst = kwant.Builder()
+    syst[lat(0, 0, 0)] = 1
+    syst[lat(1, 1, 0)] = 1
+    syst[lat(0, 1, 1)] = 1
+    syst[lat(1, 0, -1)] = 1
+    syst[lat(0, 0, 0), lat(1, 1, 0)] = -1
+    syst[lat(0, 0, 0), lat(0, 1, 1)] = -1
+    syst[lat(0, 0, 0), lat(1, 0, -1)] = -1
+
+    H = builder_to_model(syst)
+    sg, cs = symmetries(H, prettify=True)
+    assert len(sg) == 2
+    assert len(cs) == 5
+
+    # Test triangular lattice system embedded in 3D
+    sym = TranslationalSymmetry([1, 1, 0], [0, 1, 1])
+    lat = kwant.lattice.cubic(norbs=1)
+    syst = kwant.Builder(symmetry=sym)
+    syst[lat(0, 0, 0)] = 1
+    syst[lat(0, 0, 0), lat(1, 1, 0)] = -1
+    syst[lat(0, 0, 0), lat(0, 1, 1)] = -1
+    syst[lat(0, 0, 0), lat(1, 0, -1)] = -1
+
+    H = builder_to_model(syst)
+    sg, cs = symmetries(H, hexagonal(sympy_R=False), prettify=True)
+    assert len(sg) == 24
+    assert len(cs) == 0
+
+
+def test_graphene_to_kwant():
+
+    norbs = OrderedDict([('A', 1), ('B', 1)])  # A and B atom per unit cell, one orbital each
+    hopping_vectors = [('A', 'B', [1, 0])] # Hopping between neighbouring A and B atoms
+    # Atomic coordinates within the unit cell
+    atom_coords = [(0, 0), (1, 0)]
+    # We set the interatom distance to 1, so the lattice vectors have length sqrt(3)
+    lat_vecs = [(3/2, np.sqrt(3)/2), (3/2, -np.sqrt(3)/2)]
+
+    # Time reversal
+    TR = PointGroupElement(sympy.eye(2), True, False, np.eye(2))
+    # Chiral symmetry
+    C = PointGroupElement(sympy.eye(2), False, True, np.array([[1, 0], [0, -1]]))
+    # Atom A rotates into A, B into B.
+    sphi = 2*sympy.pi/3
+    RC3 = sympy.Matrix([[sympy.cos(sphi), -sympy.sin(sphi)],
+                      [sympy.sin(sphi), sympy.cos(sphi)]])
+    C3 = PointGroupElement(RC3, False, False, np.eye(2))
+
+    # Generate graphene Hamiltonian in Kwant from qsymm
+    symmetries = [C, TR, C3]
+    # Generate using a family
+    family = bloch_family(hopping_vectors, symmetries, norbs)
+    syst_from_family = model_to_builder(family, norbs, lat_vecs, atom_coords, coeffs=None)
+    # Generate using a single Model object
+    g = sympy.Symbol('g', real=True)
+    ham = hamiltonian_from_family(family, coeffs=[g])
+    ham = Model(hamiltonian=ham, momenta=family[0].momenta)
+    syst_from_model = model_to_builder(ham, norbs, lat_vecs, atom_coords)
+
+    # Make the graphene Hamiltonian using kwant only
+    atoms, orbs = zip(*[(atom, norb) for atom, norb in
+                        norbs.items()])
+    # Make the kwant lattice
+    lat = kwant.lattice.general(lat_vecs,
+                                atom_coords,
+                                norbs=orbs)
+    # Store sublattices by name
+    sublattices = {atom: sublat for atom, sublat in
+                   zip(atoms, lat.sublattices)}
+
+    sym = kwant.TranslationalSymmetry(*lat_vecs)
+    bulk = kwant.Builder(sym)
+
+    bulk[ [sublattices['A'](0, 0), sublattices['B'](0, 0)] ] = 0
+
+    def hop(site1, site2, c0):
+        return c0
+
+    bulk[lat.neighbors()] = hop
+
+    fsyst_family = kwant.wraparound.wraparound(syst_from_family).finalized()
+    fsyst_model = kwant.wraparound.wraparound(syst_from_model).finalized()
+    fsyst_kwant = kwant.wraparound.wraparound(bulk).finalized()
+
+    # Check that the energies are identical at random points in the Brillouin zone
+    coeff = 0.5 + np.random.rand()
+    for _ in range(20):
+        kx, ky = 3*np.pi*(np.random.rand(2) - 0.5)
+        params = dict(c0=coeff, k_x=kx, k_y=ky)
+        hamiltonian1 = fsyst_kwant.hamiltonian_submatrix(params=params, sparse=False)
+        hamiltonian2 = fsyst_family.hamiltonian_submatrix(params=params, sparse=False)
+        assert allclose(hamiltonian1, hamiltonian2)
+        params = dict(g=coeff, k_x=kx, k_y=ky)
+        hamiltonian3 = fsyst_model.hamiltonian_submatrix(params=params, sparse=False)
+        assert allclose(hamiltonian2, hamiltonian3)
+
+    # Include random onsites as well
+    one = sympy.numbers.One()
+    onsites = [Model({one: np.array([[1, 0], [0, 0]])}, momenta=family[0].momenta),
+               Model({one: np.array([[0, 0], [0, 1]])}, momenta=family[0].momenta)]
+    family = family + onsites
+    syst_from_family = model_to_builder(family, norbs, lat_vecs, atom_coords, coeffs=None)
+    gs = list(sympy.symbols('g0:%d'%3, real=True))
+    ham = hamiltonian_from_family(family, coeffs=gs)
+    ham = Model(hamiltonian=ham, momenta=family[0].momenta)
+    syst_from_model = model_to_builder(ham, norbs, lat_vecs, atom_coords)
+
+    def onsite_A(site, c1):
+        return c1
+
+    def onsite_B(site, c2):
+        return c2
+
+    bulk[[sublattices['A'](0, 0)]] = onsite_A
+    bulk[[sublattices['B'](0, 0)]] = onsite_B
+
+    fsyst_family = kwant.wraparound.wraparound(syst_from_family).finalized()
+    fsyst_model = kwant.wraparound.wraparound(syst_from_model).finalized()
+    fsyst_kwant = kwant.wraparound.wraparound(bulk).finalized()
+
+    # Check equivalence of the Hamiltonian at random points in the BZ
+    coeffs = 0.5 + np.random.rand(3)
+    for _ in range(20):
+        kx, ky = 3*np.pi*(np.random.rand(2) - 0.5)
+        params = dict(c0=coeffs[0], c1=coeffs[1], c2=coeffs[2], k_x=kx, k_y=ky)
+        hamiltonian1 = fsyst_kwant.hamiltonian_submatrix(params=params, sparse=False)
+        hamiltonian2 = fsyst_family.hamiltonian_submatrix(params=params, sparse=False)
+        assert allclose(hamiltonian1, hamiltonian2)
+        params = dict(g0=coeffs[0], g1=coeffs[1], g2=coeffs[2], k_x=kx, k_y=ky)
+        hamiltonian3 = fsyst_model.hamiltonian_submatrix(params=params, sparse=False)
+        assert allclose(hamiltonian2, hamiltonian3)
+
+
+def test_wraparound_convention():
+    # Test that it matches exactly kwant.wraparound convention
+    # Make the graphene Hamiltonian using kwant only
+    norbs = OrderedDict([('A', 1), ('B', 1)])  # A and B atom per unit cell, one orbital each
+    atoms, orbs = zip(*[(atom, norb) for atom, norb in
+                        norbs.items()])
+    # Atomic coordinates within the unit cell
+    atom_coords = [(0, 0), (1, 0)]
+    # We set the interatom distance to 1, so the lattice vectors have length sqrt(3)
+    lat_vecs = [(3/2, np.sqrt(3)/2), (3/2, -np.sqrt(3)/2)]
+    # Make the kwant lattice
+    lat = kwant.lattice.general(lat_vecs,
+                                atom_coords,
+                                norbs=orbs)
+    # Store sublattices by name
+    sublattices = {atom: sublat for atom, sublat in
+                   zip(atoms, lat.sublattices)}
+
+    sym = kwant.TranslationalSymmetry(*lat_vecs)
+    bulk = kwant.Builder(sym)
+
+    bulk[ [sublattices['A'](0, 0), sublattices['B'](0, 0)] ] = 0
+
+    def hop(site1, site2, c0):
+        return c0
+
+    bulk[lat.neighbors()] = hop
+
+    wrapped = kwant.wraparound.wraparound(bulk).finalized()
+    ham2 = builder_to_model(bulk, real_space=False)
+    # Check that the Hamiltonians are identical at random points in the Brillouin zone
+    H1 = wrapped.hamiltonian_submatrix
+    H2 = ham2.lambdify()
+    coeffs = 0.5 + np.random.rand(1)
+    for _ in range(20):
+        kx, ky = 3*np.pi*(np.random.rand(2) - 0.5)
+        params = dict(c0=coeffs[0], k_x=kx, k_y=ky)
+        h1, h2 = H1(params=params), H2(**params)
+        assert allclose(h1, h2), (h1, h2)
+
+
+
+def test_inverse_transform():
+    # Define family on square lattice
+    s = spin_matrices(1/2)
+    # Time reversal
+    TR = PointGroupElement(np.eye(2), True, False,
+                           spin_rotation(2 * np.pi * np.array([0, 1/2, 0]), s))
+    # Mirror symmetry
+    Mx = PointGroupElement(np.array([[-1, 0], [0, 1]]), False, False,
+                           spin_rotation(2 * np.pi * np.array([1/2, 0, 0]), s))
+    # Fourfold
+    C4 = PointGroupElement(np.array([[0, 1], [-1, 0]]), False, False,
+                           spin_rotation(2 * np.pi * np.array([0, 0, 1/4]), s))
+    symmetries = [TR, Mx, C4]
+
+    # One site per unit cell
+    norbs = OrderedDict([('A', 2)])
+    # Hopping to a neighbouring atom one primitive lattice vector away
+    hopping_vectors = [('A', 'A', [1, 0])]
+    # Make family
+    family = bloch_family(hopping_vectors, symmetries, norbs)
+    fam = hamiltonian_from_family(family, tosympy=False)
+    # Atomic coordinates within the unit cell
+    atom_coords = [(0, 0)]
+    lat_vecs = [(1, 0), (0, 1)]
+    syst = model_to_builder(fam, norbs, lat_vecs, atom_coords)
+    # Convert it back
+    ham2 = builder_to_model(syst).tomodel(nsimplify=True)
+    # Check that it's the same as the original
+    assert fam == ham2
+
+    # Check that the Hamiltonians are identical at random points in the Brillouin zone
+    sysw = kwant.wraparound.wraparound(syst).finalized()
+    H1 = sysw.hamiltonian_submatrix
+    H2 = ham2.lambdify()
+    H3 = fam.lambdify()
+    coeffs = 0.5 + np.random.rand(3)
+    for _ in range(20):
+        kx, ky = 3*np.pi*(np.random.rand(2) - 0.5)
+        params = dict(c0=coeffs[0], c1=coeffs[1], c2=coeffs[2], k_x=kx, k_y=ky)
+        assert allclose(H1(params=params), H2(**params))
+        assert allclose(H1(params=params), H3(**params))
+
+
+def test_consistency_kwant():
+    """Make a random 1D Model, convert it to a builder, and compare
+    the Bloch representation of the Model with that which Kwant uses
+    in wraparound and in Bands. Then, convert the builder back to a Model
+    and compare with the original Model.
+    For comparison, we also make the system using Kwant only.
+    """
+    orbs = 4
+    T = np.random.rand(2*orbs, 2*orbs) + 1j*np.random.rand(2*orbs, 2*orbs)
+    H = np.random.rand(2*orbs, 2*orbs) + 1j*np.random.rand(2*orbs, 2*orbs)
+    H += H.T.conj()
+
+    # Make the 1D Model manually using only qsymm features.
+    c0, c1 = sympy.symbols('c0 c1', real=True)
+    kx = _commutative_momenta[0]
+
+    Ham = Model({c0 * e**(-I*kx): T}, momenta=[0])
+    Ham += Ham.T().conj()
+    Ham += Model({c1: H}, momenta=[0])
+
+    # Two superimposed atoms, same number of orbitals on each
+    norbs = OrderedDict([('A', orbs), ('B', orbs)])
+    atom_coords = [(0.3, ), (0.3, )]
+    lat_vecs = [(1, )] # Lattice vector
+
+    # Make a Kwant builder out of the qsymm Model
+    model_syst = model_to_builder(Ham, norbs, lat_vecs, atom_coords)
+    fmodel_syst = model_syst.finalized()
+
+    # Make the same system manually using only Kwant features.
+    lat = kwant.lattice.general(np.array([[1.]]),
+                            [(0., )],
+                            norbs=2*orbs)
+    kwant_syst = kwant.Builder(kwant.TranslationalSymmetry(*lat.prim_vecs))
+
+    def onsite(site, c1):
+        return c1*H
+
+    def hopping(site1, site2, c0):
+        return c0*T
+
+    sublat = lat.sublattices[0]
+    kwant_syst[sublat(0,)] = onsite
+    hopp = kwant.builder.HoppingKind((1, ), sublat)
+    kwant_syst[hopp] = hopping
+    fkwant_syst = kwant_syst.finalized()
+
+    # Make sure we are consistent with bands calculations in kwant
+    # The Bloch Hamiltonian used in Kwant for the bands computation
+    # is h(k) = exp(-i*k)*hop + onsite + exp(i*k)*hop.T.conj.
+    # We also check that all is consistent with wraparound
+    coeffs = (0.7, 1.2)
+    params = dict(c0 = coeffs[0], c1 = coeffs[1])
+    kwant_hop = fkwant_syst.inter_cell_hopping(params=params)
+    kwant_onsite = fkwant_syst.cell_hamiltonian(params=params)
+    model_kwant_hop = fmodel_syst.inter_cell_hopping(params=params)
+    model_kwant_onsite = fmodel_syst.cell_hamiltonian(params=params)
+
+    assert allclose(model_kwant_hop, coeffs[0]*T)
+    assert allclose(model_kwant_hop, kwant_hop)
+    assert allclose(model_kwant_onsite, kwant_onsite)
+
+    h_model_kwant = (lambda k: np.exp(-1j*k)*model_kwant_hop + model_kwant_onsite +
+                     np.exp(1j*k)*model_kwant_hop.T.conj()) # As in kwant.Bands
+    h_model = Ham.lambdify()
+    wsyst = kwant.wraparound.wraparound(model_syst).finalized()
+    for _ in range(20):
+        k = (np.random.rand() - 0.5)*2*np.pi
+        assert allclose(h_model_kwant(k), h_model(coeffs[0], coeffs[1], k))
+        params['k_x'] = k
+        h_wrap = wsyst.hamiltonian_submatrix(params=params)
+        assert allclose(h_model(coeffs[0], coeffs[1], k), h_wrap)
+
+    # Get the model back from the builder
+    # From the Kwant builder based on original Model
+    Ham1 = builder_to_model(model_syst, momenta=Ham.momenta).tomodel(nsimplify=True)
+    # From the pure Kwant builder
+    Ham2 = builder_to_model(kwant_syst, momenta=Ham.momenta).tomodel(nsimplify=True)
+    assert Ham == Ham1
+    assert Ham == Ham2
+
+
+def test_find_builder_discrete_symmetries():
+    symm_class = ['AI', 'D', 'AIII', 'BDI']
+    class_dict = {'AI': ['time_reversal'],
+                'D': ['particle_hole'],
+                'AIII': ['chiral'],
+                'BDI': ['time_reversal', 'particle_hole', 'chiral']}
+    sym_dict = {'time_reversal': qsymm.PointGroupElement(np.eye(2), True, False, None),
+                'particle_hole': qsymm.PointGroupElement(np.eye(2), True, True, None),
+                'chiral': qsymm.PointGroupElement(np.eye(2), False, True, None)}
+    n = 4
+    rng = 11
+    for sym in symm_class:
+        # Random Hamiltonian in the symmetry class
+        h_ons = kwant.rmt.gaussian(n, sym, rng=rng)
+        h_hop = 10 * kwant.rmt.gaussian(2*n, sym, rng=rng)[:n, n:]
+        # Make a Kwant builder in the symmetry class and find its symmetries
+        lat = kwant.lattice.square(norbs=n)
+        bulk = kwant.Builder(TranslationalSymmetry([1, 0], [0, 1]))
+        bulk[lat(0, 0)] = h_ons
+        bulk[kwant.builder.HoppingKind((1, 0), lat)] = h_hop
+        bulk[kwant.builder.HoppingKind((0, 1), lat)] = h_hop
+        builder_symmetries = find_builder_symmetries(bulk, spatial_symmetries=True, prettify=True)
+
+        # Equality of symmetries ignores unitary part
+        fourfold_rotation = qsymm.PointGroupElement(np.array([[0, 1],[1, 0]]), False, False, None)
+        assert fourfold_rotation in builder_symmetries
+        class_symmetries = class_dict[sym]
+        for class_symmetry in class_symmetries:
+            assert sym_dict[class_symmetry] in builder_symmetries
+
+
+def random_onsite_hop(n, rng=0):
+    rng = ensure_rng(rng)
+    onsite = rng.randn(n, n) + 1j * rng.randn(n, n)
+    onsite = onsite + onsite.T.conj()
+    hop = rng.rand(n, n) + 1j * rng.rand(n, n)
+    return onsite, hop
+
+
+def test_find_cons_law():
+    sy = np.array([[0, -1j], [1j, 0]])
+
+    n = 3
+    lat = kwant.lattice.chain(norbs=2*n)
+    syst = kwant.Builder()
+    rng = 1337
+    ons, hop = random_onsite_hop(n, rng=rng)
+
+    syst[lat(0)] = np.kron(sy, ons)
+    syst[lat(1)] = np.kron(sy, ons)
+    syst[lat(1), lat(0)] = np.kron(sy, hop)
+
+    builder_symmetries = find_builder_symmetries(syst, spatial_symmetries=False, prettify=True)
+    onsites = [symm for symm in builder_symmetries if
+               isinstance(symm, qsymm.ContinuousGroupGenerator) and symm.R is None]
+    mham = builder_to_model(syst)
+    assert not any([len(symm.apply(mham)) for symm in onsites])
+
+def test_basis_ordering():
+    symm_class = ['AI', 'D', 'AIII', 'BDI']
+    n = 2
+    rng = 12
+    for sym in symm_class:
+        # Make a small finite system in the symmetry class, finalize it and check
+        # that the basis is consistent.
+        h_ons = kwant.rmt.gaussian(n, sym, rng=rng)
+        h_hop = 10 * kwant.rmt.gaussian(2*n, sym, rng=rng)[:n, n:]
+        lat = kwant.lattice.square(norbs=n)
+        bulk = kwant.Builder(TranslationalSymmetry([1, 0], [0, 1]))
+        bulk[lat(0, 0)] = h_ons
+        bulk[kwant.builder.HoppingKind((1, 0), lat)] = h_hop
+        bulk[kwant.builder.HoppingKind((0, 1), lat)] = h_hop
+
+        def rect(site):
+            x, y = site.pos
+            return (0 <= x < 2) and (0 <= y < 3)
+
+        square = kwant.Builder()
+        square.fill(bulk, lambda site: rect(site), (0, 0),
+                    max_sites=float('inf'))
+
+        # Find the symmetries of the square
+        builder_symmetries = find_builder_symmetries(square,
+                                                     spatial_symmetries=False,
+                                                     prettify=True)
+        # Finalize the square, extract Hamiltonian
+        fsquare = square.finalized()
+        ham = fsquare.hamiltonian_submatrix()
+
+        # Check manually that the found symmetries are in the same basis as the
+        # finalized system
+        for symmetry in builder_symmetries:
+            U = symmetry.U
+            if isinstance(symmetry, qsymm.ContinuousGroupGenerator):
+                assert symmetry.R is None
+                assert allclose(U.dot(ham), ham.dot(U))
+            else:
+                if symmetry.conjugate:
+                    left = U.dot(ham.conj())
+                else:
+                    left = U.dot(ham)
+                if symmetry.antisymmetry:
+                    assert allclose(left, -ham.dot(U))
+                else:
+                    assert allclose(left, ham.dot(U))