diff --git a/doc/source/reference/kwant.builder.rst b/doc/source/reference/kwant.builder.rst
index 6211991e4c1b652a09aab18de84a36b55f53469a..2554cf9acffb2342b3bfdf26e75f5f6a18ba623c 100644
--- a/doc/source/reference/kwant.builder.rst
+++ b/doc/source/reference/kwant.builder.rst
@@ -23,7 +23,6 @@ Abstract base classes
 .. autosummary::
    :toctree: generated/
 
-   Symmetry
    Lead
 
 Functions
diff --git a/doc/source/reference/kwant.system.rst b/doc/source/reference/kwant.system.rst
index 2a9bd89c2529d65e5cc3cad63a144dc31e7a03f1..ef6143a865281dcf64e8787cb833845c82988d0d 100644
--- a/doc/source/reference/kwant.system.rst
+++ b/doc/source/reference/kwant.system.rst
@@ -28,6 +28,14 @@ Sites
    SiteArray
    SiteFamily
 
+Symmetries
+----------
+.. autosummary::
+   :toctree: generated/
+
+    Symmetry
+    NoSymmetry
+
 Systems
 -------
 .. autosummary::
diff --git a/kwant/_system.pyx b/kwant/_system.pyx
index f85bed6b3863683b4174acd8a1b07b09dae6f7a4..ac7875b29dd77397102e5dde9b7967cf700312de 100644
--- a/kwant/_system.pyx
+++ b/kwant/_system.pyx
@@ -548,13 +548,12 @@ def vectorized_cell_hamiltonian(self, args=(), sparse=False, *, params=None):
     # Site array where next cell starts
     next_cell = bisect.bisect(site_offsets, self.cell_size) - 1
 
-    def inside_cell(term):
-        (to_which, from_which), _= self.subgraphs[term.subgraph]
-        return to_which < next_cell and from_which < next_cell
+    def inside_fd(term):
+        return all(s == 0 for s in term.symmetry_element)
 
     cell_terms = [
         n for n, t in enumerate(self.terms)
-        if inside_cell(t)
+        if inside_fd(t)
     ]
 
     subgraphs = [
@@ -593,34 +592,38 @@ def vectorized_cell_hamiltonian(self, args=(), sparse=False, *, params=None):
 def vectorized_inter_cell_hopping(self, args=(), sparse=False, *, params=None):
     """Hopping Hamiltonian between two cells of the infinite system.
 
+    This method returns a complex matrix that represents the hopping from
+    the *interface sites* of unit cell ``n - 1`` to *all* the sites of
+    unit cell ``n``. It is therefore generally a *rectangular* matrix of
+    shape ``(N_uc, N_iface)`` where ``N_uc`` is the number of orbitals
+    in the unit cell, and ``N_iface`` is the number of orbitals on the
+    *interface* sites (i.e. the sites with hoppings *to* the next unit cell).
+
     Providing positional arguments via 'args' is deprecated,
     instead, provide named parameters as a dictionary via 'params'.
     """
 
-    # Take advantage of the fact that there are distinct entries in
-    # onsite_terms for sites inside the unit cell, and distinct entries
-    # in onsite_terms for hoppings between the unit cell sites and
-    # interface sites. This way we only need to check the first entry
-    # in onsite/hopping_terms
-
     site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
-    # Site array where next cell starts
-    next_cell = bisect.bisect(site_offsets, self.cell_size) - 1
 
+    # This method is only meaningful for systems with a 1D translational
+    # symmetry, and we use this fact in several places
+    assert all(len(t.symmetry_element) == 1 for t in self.terms)
+
+    # Symmetry element -1 means hoppings *from* the *previous*
+    # unit cell. These are directly the hoppings we wish to return.
     inter_cell_hopping_terms = [
         n for n, t in enumerate(self.terms)
-        # *from* site is in interface
-        if (self.subgraphs[t.subgraph][0][1] >= next_cell
-            and self.subgraphs[t.subgraph][0][0] < next_cell)
+        if t.symmetry_element[0] == -1
     ]
+    # Symmetry element +1 means hoppings *from* the *next* unit cell.
+    # These are related by translational symmetry to hoppings *to*
+    # the *previous* unit cell. We therefore need the *reverse*
+    # (and conjugate) of these hoppings.
     reversed_inter_cell_hopping_terms = [
         n for n, t in enumerate(self.terms)
-        # *to* site is in interface
-        if (self.subgraphs[t.subgraph][0][0] >= next_cell
-            and self.subgraphs[t.subgraph][0][1] < next_cell)
+        if t.symmetry_element[0] == +1
     ]
 
-    # Evaluate inter-cell hoppings only
     inter_cell_hams = [
         self.hamiltonian_term(n, args=args, params=params)
         for n in inter_cell_hopping_terms
@@ -642,18 +645,21 @@ def vectorized_inter_cell_hopping(self, args=(), sparse=False, *, params=None):
         for n in reversed_inter_cell_hopping_terms
     ]
 
-    # All the 'from' sites are in the previous domain, but to build a
-    # matrix we need to get the equivalent sites in the fundamental domain.
-    # We set the 'from' site array to the one from the fundamental domain.
-    subgraphs = [
-        ((to_sa, from_sa - next_cell), (to_off, from_off))
-        for (to_sa, from_sa), (to_off, from_off) in subgraphs
-    ]
-
     _, norbs, orb_offsets = self.site_ranges.transpose()
 
-    shape = (orb_offsets[next_cell],
-             orb_offsets[len(self.site_arrays) - next_cell])
+    # SiteArrays containing interface sites appear before SiteArrays
+    # containing non-interface sites, so the max of the site array
+    # indices that appear in the 'from' site arrays of the inter-cell
+    # hoppings allows us to get the number of interface orbitals.
+    last_iface_site_array = max(
+        from_site_array for (_, from_site_array), _ in subgraphs
+    )
+    iface_norbs = orb_offsets[last_iface_site_array + 1]
+    fd_norbs = orb_offsets[-1]
+
+    # TODO: return a square matrix when we no longer need to maintain
+    #       backwards compatibility with unvectorized systems.
+    shape = (fd_norbs, iface_norbs)
     func = _vectorized_make_sparse if sparse else _vectorized_make_dense
     mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape)
     return mat
diff --git a/kwant/builder.py b/kwant/builder.py
index 26e1d60525084f10f42733fe92c3306cdd14b3a2..252c22414fc4e56b35ae39da6869ae5e7a313559 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -21,7 +21,7 @@ import tinyarray as ta
 import numpy as np
 from scipy import sparse
 from . import system, graph, KwantDeprecationWarning, UserCodeError
-from .system import Site, SiteArray, SiteFamily
+from .system import Site, SiteArray, SiteFamily, Symmetry, NoSymmetry
 from .linalg import lll
 from .operator import Density
 from .physics import DiscreteSymmetry, magnetic_gauge
@@ -29,7 +29,7 @@ from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
                       interleave, deprecate_args, memoize)
 
 
-__all__ = ['Builder', 'Symmetry', 'HoppingKind', 'Lead',
+__all__ = ['Builder', 'HoppingKind', 'Lead',
            'BuilderLead', 'SelfEnergyLead', 'ModesLead', 'add_peierls_phase']
 
 
@@ -64,182 +64,6 @@ def validate_hopping(hopping):
         raise ValueError("A hopping connects the following site to itself:\n"
                          "{0}".format(a))
 
-
-
-################ Symmetries
-
-class Symmetry(metaclass=abc.ABCMeta):
-    """Abstract base class for spatial symmetries.
-
-    Many physical systems possess a discrete spatial symmetry, which results in
-    special properties of these systems.  This class is the standard way to
-    describe discrete spatial symmetries in Kwant.  An instance of this class
-    can be passed to a `Builder` instance at its creation.  The most important
-    kind of symmetry is translational symmetry, used to define scattering
-    leads.
-
-    Each symmetry has a fundamental domain -- a set of sites and hoppings,
-    generating all the possible sites and hoppings upon action of symmetry
-    group elements.  A class derived from `Symmetry` has to implement mapping
-    of any site or hopping into the fundamental domain, applying a symmetry
-    group element to a site or a hopping, and a method `which` to determine the
-    group element bringing some site from the fundamental domain to the
-    requested one.  Additionally, it has to have a property `num_directions`
-    returning the number of independent symmetry group generators (number of
-    elementary periods for translational symmetry).
-
-    A ``ValueError`` must be raised by the symmetry class whenever a symmetry
-    is used together with sites whose site family is not compatible with it.  A
-    typical example of this is when the vector defining a translational
-    symmetry is not a lattice vector.
-
-    The type of the domain objects as handled by the methods of this class is
-    not specified.  The only requirement is that it must support the unary
-    minus operation.  The reference implementation of `to_fd()` is hence
-    `self.act(-self.which(a), a, b)`.
-    """
-
-    @abc.abstractproperty
-    def num_directions(self):
-        """Number of elementary periods of the symmetry."""
-        pass
-
-    @abc.abstractmethod
-    def which(self, site):
-        """Calculate the domain of the site.
-
-        Parameters
-        ----------
-        site : `~kwant.system.Site` or `~kwant.system.SiteArray`
-
-        Returns
-        -------
-        group_element : tuple or sequence of tuples
-            A single tuple if ``site`` is a Site, or a sequence of tuples if
-            ``site`` is a SiteArray.  The group element(s) whose action
-            on a certain site(s) from the fundamental domain will result
-            in the given ``site``.
-        """
-        pass
-
-    @abc.abstractmethod
-    def act(self, element, a, b=None):
-        """Act with symmetry group element(s) on site(s) or hopping(s).
-
-        Parameters
-        ----------
-        element : tuple or sequence of tuples
-            Group element(s) with which to act on the provided site(s)
-            or hopping(s)
-        a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
-            If Site then ``element`` is a single tuple, if SiteArray then
-            ``element`` is a sequence of tuples. If only ``a`` is provided then
-            ``element`` acts on the site(s) of ``a``. If ``b`` is also provided
-            then ``element`` acts on the hopping(s) ``(a, b)``.
-        """
-        pass
-
-    def to_fd(self, a, b=None):
-        """Map a site or hopping to the fundamental domain.
-
-        Parameters
-        ----------
-        a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
-
-        If ``b`` is None, return a site equivalent to ``a`` within the
-        fundamental domain.  Otherwise, return a hopping equivalent to ``(a,
-        b)`` but where the first element belongs to the fundamental domain.
-
-        Equivalent to `self.act(-self.which(a), a, b)`.
-        """
-        return self.act(-self.which(a), a, b)
-
-    def in_fd(self, site):
-        """Tell whether ``site`` lies within the fundamental domain.
-
-        Parameters
-        ----------
-        site : `~kwant.system.Site` or `~kwant.system.SiteArray`
-
-        Returns
-        -------
-        in_fd : bool or sequence of bool
-            single bool if ``site`` is a Site, or a sequence of
-            bool if ``site`` is a SiteArray. In the latter case
-            we return whether each site in the SiteArray is in
-            the fundamental domain.
-        """
-        if isinstance(site, Site):
-            for d in self.which(site):
-                if d != 0:
-                    return False
-            return True
-        elif isinstance(site, SiteArray):
-            which = self.which(site)
-            return np.logical_and.reduce(which != 0, axis=1)
-        else:
-            raise TypeError("'site' must be a Site or SiteArray")
-
-    @abc.abstractmethod
-    def subgroup(self, *generators):
-        """Return the subgroup generated by a sequence of group elements."""
-        pass
-
-    @abc.abstractmethod
-    def has_subgroup(self, other):
-        """Test whether `self` has the subgroup `other`...
-
-        or, in other words, whether `other` is a subgroup of `self`.  The
-        reason why this is the abstract method (and not `is_subgroup`) is that
-        in general it's not possible for a subgroup to know its supergroups.
-
-        """
-        pass
-
-
-class NoSymmetry(Symmetry):
-    """A symmetry with a trivial symmetry group."""
-
-    def __eq__(self, other):
-        return isinstance(other, NoSymmetry)
-
-    def __ne__(self, other):
-        return not self.__eq__(other)
-
-    def __repr__(self):
-        return 'NoSymmetry()'
-
-    @property
-    def num_directions(self):
-        return 0
-
-    periods = ()
-
-    _empty_array = ta.array((), int)
-
-    def which(self, site):
-        return self._empty_array
-
-    def act(self, element, a, b=None):
-        if element:
-            raise ValueError('`element` must be empty for NoSymmetry.')
-        return a if b is None else (a, b)
-
-    def to_fd(self, a, b=None):
-        return a if b is None else (a, b)
-
-    def in_fd(self, site):
-        return True
-
-    def subgroup(self, *generators):
-        if any(generators):
-            raise ValueError('Generators must be empty for NoSymmetry.')
-        return NoSymmetry(generators)
-
-    def has_subgroup(self, other):
-        return isinstance(other, NoSymmetry)
-
-
 
 ################ Hopping kinds
 
@@ -658,7 +482,7 @@ class Builder:
 
     Parameters
     ----------
-    symmetry : `Symmetry` or `None`
+    symmetry : `~kwant.system.Symmetry` or `None`
         The spatial symmetry of the system.
     conservation_law : 2D array, dictionary, function, or `None`
         An onsite operator with integer eigenvalues that commutes with the
@@ -707,8 +531,8 @@ class Builder:
     that if ``builder[a, b]`` has been set, there is no need to set
     ``builder[b, a]``.
 
-    Builder instances can be made to automatically respect a `Symmetry` that is
-    passed to them during creation.  The behavior of builders with a symmetry
+    Builder instances can be made to automatically respect a `~kwant.system.Symmetry`
+    that is passed to them during creation.  The behavior of builders with a symmetry
     is slightly more sophisticated: all keys are mapped to the fundamental
     domain of the symmetry before storing them.  This may produce confusing
     results when neighbors of a site are queried.
@@ -1657,7 +1481,7 @@ class Builder:
         system to be returned.
 
         Currently, only Builder instances without or with a 1D translational
-        `Symmetry` can be finalized.
+        `~kwant.system.Symmetry` can be finalized.
         """
         if self.symmetry.num_directions == 0:
             if self.vectorize:
@@ -2025,7 +1849,7 @@ class _VectorizedFinalizedBuilderMixin(_FinalizedBuilderMixin):
             "elements at once using 'hamiltonian_term'.",
             KwantDeprecationWarning
         )
-        site_offsets = np.cumsum([0] + [len(s) for s in self.site_arrays])
+        site_offsets, _, _ = self.site_ranges.transpose()
         if i == j:
             which_term = self._onsite_term_by_site_id[i]
             (w, _), (off, _) = self.subgraphs[self.terms[which_term].subgraph]
@@ -2036,6 +1860,15 @@ class _VectorizedFinalizedBuilderMixin(_FinalizedBuilderMixin):
             return onsite[0]
         else:
             edge_id = self.graph.first_edge_id(i, j)
+            # Map sites in previous cell to fundamental domain; vectorized
+            # infinite systems only store sites in the FD. We already know
+            # that this hopping is between unit cells because this is encoded
+            # in the edge_id and term id.
+            # Using 'is_infinite' is a bit of a hack, but 'hamiltonian' is
+            # deprecated anyway, and refactoring everything to avoid this check
+            # is not worth it.
+            if system.is_infinite(self):
+                i, j = i % self.cell_size, j % self.cell_size
             which_term = self._hopping_term_by_edge_id[edge_id]
             herm_conj = which_term < 0
             if herm_conj:
@@ -2087,10 +1920,15 @@ class _VectorizedFinalizedBuilderMixin(_FinalizedBuilderMixin):
         to_family = self.site_arrays[to_which].family
         to_tags = self.site_arrays[to_which].tags
         to_site_array = SiteArray(to_family, to_tags[to_off])
+        site_arrays = (to_site_array,)
         if not is_onsite:
             from_family = self.site_arrays[from_which].family
             from_tags = self.site_arrays[from_which].tags
-            from_site_array = SiteArray(from_family, from_tags[from_off])
+            from_site_array = self.symmetry.act(
+                term.symmetry_element,
+                SiteArray(from_family, from_tags[from_off])
+            )
+            site_arrays = (to_site_array, from_site_array)
 
         # Construct args from params
         if params:
@@ -2107,20 +1945,10 @@ class _VectorizedFinalizedBuilderMixin(_FinalizedBuilderMixin):
                        ', '.join(map('"{}"'.format, missing)))
                 raise TypeError(''.join(msg))
 
-        if is_onsite:
-            try:
-                ham = val(to_site_array, *args)
-            except Exception as exc:
-                _raise_user_error(exc, val)
-        else:
-            try:
-                ham = val(
-                        *self.symmetry.to_fd(
-                            to_site_array,
-                            from_site_array),
-                        *args)
-            except Exception as exc:
-                _raise_user_error(exc, val)
+        try:
+            ham = val(*site_arrays, *args)
+        except Exception as exc:
+            _raise_user_error(exc, val)
 
         expected_shape = (
             len(to_site_array),
@@ -2250,8 +2078,8 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
 
         graph = _make_graph(builder.H, id_by_site)
 
-        finalized_leads, lead_interfaces, lead_paddings =\
-            _finalize_leads(builder.leads, id_by_site)
+        finalized_leads, lead_interfaces, lead_paddings = _finalize_leads(
+            builder.leads, id_by_site)
 
         # Because many onsites/hoppings share the same (value, parameter)
         # pairs, we keep them in a cache so that we only store a given pair
@@ -2485,9 +2313,15 @@ def _make_onsite_terms(builder, sites, site_arrays, term_offset):
         err if isinstance(err, Exception) else None
         for err in onsite_term_parameters
     ]
+
+    # We store sites in the fundamental domain, so the symmetry element
+    # is the same for all onsite terms; it's just the identity.
+    identity = ta.array([0] * builder.symmetry.num_directions)
+
     onsite_terms = [
         system.Term(
             subgraph=term_offset + i,
+            symmetry_element=identity,
             hermitian=False,
             parameters=(
                 params if not isinstance(params, Exception) else None
@@ -2511,6 +2345,11 @@ def _make_hopping_terms(builder, graph, sites, site_arrays, cell_size, term_offs
     #   Maps hopping edge IDs to the number of the term that the hopping
     #   is a part of. For Hermitian conjugate hoppings "-term_number -1"
     #   is stored instead.
+    #
+    # NOTE: 'graph' contains site indices >= cell_size, however 'sites'
+    #       and 'site_arrays' only contain information about the sites
+    #       in the fundamental domain.
+    sym = builder.symmetry
 
     site_offsets = np.cumsum([0] + [len(sa) for sa in site_arrays])
 
@@ -2523,11 +2362,27 @@ def _make_hopping_terms(builder, graph, sites, site_arrays, cell_size, term_offs
     hopping_to_term_nr = collections.OrderedDict()
     _hopping_term_by_edge_id = []
     for tail, head in graph:
-        tail_site, head_site = sites[tail], sites[head]
-        if tail >= cell_size:
-            # The tail belongs to the previous domain.  Find the
-            # corresponding hopping with the tail in the fund. domain.
-            tail_site, head_site = builder.symmetry.to_fd(tail_site, head_site)
+        # map 'tail' and 'head' so that they are both in the FD, and
+        # assign the symmetry element to apply to 'sites[head]'
+        # to make the actual hopping. Here we assume that the symmetry
+        # may only be NoSymmetry or 1D TranslationalSymmetry.
+        if isinstance(sym, NoSymmetry):
+            tail_site = sites[tail]
+            head_site = sites[head]
+            sym_el = ta.array([])
+        else:
+            if tail < cell_size and head < cell_size:
+                sym_el = ta.array([0])
+            elif head >= cell_size:
+                assert tail < cell_size
+                head = head - cell_size
+                sym_el = ta.array([-1])
+            else:
+                tail = tail - cell_size
+                sym_el = ta.array([1])
+            tail_site = sites[tail]
+            head_site = sym.act(sym_el, sites[head])
+
         val = builder._get_edge(tail_site, head_site)
         herm_conj = val is Other
         if herm_conj:
@@ -2538,11 +2393,11 @@ def _make_hopping_terms(builder, graph, sites, site_arrays, cell_size, term_offs
         head_site_array = bisect.bisect(site_offsets, head) - 1
         # "key" uniquely identifies a hopping term.
         if const_val:
-            key = (None, tail_site_array, head_site_array)
+            key = (None, sym_el, tail_site_array, head_site_array)
         else:
-            key = (id(val), tail_site_array, head_site_array)
+            key = (id(val), sym_el, tail_site_array, head_site_array)
         if herm_conj:
-            key = (key[0], head_site_array, tail_site_array)
+            key = (key[0], -sym_el, head_site_array, tail_site_array)
 
         if key not in hopping_to_term_nr:
             # start a new term
@@ -2584,7 +2439,7 @@ def _make_hopping_terms(builder, graph, sites, site_arrays, cell_size, term_offs
                         site_arrays[from_sa].family.norbs,
                     ),
                 )
-            for (_, to_sa, from_sa), val in
+            for (_, _, to_sa, from_sa), val in
             zip(hopping_to_term_nr.keys(), hopping_term_values)
         ]
         # Sort the hoppings in each term, and also sort the values
@@ -2595,8 +2450,8 @@ def _make_hopping_terms(builder, graph, sites, site_arrays, cell_size, term_offs
             zip(hopping_subgraphs, hopping_term_values)))
     # Store site array index and site offsets rather than sites themselves
     tmp = []
-    for (_, tail_which, head_which), h in zip(hopping_to_term_nr,
-                                              hopping_subgraphs):
+    for (_, _, tail_which, head_which), h in zip(hopping_to_term_nr,
+                                                 hopping_subgraphs):
         start = (site_offsets[tail_which], site_offsets[head_which])
         # Transpose to get a pair of arrays rather than array of pairs
         # We use the fact that the underlying array is stored in
@@ -2613,15 +2468,20 @@ def _make_hopping_terms(builder, graph, sites, site_arrays, cell_size, term_offs
         err if isinstance(err, Exception) else None
         for err in hopping_term_parameters
     ]
+    term_symmetry_elements = [
+        sym_el for (_, sym_el, *_) in hopping_to_term_nr.keys()
+    ]
     hopping_terms = [
         system.Term(
             subgraph=term_offset + i,
+            symmetry_element=sym_el,
             hermitian=True,  # Builders are always Hermitian
             parameters=(
                 params if not isinstance(params, Exception) else None
             ),
         )
-        for i, params in enumerate(hopping_term_parameters)
+        for i, (sym_el, params) in
+        enumerate(zip(term_symmetry_elements, hopping_term_parameters))
     ]
     _hopping_term_by_edge_id = np.array(_hopping_term_by_edge_id)
 
@@ -2662,21 +2522,21 @@ class FiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.FiniteVect
 
         graph = _make_graph(builder.H, id_by_site)
 
-        finalized_leads, lead_interfaces, lead_paddings =\
-            _finalize_leads(builder.leads, id_by_site)
+        finalized_leads, lead_interfaces, lead_paddings = _finalize_leads(
+            builder.leads, id_by_site)
 
         del id_by_site  # cleanup due to large size
 
         site_arrays = _make_site_arrays(builder.H)
 
         (onsite_subgraphs, onsite_terms, onsite_term_values,
-         onsite_term_errors, _onsite_term_by_site_id) =\
-            _make_onsite_terms(builder, sites, site_arrays, term_offset=0)
+         onsite_term_errors, _onsite_term_by_site_id) = _make_onsite_terms(
+             builder, sites, site_arrays, term_offset=0)
 
         (hopping_subgraphs, hopping_terms, hopping_term_values,
-         hopping_term_errors, _hopping_term_by_edge_id) =\
-            _make_hopping_terms(builder, graph, sites, site_arrays,
-                                len(sites), term_offset=len(onsite_terms))
+         hopping_term_errors, _hopping_term_by_edge_id) = _make_hopping_terms(
+             builder, graph, sites, site_arrays, len(sites),
+             term_offset=len(onsite_terms))
 
         # Construct the combined onsite/hopping term datastructures
         subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
@@ -2834,8 +2694,8 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
         sym = builder.symmetry
         assert sym.num_directions == 1
 
-        lsites_with, lsites_without, interface =\
-            _make_lead_sites(builder, interface_order)
+        lsites_with, lsites_without, interface = _make_lead_sites(
+            builder, interface_order)
         cell_size = len(lsites_with) + len(lsites_without)
 
         # we previously sorted the interface, so don't sort it again
@@ -2943,12 +2803,13 @@ class InfiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.Infinite
         assert sym.num_directions == 1
         assert builder.vectorize
 
-        lsites_with, lsites_without, interface =\
-            _make_lead_sites(builder, interface_order)
+        lsites_with, lsites_without, interface = _make_lead_sites(
+            builder, interface_order)
         cell_size = len(lsites_with) + len(lsites_without)
 
 
-        sites = lsites_with + lsites_without + interface
+        fd_sites = lsites_with + lsites_without
+        sites = fd_sites + interface
         id_by_site = {}
         for site_id, site in enumerate(sites):
             id_by_site[site] = site_id
@@ -2967,25 +2828,24 @@ class InfiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.Infinite
         del id_by_site  # cleanup due to large size
 
         # In order to conform to the kwant.system.InfiniteVectorizedSystem
-        # interface we need to put the sites that connect to the previous
-        # cell *first*, then the sites that do not couple to the previous
-        # cell, then the sites in the previous cell. Because sites in
-        # a SiteArray are sorted by tag this means that the sites in these
-        # 3 different sets need to be in different SiteArrays.
+        # interface we need to put the interface sites (which have hoppings
+        # to the next unit cell) before non-interface sites.
+        # Note that we do not *store* the interface sites of the previous
+        # unit cell, but we still need to *handle* site indices >= cell_size
+        # because those are used e.g. in the graph.
         site_arrays = (
             _make_site_arrays(lsites_with)
             + _make_site_arrays(lsites_without)
-            + _make_site_arrays(interface)
         )
 
         (onsite_subgraphs, onsite_terms, onsite_term_values,
-         onsite_term_errors, _onsite_term_by_site_id) =\
-            _make_onsite_terms(builder, sites, site_arrays, term_offset=0)
+         onsite_term_errors, _onsite_term_by_site_id) = _make_onsite_terms(
+             builder, fd_sites, site_arrays, term_offset=0)
 
         (hopping_subgraphs, hopping_terms, hopping_term_values,
-         hopping_term_errors, _hopping_term_by_edge_id) =\
-            _make_hopping_terms(builder, graph, sites, site_arrays,
-                                cell_size, term_offset=len(onsite_terms))
+         hopping_term_errors, _hopping_term_by_edge_id) = _make_hopping_terms(
+             builder, graph, fd_sites, site_arrays, cell_size,
+             term_offset=len(onsite_terms))
 
         # Construct the combined onsite/hopping term datastructures
         subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
@@ -3001,8 +2861,12 @@ class InfiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.Infinite
         parameters = frozenset(parameters)
 
         self.site_arrays = site_arrays
-        self.sites = _Sites(self.site_arrays)
-        self.id_by_site = _IdBySite(self.site_arrays)
+        # 'sites' and 'id_by_site' have to be backwards compatible with the
+        # unvectorized interface, so we have to be able to index the interface
+        # sites in the previous unit cell also
+        _extended_site_arrays = self.site_arrays + _make_site_arrays(interface)
+        self.sites = _Sites(_extended_site_arrays)
+        self.id_by_site = _IdBySite(_extended_site_arrays)
         self.graph = graph
         self.subgraphs = subgraphs
         self.terms = terms
diff --git a/kwant/continuum/landau_levels.py b/kwant/continuum/landau_levels.py
index c653c55ef2634cd13e932f0c3f9c5fd026c10ff9..b2b8fc802dfa03343c4d7d314d84203c7c581302 100644
--- a/kwant/continuum/landau_levels.py
+++ b/kwant/continuum/landau_levels.py
@@ -18,6 +18,7 @@ import sympy
 
 import kwant.lattice
 import kwant.builder
+import kwant.system
 
 import kwant.continuum
 import kwant.continuum._common
@@ -144,7 +145,7 @@ def discretize_landau(hamiltonian, N, momenta=None, grid_spacing=1):
     if _has_coordinate(normal_coordinate, hamiltonian):
         sym = kwant.lattice.TranslationalSymmetry([grid_spacing, 0])
     else:
-        sym = kwant.builder.NoSymmetry()
+        sym = kwant.system.NoSymmetry()
     lat = LandauLattice(grid_spacing, norbs=norbs)
     syst = kwant.Builder(sym)
 
diff --git a/kwant/continuum/tests/test_landau_levels.py b/kwant/continuum/tests/test_landau_levels.py
index 8a568e6f4a4a716f1daa2c82542f47ea237a0f9d..5878724fbf997949d791dc33ae49a90e55b24a39 100644
--- a/kwant/continuum/tests/test_landau_levels.py
+++ b/kwant/continuum/tests/test_landau_levels.py
@@ -13,8 +13,8 @@ import sympy
 import pytest
 import itertools
 
-import kwant.builder
 import kwant.lattice
+import kwant.system
 
 from .._common import position_operators, momentum_operators, sympify
 from ..landau_levels import (
@@ -118,7 +118,7 @@ def test_discretize_landau():
     # test a basic Hamiltonian with no normal coordinate dependence
     syst = discretize_landau("k_x**2 + k_y**2", N=n_levels)
     lat = LandauLattice(1, norbs=1)
-    assert isinstance(syst.symmetry, kwant.builder.NoSymmetry)
+    assert isinstance(syst.symmetry, kwant.system.NoSymmetry)
     syst = syst.finalized()
     assert set(syst.sites) == {lat(0, j) for j in range(n_levels)}
     assert np.allclose(
diff --git a/kwant/lattice.py b/kwant/lattice.py
index e2557796f127c14ffde5a14ac7e2681a98ba7b60..976213d3578899a497e06da357aa512d8289f24d 100644
--- a/kwant/lattice.py
+++ b/kwant/lattice.py
@@ -153,7 +153,7 @@ class Polyatomic:
         algorithm finds and yields all the lattice sites inside the specified
         shape starting from the specified position.
 
-        A `~kwant.builder.Symmetry` or `~kwant.builder.Builder` may be passed as
+        A `~kwant.system.Symmetry` or `~kwant.builder.Builder` may be passed as
         sole argument when calling the function returned by this method.  This
         will restrict the flood-fill to the fundamental domain of the symmetry
         (or the builder's symmetry).  Note that unless the shape function has
@@ -174,8 +174,8 @@ class Polyatomic:
             Site = system.Site
 
             if symmetry is None:
-                symmetry = builder.NoSymmetry()
-            elif not isinstance(symmetry, builder.Symmetry):
+                symmetry = system.NoSymmetry()
+            elif not isinstance(symmetry, system.Symmetry):
                 symmetry = symmetry.symmetry
 
             def fd_site(lat, tag):
@@ -259,7 +259,7 @@ class Polyatomic:
         center = ta.array(center, float)
 
         def wire_sites(sym):
-            if not isinstance(sym, builder.Symmetry):
+            if not isinstance(sym, system.Symmetry):
                 sym = sym.symmetry
             if not isinstance(sym, TranslationalSymmetry):
                 raise ValueError('wire shape only works with '
@@ -514,7 +514,7 @@ class Monatomic(system.SiteFamily, Polyatomic):
 # The following class is designed such that it should avoid floating
 # point precision issues.
 
-class TranslationalSymmetry(builder.Symmetry):
+class TranslationalSymmetry(system.Symmetry):
     """A translational symmetry defined in real space.
 
     An alias exists for this common name: ``kwant.TranslationalSymmetry``.
@@ -568,7 +568,7 @@ class TranslationalSymmetry(builder.Symmetry):
         return TranslationalSymmetry(*ta.dot(generators, self.periods))
 
     def has_subgroup(self, other):
-        if isinstance(other, builder.NoSymmetry):
+        if isinstance(other, system.NoSymmetry):
             return True
         elif not isinstance(other, TranslationalSymmetry):
             raise ValueError("Unknown symmetry type.")
@@ -718,8 +718,9 @@ class TranslationalSymmetry(builder.Symmetry):
             raise ValueError("must provide a single group element when "
                              "acting on single sites.")
         if (len(element.shape) == 1 and not is_site):
-            raise ValueError("must provide a sequence of group elements "
-                             "when acting on site arrays.")
+            # We can act on a whole SiteArray with a single group element
+            # using numpy broadcasting.
+            element = element.reshape(1, -1)
         m_part = self._get_site_family_data(a.family)[0]
         try:
             delta = array_mod.dot(m_part, element)
diff --git a/kwant/physics/tests/test_gauge.py b/kwant/physics/tests/test_gauge.py
index 067c3244d44363ffcf6e9e4f993485b465e10e34..40e5f66d7560ad8c6f70111ac1dd73c7b89411bb 100644
--- a/kwant/physics/tests/test_gauge.py
+++ b/kwant/physics/tests/test_gauge.py
@@ -6,7 +6,8 @@ import pytest
 
 import kwant
 from ... import lattice
-from ...builder import HoppingKind, Builder, NoSymmetry, Site
+from ...builder import HoppingKind, Builder, Site
+from ...system import NoSymmetry
 from .. import gauge
 
 
diff --git a/kwant/system.py b/kwant/system.py
index 9674d8da1d1e5e850e6d823119597c3399db938e..4333e625980d1089f1bcfa03fec1ed114c7e294d 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -9,7 +9,7 @@
 """Low-level interface of systems"""
 
 __all__ = [
-    'Site', 'SiteArray', 'SiteFamily',
+    'Site', 'SiteArray', 'SiteFamily', 'Symmetry', 'NoSymmetry',
     'System', 'VectorizedSystem', 'FiniteSystem', 'FiniteVectorizedSystem',
     'InfiniteSystem', 'InfiniteVectorizedSystem',
     'is_finite', 'is_infinite', 'is_vectorized',
@@ -22,6 +22,7 @@ from copy import copy
 from collections import namedtuple
 from functools import total_ordering, lru_cache
 import numpy as np
+import tinyarray as ta
 from . import _system
 from ._common  import deprecate_args, KwantDeprecationWarning
 
@@ -269,6 +270,180 @@ class SiteFamily:
                              'site_family((1, 2))!')
         return Site(self, tag)
 
+
+################ Symmetries
+
+class Symmetry(metaclass=abc.ABCMeta):
+    """Abstract base class for spatial symmetries.
+
+    Many physical systems possess a discrete spatial symmetry, which results in
+    special properties of these systems.  This class is the standard way to
+    describe discrete spatial symmetries in Kwant.  An instance of this class
+    can be passed to a `Builder` instance at its creation.  The most important
+    kind of symmetry is translational symmetry, used to define scattering
+    leads.
+
+    Each symmetry has a fundamental domain -- a set of sites and hoppings,
+    generating all the possible sites and hoppings upon action of symmetry
+    group elements.  A class derived from `Symmetry` has to implement mapping
+    of any site or hopping into the fundamental domain, applying a symmetry
+    group element to a site or a hopping, and a method `which` to determine the
+    group element bringing some site from the fundamental domain to the
+    requested one.  Additionally, it has to have a property `num_directions`
+    returning the number of independent symmetry group generators (number of
+    elementary periods for translational symmetry).
+
+    A ``ValueError`` must be raised by the symmetry class whenever a symmetry
+    is used together with sites whose site family is not compatible with it.  A
+    typical example of this is when the vector defining a translational
+    symmetry is not a lattice vector.
+
+    The type of the domain objects as handled by the methods of this class is
+    not specified.  The only requirement is that it must support the unary
+    minus operation.  The reference implementation of `to_fd()` is hence
+    `self.act(-self.which(a), a, b)`.
+    """
+
+    @abc.abstractproperty
+    def num_directions(self):
+        """Number of elementary periods of the symmetry."""
+        pass
+
+    @abc.abstractmethod
+    def which(self, site):
+        """Calculate the domain of the site.
+
+        Parameters
+        ----------
+        site : `~kwant.system.Site` or `~kwant.system.SiteArray`
+
+        Returns
+        -------
+        group_element : tuple or sequence of tuples
+            A single tuple if ``site`` is a Site, or a sequence of tuples if
+            ``site`` is a SiteArray.  The group element(s) whose action
+            on a certain site(s) from the fundamental domain will result
+            in the given ``site``.
+        """
+        pass
+
+    @abc.abstractmethod
+    def act(self, element, a, b=None):
+        """Act with symmetry group element(s) on site(s) or hopping(s).
+
+        Parameters
+        ----------
+        element : tuple or sequence of tuples
+            Group element(s) with which to act on the provided site(s)
+            or hopping(s)
+        a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
+            If Site then ``element`` is a single tuple, if SiteArray then
+            ``element`` is a single tuple or a sequence of tuples.
+            If only ``a`` is provided then ``element`` acts on the site(s)
+            of ``a``. If ``b`` is also provided then ``element`` acts
+            on the hopping(s) ``(a, b)``.
+        """
+        pass
+
+    def to_fd(self, a, b=None):
+        """Map a site or hopping to the fundamental domain.
+
+        Parameters
+        ----------
+        a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
+
+        If ``b`` is None, return a site equivalent to ``a`` within the
+        fundamental domain.  Otherwise, return a hopping equivalent to ``(a,
+        b)`` but where the first element belongs to the fundamental domain.
+
+        Equivalent to `self.act(-self.which(a), a, b)`.
+        """
+        return self.act(-self.which(a), a, b)
+
+    def in_fd(self, site):
+        """Tell whether ``site`` lies within the fundamental domain.
+
+        Parameters
+        ----------
+        site : `~kwant.system.Site` or `~kwant.system.SiteArray`
+
+        Returns
+        -------
+        in_fd : bool or sequence of bool
+            single bool if ``site`` is a Site, or a sequence of
+            bool if ``site`` is a SiteArray. In the latter case
+            we return whether each site in the SiteArray is in
+            the fundamental domain.
+        """
+        if isinstance(site, Site):
+            for d in self.which(site):
+                if d != 0:
+                    return False
+            return True
+        elif isinstance(site, SiteArray):
+            which = self.which(site)
+            return np.logical_and.reduce(which != 0, axis=1)
+        else:
+            raise TypeError("'site' must be a Site or SiteArray")
+
+    @abc.abstractmethod
+    def subgroup(self, *generators):
+        """Return the subgroup generated by a sequence of group elements."""
+        pass
+
+    @abc.abstractmethod
+    def has_subgroup(self, other):
+        """Test whether `self` has the subgroup `other`...
+
+        or, in other words, whether `other` is a subgroup of `self`.  The
+        reason why this is the abstract method (and not `is_subgroup`) is that
+        in general it's not possible for a subgroup to know its supergroups.
+
+        """
+        pass
+
+
+class NoSymmetry(Symmetry):
+    """A symmetry with a trivial symmetry group."""
+
+    def __eq__(self, other):
+        return isinstance(other, NoSymmetry)
+
+    def __ne__(self, other):
+        return not self.__eq__(other)
+
+    def __repr__(self):
+        return 'NoSymmetry()'
+
+    @property
+    def num_directions(self):
+        return 0
+
+    periods = ()
+
+    _empty_array = ta.array((), int)
+
+    def which(self, site):
+        return self._empty_array
+
+    def act(self, element, a, b=None):
+        if element:
+            raise ValueError('`element` must be empty for NoSymmetry.')
+        return a if b is None else (a, b)
+
+    def to_fd(self, a, b=None):
+        return a if b is None else (a, b)
+
+    def in_fd(self, site):
+        return True
+
+    def subgroup(self, *generators):
+        if any(generators):
+            raise ValueError('Generators must be empty for NoSymmetry.')
+        return NoSymmetry(generators)
+
+    def has_subgroup(self, other):
+        return isinstance(other, NoSymmetry)
 
 
 ################ Systems
@@ -354,7 +529,7 @@ class System(metaclass=abc.ABCMeta):
 
 Term = namedtuple(
     "Term",
-    ["subgraph", "hermitian", "parameters"],
+    ["subgraph", "symmetry_element", "hermitian", "parameters"],
 )
 
 
@@ -363,6 +538,8 @@ class VectorizedSystem(System, metaclass=abc.ABCMeta):
 
     Attributes
     ----------
+    symmetry : kwant.system.Symmetry
+        The symmetry of the system.
     graph : kwant.graph.CGraph
         The system graph.
     subgraphs : sequence of tuples
@@ -371,9 +548,11 @@ class VectorizedSystem(System, metaclass=abc.ABCMeta):
         indexed by 'idx1' and 'idx2'.
     terms : sequence of tuples
         Each tuple has the following structure:
-        (subgraph: int, hermitian: bool, parameters: List(str))
+        (subgraph: int, symmetry_element: tuple, hermitian: bool, parameters: List(str))
         'subgraph' indexes 'subgraphs' and supplies the to/from sites of this
-        term. 'hermitian' is 'True' if the term needs its Hermitian
+        term. 'symmetry_element' is the symmetry group element that should be
+        applied to the 'to-sites' of this term.
+        'hermitian' is 'True' if the term needs its Hermitian
         conjugate to be added when evaluating the Hamiltonian, and 'parameters'
         contains a list of parameter names used when evaluating this term.
     site_arrays : sequence of SiteArray
@@ -576,50 +755,7 @@ def is_finite(syst):
 
 
 class InfiniteSystemMixin(metaclass=abc.ABCMeta):
-    """Abstract infinite low-level system.
-
-    An infinite system consists of an infinite series of identical cells.
-    Adjacent cells are connected by identical inter-cell hoppings.
-
-    Attributes
-    ----------
-    cell_size : integer
-        The number of sites in a single cell of the system.
-
-    Notes
-    -----
-    The system graph of an infinite systems contains a single cell, as well as
-    the part of the previous cell which is connected to it.  The first
-    `cell_size` sites form one complete single cell.  The remaining ``N`` sites
-    of the graph (``N`` equals ``graph.num_nodes - cell_size``) belong to the
-    previous cell.  They are included so that hoppings between cells can be
-    represented.  The N sites of the previous cell correspond to the first
-    ``N`` sites of the fully included cell.  When an ``InfiniteSystem`` is used
-    as a lead, ``N`` acts also as the number of interface sites to which it
-    must be connected.
-
-    The drawing shows three cells of an infinite system.  Each cell consists
-    of three sites.  Numbers denote sites which are included into the system
-    graph.  Stars denote sites which are not included.  Hoppings are included
-    in the graph if and only if they occur between two sites which are part of
-    the graph::
-
-            * 2 *
-        ... | | | ...
-            * 0 3
-            |/|/|
-            *-1-4
-
-        <-- order of cells
 
-    The numbering of sites in the drawing is one of the two valid ones for that
-    infinite system.  The other scheme has the numbers of site 0 and 1
-    exchanged, as well as of site 3 and 4.
-
-    Sites in the fundamental domain cell must belong to a different site array
-    than the sites in the previous cell. In the above example this means that
-    sites '(0, 1, 2)' and '(3, 4)' must belong to different site arrays.
-    """
     @deprecate_args
     def modes(self, energy=0, args=(), *, params=None):
         """Return mode decomposition of the lead
@@ -704,6 +840,46 @@ class InfiniteSystemMixin(metaclass=abc.ABCMeta):
 
 
 class InfiniteSystem(System, InfiniteSystemMixin, metaclass=abc.ABCMeta):
+    """Abstract infinite low-level system.
+
+    An infinite system consists of an infinite series of identical cells.
+    Adjacent cells are connected by identical inter-cell hoppings.
+
+    Attributes
+    ----------
+    cell_size : integer
+        The number of sites in a single cell of the system.
+
+    Notes
+    -----
+    The system graph of an infinite systems contains a single cell, as well as
+    the part of the previous cell which is connected to it.  The first
+    `cell_size` sites form one complete single cell.  The remaining ``N`` sites
+    of the graph (``N`` equals ``graph.num_nodes - cell_size``) belong to the
+    previous cell.  They are included so that hoppings between cells can be
+    represented.  The N sites of the previous cell correspond to the first
+    ``N`` sites of the fully included cell.  When an ``InfiniteSystem`` is used
+    as a lead, ``N`` acts also as the number of interface sites to which it
+    must be connected.
+
+    The drawing shows three cells of an infinite system.  Each cell consists
+    of three sites.  Numbers denote sites which are included into the system
+    graph.  Stars denote sites which are not included.  Hoppings are included
+    in the graph if and only if they occur between two sites which are part of
+    the graph::
+
+            * 2 *
+        ... | | | ...
+            * 0 3
+            |/|/|
+            *-1-4
+
+        <-- order of cells
+
+    The numbering of sites in the drawing is one of the two valid ones for that
+    infinite system.  The other scheme has the numbers of site 0 and 1
+    exchanged, as well as of site 3 and 4.
+    """
 
     @deprecate_args
     def cell_hamiltonian(self, args=(), sparse=False, *, params=None):
@@ -730,9 +906,46 @@ class InfiniteSystem(System, InfiniteSystemMixin, metaclass=abc.ABCMeta):
 
 
 class InfiniteVectorizedSystem(VectorizedSystem, InfiniteSystemMixin, metaclass=abc.ABCMeta):
+    """Abstract vectorized infinite low-level system.
+
+    An infinite system consists of an infinite series of identical cells.
+    Adjacent cells are connected by identical inter-cell hoppings.
+
+    Attributes
+    ----------
+    cell_size : integer
+        The number of sites in a single cell of the system.
+
+    Notes
+    -----
+    Unlike `~kwant.system.InfiniteSystem`, vectorized infinite systems do
+    not explicitly store the sites in the previous unit cell; only the
+    sites in the fundamental domain are stored. Nevertheless, the
+    SiteArrays of `~kwant.system.InfiniteVectorizedSystem` are ordered
+    in an analogous way, in order to facilitate the representation of
+    inter-cell hoppings. The ordering is as follows. The *interface sites*
+    of a unit cell are the sites that have hoppings to the *next* unit cell
+    (along the symmetry direction). Interface sites are always in different
+    SiteArrays than non-interface sites, i.e. the sites in a given SiteArray
+    are either all interface sites, or all non-interface sites.
+    The SiteArrays consisting of interface sites always appear *before* the
+    SiteArrays consisting of non-interface sites in ``self.site_arrays``.
+    This is backwards compatible with `kwant.system.InfiniteSystem`.
+
+    For backwards compatibility, `~kwant.system.InfiniteVectorizedSystem`
+    maintains a ``graph``, that includes nodes for the sites
+    in the previous unit cell.
+    """
     cell_hamiltonian = _system.vectorized_cell_hamiltonian
     inter_cell_hopping = _system.vectorized_inter_cell_hopping
 
+    def hamiltonian_submatrix(self, args=(), sparse=False,
+                              return_norb=False, *, params=None):
+        raise ValueError(
+            "'hamiltonian_submatrix' is not meaningful for infinite"
+            "systems. Use 'cell_hamiltonian' or 'inter_cell_hopping."
+        )
+
 
 def is_infinite(syst):
     return isinstance(syst, (InfiniteSystem, InfiniteVectorizedSystem))
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index 2ee6e1e1ec855767f29826e78bf4502be2d35e52..b4c26c5b78c569e157c76081592e3f133273fe4c 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -153,7 +153,7 @@ def test_site_families():
     assert fam1 < fam2  # string '1' is lexicographically less than '2'
 
 
-class VerySimpleSymmetry(builder.Symmetry):
+class VerySimpleSymmetry(system.Symmetry):
     def __init__(self, period):
         self.period = period
 
@@ -162,7 +162,7 @@ class VerySimpleSymmetry(builder.Symmetry):
         return 1
 
     def has_subgroup(self, other):
-        if isinstance(other, builder.NoSymmetry):
+        if isinstance(other, system.NoSymmetry):
             return True
         elif isinstance(other, VerySimpleSymmetry):
             return not other.period % self.period
@@ -321,10 +321,10 @@ def random_hopping_integral(rng):
 
 
 def check_onsite(fsyst, sites, subset=False, check_values=True):
-    vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))
+    vectorized = system.is_vectorized(fsyst)
 
     if vectorized:
-        site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])
+        site_offsets, _, _ = fsyst.site_ranges.transpose()
 
     freq = {}
     for node in range(fsyst.graph.num_nodes):
@@ -353,10 +353,10 @@ def check_onsite(fsyst, sites, subset=False, check_values=True):
 
 
 def check_hoppings(fsyst, hops):
-    vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))
+    vectorized = system.is_vectorized(fsyst)
 
     if vectorized:
-        site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])
+        site_offsets, _, _ = fsyst.site_ranges.transpose()
 
     assert fsyst.graph.num_edges == 2 * len(hops)
     for edge_id, edge in enumerate(fsyst.graph):
@@ -366,6 +366,12 @@ def check_hoppings(fsyst, hops):
 
         if vectorized:
             term = fsyst._hopping_term_by_edge_id[edge_id]
+            # Map sites in previous cell to fundamental domain; vectorized
+            # infinite systems only store sites in the FD. We already know
+            # that this hopping is between unit cells because this is encoded
+            # in the edge_id and term id.
+            if system.is_infinite(fsyst):
+                i, j = i % fsyst.cell_size, j % fsyst.cell_size
             if term < 0:  # Hermitian conjugate
                 assert (head, tail) in hops
             else:
@@ -464,7 +470,8 @@ def test_finalization(vectorize):
     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)))
+    lead = builder.Builder(kwant.TranslationalSymmetry((size, 0)),
+                           vectorize=vectorize)
     for site, value in lead_sites.items():
         shift = rng.randrange(-5, 6) * size
         site = site[0] + shift, site[1]
@@ -757,10 +764,6 @@ def test_vectorized_hamiltonian_evaluation():
     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(),
@@ -812,7 +815,7 @@ def test_vectorized_value_normalization():
 
 
 @pytest.mark.parametrize("sym", [
-    builder.NoSymmetry(),
+    system.NoSymmetry(),
     kwant.TranslationalSymmetry([-1]),
 ])
 def test_vectorized_requires_norbs(sym):
@@ -921,7 +924,7 @@ def test_fill():
     ## 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(),
+                      (system.NoSymmetry(),
                        lambda pos: ta.dot(pos, pos) < 17)]:
         cubic = kwant.lattice.general(ta.identity(3), norbs=1)
 
@@ -1642,7 +1645,7 @@ def test_subs():
 
     lat = kwant.lattice.chain(norbs=1)
 
-    def make_system(sym=kwant.builder.NoSymmetry(), n=3):
+    def make_system(sym=system.NoSymmetry(), n=3):
         syst = kwant.Builder(sym)
         syst[(lat(i) for i in range(n))] = onsite
         syst[lat.neighbors()] = hopping
diff --git a/kwant/tests/test_lattice.py b/kwant/tests/test_lattice.py
index 5e7a9976ed74ef881ece72ada167c48c39f8ff13..4850e45137076dbe05615225530e2d8b1e671eed 100644
--- a/kwant/tests/test_lattice.py
+++ b/kwant/tests/test_lattice.py
@@ -11,7 +11,7 @@ from math import sqrt
 import numpy as np
 import tinyarray as ta
 from pytest import raises
-from kwant import lattice, builder
+from kwant import lattice, builder, system
 from kwant._common import ensure_rng
 import pytest
 
@@ -285,7 +285,7 @@ def test_symmetry_has_subgroup():
     ## test whether actual subgroups are detected as such
     vecs = rng.randn(3, 3)
     sym1 = lattice.TranslationalSymmetry(*vecs)
-    ns = builder.NoSymmetry()
+    ns = system.NoSymmetry()
     assert ns.has_subgroup(ns)
     assert sym1.has_subgroup(sym1)
     assert sym1.has_subgroup(ns)