diff --git a/TODO b/TODO
index adca04a7efcbf2f7e7c0c6296a3c44e2e3016e18..7a00941b3893a06e7c0a9a58efbcb2265730c2d2 100644
--- a/TODO
+++ b/TODO
@@ -59,11 +59,11 @@ Roughly in order of importance.                                     -*-org-*-
 * Use sparse linear algebra to calculate bands
   However, SciPy's sparse eigenvalues don't seem to work well.
 
-* Allow attaching leads with further than nearest slice hoppings.
+* Allow attaching leads with further than nearest lead unit cell hoppings.
   The most easy way to do this is increasing the period of the lead.
   Alternatively, generalize modes and InfiniteSystem format.
 
-* In finalized leads, only keep the sites of the slice.
+* In finalized leads, only keep the sites of a single lead unit cell.
 
 * Add a test of kwant that verifies QHE conductance quantization
 
diff --git a/doc/source/whatsnew/1.0.rst b/doc/source/whatsnew/1.0.rst
index 071736d7e14f11d74272c9bed2aa8189c8b01c15..13c8f574719f0f9d722c0604db9c0da3b59be286 100644
--- a/doc/source/whatsnew/1.0.rst
+++ b/doc/source/whatsnew/1.0.rst
@@ -47,6 +47,8 @@ Some renames
 ------------
 * site groups are now called site families.  This affects all the names that
   used to contain "group" or "groups".
+* lead slices are now referred to as lead cells:  This affects all names that
+  used to contain "slice" or "slices" in the context of leads.
 * ``self_energy`` has been renamed to ``selfenergy`` in all cases, most notably
   in `kwant.physics.selfenergy`.
 * ``wave_func`` has been renamed to `~kwant.solvers.default.wave_function`,
diff --git a/kwant/builder.py b/kwant/builder.py
index 72084c9e400728502c4bc7d794336b25ead2f604..beabe120900a2ff6875d7602f373998164c726ff 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -412,17 +412,17 @@ class BuilderLead(Lead):
     builder : `Builder`
         The tight-binding system of a lead. It has to possess appropriate
         symmetry, and it may not contain hoppings between further than
-        neighboring lead slices.
+        neighboring lead unit cells.
     interface : sequence of `Site` instances
         Sequence of sites in the scattering region to which the lead is
         attached.
 
     Notes
     -----
-    The hopping from the scattering region to the lead is assumed to be
-    equal to the hopping from a lead slice to the next one in the direction of
-    the symmetry vector (i.e. the lead is 'leaving' the system and starts
-    with a hopping).
+    The hopping from the scattering region to the lead is assumed to be equal to
+    the hopping from a lead unit cell to the next one in the direction of the
+    symmetry vector (i.e. the lead is 'leaving' the system and starts with a
+    hopping).
 
     The given order of interface sites is preserved throughout finalization.
 
@@ -971,7 +971,7 @@ class Builder(object):
         ------
         ValueError
             If `lead_builder` does not have proper symmetry, has hoppings with
-            range of more than one slice, or if it is not completely
+            range of more than one lead unit cell, or if it is not completely
             interrupted by the system.
 
         Notes
@@ -986,8 +986,8 @@ class Builder(object):
             raise ValueError('Only builders with a 1D symmetry are allowed.')
         for hopping in lead_builder.hoppings():
             if not -1 <= sym.which(hopping[1])[0] <= 1:
-                msg = 'Hopping {0} connects non-neighboring slices. Only ' +\
-                      'nearest-slice hoppings are allowed ' +\
+                msg = 'Hopping {0} connects non-neighboring lead unit cells.' +\
+                      'Only nearest-cell hoppings are allowed ' +\
                       '(consider increasing the lead period).'
                 raise ValueError(msg.format(hopping))
         if not H:
@@ -1169,10 +1169,10 @@ class Builder(object):
             else:
                 # Tail is a fund. domain site not connected to prev. domain.
                 lsites_without.append(tail)
-        slice_size = len(lsites_with) + len(lsites_without)
+        cell_size = len(lsites_with) + len(lsites_without)
 
         if not lsites_with:
-            warnings.warn('Infinite system with disconnected slices.',
+            warnings.warn('Infinite system with disconnected cells.',
                           RuntimeWarning, stacklevel=3)
 
         ### Create list of sites and a lookup table
@@ -1200,7 +1200,7 @@ class Builder(object):
                     if (-sym.which(shifted_iface_site)[0] - 1,) != shift:
                         raise ValueError(
                             'The sites in interface_order do not all '
-                            'belong to the same lead slice.')
+                            'belong to the same lead cell.')
                     else:
                         raise ValueError('A site in interface_order is not an '
                                          'interface site:\n' + str(iface_site))
@@ -1223,7 +1223,7 @@ class Builder(object):
         g = graph.Graph()
         g.num_nodes = len(sites)  # Some sites could not appear in any edge.
         onsite_hamiltonians = []
-        for tail_id, tail in enumerate(sites[:slice_size]):
+        for tail_id, tail in enumerate(sites[:cell_size]):
             onsite_hamiltonians.append(self.H[tail][1])
             for head in self._out_neighbors(tail):
                 head_id = id_by_site.get(head)
@@ -1234,11 +1234,11 @@ class Builder(object):
                     # to this one has been added already or will be added.
                     fd = sym.which(head)[0]
                     if fd != 1:
-                        msg = 'Further-than-nearest-neighbor slices ' \
+                        msg = 'Further-than-nearest-neighbor cells ' \
                               'are connected by hopping\n{0}.'
                         raise ValueError(msg.format((tail, head)))
                     continue
-                if head_id >= slice_size:
+                if head_id >= cell_size:
                     # Head belongs to previous domain.  The edge added here
                     # correspond to one left out just above.
                     g.add_edge(head_id, tail_id)
@@ -1251,7 +1251,7 @@ class Builder(object):
         for tail_id, head_id in g:
             tail = sites[tail_id]
             head = sites[head_id]
-            if tail_id >= slice_size:
+            if tail_id >= cell_size:
                 # The tail belongs to the previous domain.  Find the
                 # corresponding hopping with the tail in the fund. domain.
                 tail, head = sym.to_fd(tail, head)
@@ -1259,7 +1259,7 @@ class Builder(object):
 
         #### Assemble and return result.
         result = InfiniteSystem()
-        result.slice_size = slice_size
+        result.cell_size = cell_size
         result.sites = sites
         result.graph = g
         result.hoppings = hoppings
@@ -1307,8 +1307,8 @@ class InfiniteSystem(system.InfiniteSystem):
     """Finalized infinite system, extracted from a `Builder`."""
     def hamiltonian(self, i, j, *args):
         if i == j:
-            if i >= self.slice_size:
-                i -= self.slice_size
+            if i >= self.cell_size:
+                i -= self.cell_size
             value = self.onsite_hamiltonians[i]
             if callable(value):
                 value = value(self.symmetry.to_fd(self.sites[i]), *args)
diff --git a/kwant/physics/dispersion.py b/kwant/physics/dispersion.py
index a951424a8649c677eeb5a12935d8ad36cdf30f9f..be9277ff4386de64eee2d8ac448517c79a633d2a 100644
--- a/kwant/physics/dispersion.py
+++ b/kwant/physics/dispersion.py
@@ -40,10 +40,10 @@ class Bands(object):
     """
 
     def __init__(self, sys, args=()):
-        self.ham = sys.slice_hamiltonian(args=args)
+        self.ham = sys.cell_hamiltonian(args=args)
         if not np.allclose(self.ham, self.ham.T.conj()):
-            raise ValueError('The slice Hamiltonian is not Hermitian.')
-        hop = sys.inter_slice_hopping(args=args)
+            raise ValueError('The cell Hamiltonian is not Hermitian.')
+        hop = sys.inter_cell_hopping(args=args)
         self.hop = np.empty(self.ham.shape, dtype=complex)
         self.hop[:, : hop.shape[1]] = hop
         self.hop[:, hop.shape[1]:] = 0
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index f251602809ecc6821f7503d4d3a2a1529bea8775..1326a0581b7ced1fd6c696d32f54042fb80248a3 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -21,16 +21,16 @@ __all__ = ['selfenergy', 'modes', 'ModesTuple', 'selfenergy_from_modes']
 
 Linsys = namedtuple('Linsys', ['eigenproblem', 'v', 'extract', 'project'])
 
-def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
+def setup_linsys(h_cell, h_hop, tol=1e6, algorithm=None):
     """
     Make an eigenvalue problem for eigenvectors of translation operator.
 
     Parameters
     ----------
-    h_onslice : numpy array with shape (n, n)
-        Hamiltonian of a single lead slice.
+    h_cell : numpy array with shape (n, n)
+        Hamiltonian of a single lead unit cell.
     h_hop : numpy array with shape (n, m), m <= n
-        Hopping Hamiltonian from the slice to the next one.
+        Hopping Hamiltonian from a cell to the next one.
     tol : float
         Numbers are considered zero when they are smaller than `tol` times
         the machine precision.
@@ -40,7 +40,7 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
         hopping svd, or lattice basis. If the real space basis is chosen, the
         following two options do not apply.
         The second value selects whether to add an anti-Hermitian term to the
-        slice Hamiltonian before inverting. Finally the third value selects
+        cell Hamiltonian before inverting. Finally the third value selects
         whether to reduce a generalized eigenvalue problem to the regular one.
         The default value, `None`, results in kwant selecting the algorithm
         that is the most efficient without sacrificing stability. Manual
@@ -61,22 +61,22 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
     The lead problem with degenerate hopping is rather complicated, and the
     details of the algorithm will be published elsewhere.
     """
-    n = h_onslice.shape[0]
+    n = h_cell.shape[0]
     m = h_hop.shape[1]
 
     if not (np.any(h_hop.real) or np.any(h_hop.imag)):
-        # Inter-slice hopping is zero.  The current algorithm is not suited to
+        # Inter-cell hopping is zero.  The current algorithm is not suited to
         # treat this extremely singular case.
         # Note: np.any(h_hop) returns (at least from numpy 1.6.*)
         #       False if h_hop is purely imaginary
-        raise ValueError("Inter-slice hopping is exactly zero.")
+        raise ValueError("Inter-cell hopping is exactly zero.")
 
     # If both h and t are real, it may be possible to use the real eigenproblem.
-    if (not np.any(h_hop.imag)) and (not np.any(h_onslice.imag)):
+    if (not np.any(h_hop.imag)) and (not np.any(h_cell.imag)):
         h_hop = h_hop.real
-        h_onslice = h_onslice.real
+        h_cell = h_cell.real
 
-    eps = np.finfo(np.common_type(h_onslice, h_hop)).eps * tol
+    eps = np.finfo(np.common_type(h_cell, h_hop)).eps * tol
 
     # First check if the hopping matrix has singular values close to 0.
     # (Close to zero is defined here as |x| < eps * tol * s[0] , where
@@ -92,8 +92,8 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
         # Hence the regular transfer matrix may be used.
         hop_inv = la.inv(h_hop)
 
-        A = np.zeros((2*n, 2*n), dtype=np.common_type(h_onslice, h_hop))
-        A[:n, :n] = dot(hop_inv, -h_onslice)
+        A = np.zeros((2*n, 2*n), dtype=np.common_type(h_cell, h_hop))
+        A[:n, :n] = dot(hop_inv, -h_cell)
         A[:n, n:] = -hop_inv
         A[n:, :n] = h_hop.T.conj()
 
@@ -140,9 +140,9 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
         # always if the original Hamiltonian is complex, and check for
         # invertibility first if it is real
 
-        h = h_onslice
+        h = h_cell
         sol = kla.lu_factor(h)
-        if issubclass(np.common_type(h_onslice, h_hop), np.floating) \
+        if issubclass(np.common_type(h_cell, h_hop), np.floating) \
            and need_to_stabilize is None:
             # Check if stabilization is needed.
             rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))
@@ -157,7 +157,7 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
             # Matrices are complex or need self-energy-like term to be
             # stabilized.
             temp = dot(u, u.T.conj()) + dot(v, v.T.conj())
-            h = h_onslice + 1j * temp
+            h = h_cell + 1j * temp
 
             sol = kla.lu_factor(h)
             rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))
@@ -358,7 +358,7 @@ def unified_eigenproblem(a, b=None, tol=1e6):
         An array of eigenvalues (can contain NaNs and Infs, but those
         are not accessed in `modes()`) The number of eigenvalues equals
         twice the number of nonzero singular values of
-        `h_hop` (so `2*h_onslice.shape[0]` if `h_hop` is invertible).
+        `h_hop` (so `2*h_cell.shape[0]` if `h_hop` is invertible).
     evanselect : numpy integer array
         Index array of right-decaying modes.
     propselect : numpy integer array
@@ -428,16 +428,16 @@ d.update(__doc__=modes_docstring)
 ModesTuple = type('ModesTuple', _Modes.__bases__, d)
 del _Modes, d, modes_docstring
 
-def modes(h_onslice, h_hop, tol=1e6, algorithm=None):
+def modes(h_cell, h_hop, tol=1e6, algorithm=None):
     """
     Compute the eigendecomposition of a translation operator of a lead.
 
     Parameters
     ----------
-    h_onslice : numpy array, real or complex, shape (N,N) The unit cell
-        Hamiltonian of the lead slice.
+    h_cell : numpy array, real or complex, shape (N,N) The unit cell
+        Hamiltonian of the lead unit cell.
     h_hop : numpy array, real or complex, shape (N,M)
-        The hopping matrix from a lead slice to the one on which self-energy
+        The hopping matrix from a lead cell to the one on which self-energy
         has to be calculated (and any other hopping in the same direction).
     tol : float
         Numbers and differences are considered zero when they are smaller
@@ -450,7 +450,7 @@ def modes(h_onslice, h_hop, tol=1e6, algorithm=None):
         mostly required for testing purposes.  The first value selects whether
         to work in the basis of the hopping svd, or lattice basis. If the real
         space basis is chosen, the following two options do not apply.  The
-        second value selects whether to add an anti-Hermitian term to the slice
+        second value selects whether to add an anti-Hermitian term to the cell
         Hamiltonian before inverting. Finally the third value selects whether
         to reduce a generalized eigenvalue problem to the regular one.
 
@@ -493,9 +493,9 @@ def modes(h_onslice, h_hop, tol=1e6, algorithm=None):
     """
     m = h_hop.shape[1]
 
-    if (h_onslice.shape[0] != h_onslice.shape[1] or
-        h_onslice.shape[0] != h_hop.shape[0]):
-        raise ValueError("Incompatible matrix sizes for h_onslice and h_hop.")
+    if (h_cell.shape[0] != h_cell.shape[1] or
+        h_cell.shape[0] != h_hop.shape[0]):
+        raise ValueError("Incompatible matrix sizes for h_cell and h_hop.")
 
     # Note: np.any(h_hop) returns (at least from numpy 1.6.1 - 1.8-devel)
     #       False if h_hop is purely imaginary
@@ -505,15 +505,14 @@ def modes(h_onslice, h_hop, tol=1e6, algorithm=None):
         return ModesTuple(np.empty((0, 0)), np.empty((0, 0)), 0, v)
 
     # Defer most of the calculation to helper routines.
-    matrices, v, extract, project = setup_linsys(h_onslice, h_hop, tol,
-                                                 algorithm)
+    matrices, v, extract, project = setup_linsys(h_cell, h_hop, tol, algorithm)
     ev, evanselect, propselect, vec_gen, ord_schur =\
          unified_eigenproblem(*(matrices + (tol,)))
 
     if v is not None:
         n = v.shape[1]
     else:
-        n = h_onslice.shape[0]
+        n = h_cell.shape[0]
 
     nprop = np.sum(propselect)
     nevan = n - nprop // 2
@@ -579,9 +578,9 @@ def selfenergy_from_modes(lead_modes):
     Returns
     -------
     Sigma : numpy array, real or complex, shape (M,M)
-        The computed self-energy. Note that even if `h_onslice` and `h_hop`
-        are both real, `Sigma` will typically be complex. (More precisely, if
-        there is a propagating mode, `Sigma` will definitely be complex.)
+        The computed self-energy. Note that even if `h_cell` and `h_hop` are
+        both real, `Sigma` will typically be complex. (More precisely, if there
+        is a propagating mode, `Sigma` will definitely be complex.)
 
     Notes
     -----
@@ -597,16 +596,16 @@ def selfenergy_from_modes(lead_modes):
         return la.solve(vecslmbdainv.T, vecs.T).T
 
 
-def selfenergy(h_onslice, h_hop, tol=1e6):
+def selfenergy(h_cell, h_hop, tol=1e6):
     """
     Compute the self-energy generated by the lead.
 
     Parameters
     ----------
-    h_onslice : numpy array, real or complex, shape (N,N) The unit cell
-        Hamiltonian of the lead slice.
+    h_cell : numpy array, real or complex, shape (N,N) The unit cell Hamiltonian
+        of the lead unit cell.
     h_hop : numpy array, real or complex, shape (N,M)
-        The hopping matrix from a lead slice to the one on which self-energy
+        The hopping matrix from a lead cell to the one on which self-energy
         has to be calculated (and any other hopping in the same direction).
     tol : float
         Numbers are considered zero when they are smaller than `tol` times
@@ -615,16 +614,16 @@ def selfenergy(h_onslice, h_hop, tol=1e6):
     Returns
     -------
     Sigma : numpy array, real or complex, shape (M,M)
-        The computed self-energy. Note that even if `h_onslice` and `h_hop`
-        are both real, `Sigma` will typically be complex. (More precisely, if
-        there is a propagating mode, `Sigma` will definitely be complex.)
+        The computed self-energy. Note that even if `h_cell` and `h_hop` are
+        both real, `Sigma` will typically be complex. (More precisely, if there
+        is a propagating mode, `Sigma` will definitely be complex.)
 
     Notes
     -----
     For simplicity this function internally calculates the modes first.
     This may cause a small slowdown, and can be improved if necessary.
     """
-    return selfenergy_from_modes(modes(h_onslice, h_hop, tol))
+    return selfenergy_from_modes(modes(h_cell, h_hop, tol))
 
 
 def square_selfenergy(width, hopping, fermi_energy):
@@ -644,7 +643,7 @@ def square_selfenergy(width, hopping, fermi_energy):
     # http://www.physik.uni-regensburg.de/forschung/\
     # richter/richter/media/research/publications2004/wimmer-Diplomarbeit.pdf
 
-    # p labels transversal modes.  i and j label the sites of a slice.
+    # p labels transversal modes.  i and j label the sites of a cell.
 
     # Precalculate the transverse wave function.
     psi_p_i = np.empty((width, width))
diff --git a/kwant/physics/tests/test_leads.py b/kwant/physics/tests/test_leads.py
index 88e0ea809f90ad9ecb2889ac4d449297d8a4d1b7..5568182edbdefabb7508b5fe442f0d9484283c81 100644
--- a/kwant/physics/tests/test_leads.py
+++ b/kwant/physics/tests/test_leads.py
@@ -15,7 +15,7 @@ import kwant
 
 modes_se = leads.selfenergy
 
-def h_slice(t, w, e):
+def h_cell_s_func(t, w, e):
     h = (4 * t - e) * np.identity(w)
     h.flat[1 :: w + 1] = -t
     h.flat[w :: w + 1] = -t
@@ -28,7 +28,7 @@ def test_analytic_numeric():
     e = 1.3                     # Fermi energy
 
     assert_almost_equal(leads.square_selfenergy(w, t, e),
-                        modes_se(h_slice(t, w, e), -t * np.identity(w)))
+                        modes_se(h_cell_s_func(t, w, e), -t * np.identity(w)))
 
 
 def test_regular_fully_degenerate():
@@ -39,21 +39,21 @@ def test_regular_fully_degenerate():
     e = 1.3                     # Fermi energy
 
     h_hop_s = -t * np.identity(w)
-    h_onslice_s = h_slice(t, w, e)
+    h_cell_s = h_cell_s_func(t, w, e)
 
     h_hop = np.zeros((2*w, 2*w))
     h_hop[:w, :w] = h_hop_s
     h_hop[w:, w:] = h_hop_s
 
-    h_onslice = np.zeros((2*w, 2*w))
-    h_onslice[:w, :w] = h_onslice_s
-    h_onslice[w:, w:] = h_onslice_s
+    h_cell = np.zeros((2*w, 2*w))
+    h_cell[:w, :w] = h_cell_s
+    h_cell[w:, w:] = h_cell_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
     g[:w, :w] = leads.square_selfenergy(w, t, e)
     g[w:, w:] = leads.square_selfenergy(w, t, e)
 
-    assert_almost_equal(g, modes_se(h_onslice, h_hop))
+    assert_almost_equal(g, modes_se(h_cell, h_hop))
 
 
 def test_regular_degenerate_with_crossing():
@@ -70,21 +70,21 @@ def test_regular_degenerate_with_crossing():
 
     global h_hop
     h_hop_s = -t * np.identity(w)
-    h_onslice_s = h_slice(t, w, e)
+    h_cell_s = h_cell_s_func(t, w, e)
 
     hop = np.zeros((2*w, 2*w))
     hop[:w, :w] = h_hop_s
     hop[w:, w:] = -h_hop_s
 
-    h_onslice = np.zeros((2*w, 2*w))
-    h_onslice[:w, :w] = h_onslice_s
-    h_onslice[w:, w:] = -h_onslice_s
+    h_cell = np.zeros((2*w, 2*w))
+    h_cell[:w, :w] = h_cell_s
+    h_cell[w:, w:] = -h_cell_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
     g[:w, :w] = leads.square_selfenergy(w, t, e)
     g[w:, w:] = -np.conj(leads.square_selfenergy(w, t, e))
 
-    assert_almost_equal(g, modes_se(h_onslice, hop))
+    assert_almost_equal(g, modes_se(h_cell, hop))
 
 
 def test_singular():
@@ -98,19 +98,19 @@ def test_singular():
     e = 0.4                     # Fermi energy
 
     h_hop_s = -t * np.identity(w)
-    h_onslice_s = h_slice(t, w, e)
+    h_cell_s = h_cell_s_func(t, w, e)
 
     h_hop = np.zeros((2*w, w))
     h_hop[w:, :w] = h_hop_s
 
-    h_onslice = np.zeros((2*w, 2*w))
-    h_onslice[:w, :w] = h_onslice_s
-    h_onslice[:w, w:] = h_hop_s
-    h_onslice[w:, :w] = h_hop_s
-    h_onslice[w:, w:] = h_onslice_s
+    h_cell = np.zeros((2*w, 2*w))
+    h_cell[:w, :w] = h_cell_s
+    h_cell[:w, w:] = h_hop_s
+    h_cell[w:, :w] = h_hop_s
+    h_cell[w:, w:] = h_cell_s
     g = leads.square_selfenergy(w, t, e)
 
-    assert_almost_equal(g, modes_se(h_onslice, h_hop))
+    assert_almost_equal(g, modes_se(h_cell, h_hop))
 
 
 def test_singular_but_square():
@@ -124,20 +124,20 @@ def test_singular_but_square():
     e = 2.38                     # Fermi energy
 
     h_hop_s = -t * np.identity(w)
-    h_onslice_s = h_slice(t, w, e)
+    h_cell_s = h_cell_s_func(t, w, e)
 
     h_hop = np.zeros((2*w, 2*w))
     h_hop[w:, :w] = h_hop_s
 
-    h_onslice = np.zeros((2*w, 2*w))
-    h_onslice[:w, :w] = h_onslice_s
-    h_onslice[:w, w:] = h_hop_s
-    h_onslice[w:, :w] = h_hop_s
-    h_onslice[w:, w:] = h_onslice_s
+    h_cell = np.zeros((2*w, 2*w))
+    h_cell[:w, :w] = h_cell_s
+    h_cell[:w, w:] = h_hop_s
+    h_cell[w:, :w] = h_hop_s
+    h_cell[w:, w:] = h_cell_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
     g[:w, :w] = leads.square_selfenergy(w, t, e)
-    assert_almost_equal(g, modes_se(h_onslice, h_hop))
+    assert_almost_equal(g, modes_se(h_cell, h_hop))
 
 
 def test_singular_fully_degenerate():
@@ -151,27 +151,27 @@ def test_singular_fully_degenerate():
     e = 3.3                     # Fermi energy
 
     h_hop_s = -t * np.identity(w)
-    h_onslice_s = h_slice(t, w, e)
+    h_cell_s = h_cell_s_func(t, w, e)
 
     h_hop = np.zeros((4*w, 2*w))
     h_hop[2*w:3*w, :w] = h_hop_s
     h_hop[3*w:4*w, w:2*w] = h_hop_s
 
-    h_onslice = np.zeros((4*w, 4*w))
-    h_onslice[:w, :w] = h_onslice_s
-    h_onslice[:w, 2*w:3*w] = h_hop_s
-    h_onslice[w:2*w, w:2*w] = h_onslice_s
-    h_onslice[w:2*w, 3*w:4*w] = h_hop_s
-    h_onslice[2*w:3*w, :w] = h_hop_s
-    h_onslice[2*w:3*w, 2*w:3*w] = h_onslice_s
-    h_onslice[3*w:4*w, w:2*w] = h_hop_s
-    h_onslice[3*w:4*w, 3*w:4*w] = h_onslice_s
+    h_cell = np.zeros((4*w, 4*w))
+    h_cell[:w, :w] = h_cell_s
+    h_cell[:w, 2*w:3*w] = h_hop_s
+    h_cell[w:2*w, w:2*w] = h_cell_s
+    h_cell[w:2*w, 3*w:4*w] = h_hop_s
+    h_cell[2*w:3*w, :w] = h_hop_s
+    h_cell[2*w:3*w, 2*w:3*w] = h_cell_s
+    h_cell[3*w:4*w, w:2*w] = h_hop_s
+    h_cell[3*w:4*w, 3*w:4*w] = h_cell_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
     g[:w, :w] = leads.square_selfenergy(w, t, e)
     g[w:, w:] = leads.square_selfenergy(w, t, e)
 
-    assert_almost_equal(g, modes_se(h_onslice, h_hop))
+    assert_almost_equal(g, modes_se(h_cell, h_hop))
 
 
 def test_singular_degenerate_with_crossing():
@@ -186,27 +186,27 @@ def test_singular_degenerate_with_crossing():
     e = 3.3                     # Fermi energy
 
     h_hop_s = -t * np.identity(w)
-    h_onslice_s = h_slice(t, w, e)
+    h_cell_s = h_cell_s_func(t, w, e)
 
     h_hop = np.zeros((4*w, 2*w))
     h_hop[2*w:3*w, :w] = h_hop_s
     h_hop[3*w:4*w, w:2*w] = -h_hop_s
 
-    h_onslice = np.zeros((4*w, 4*w))
-    h_onslice[:w, :w] = h_onslice_s
-    h_onslice[:w, 2*w:3*w] = h_hop_s
-    h_onslice[w:2*w, w:2*w] = -h_onslice_s
-    h_onslice[w:2*w, 3*w:4*w] = -h_hop_s
-    h_onslice[2*w:3*w, :w] = h_hop_s
-    h_onslice[2*w:3*w, 2*w:3*w] = h_onslice_s
-    h_onslice[3*w:4*w, w:2*w] = -h_hop_s
-    h_onslice[3*w:4*w, 3*w:4*w] = -h_onslice_s
+    h_cell = np.zeros((4*w, 4*w))
+    h_cell[:w, :w] = h_cell_s
+    h_cell[:w, 2*w:3*w] = h_hop_s
+    h_cell[w:2*w, w:2*w] = -h_cell_s
+    h_cell[w:2*w, 3*w:4*w] = -h_hop_s
+    h_cell[2*w:3*w, :w] = h_hop_s
+    h_cell[2*w:3*w, 2*w:3*w] = h_cell_s
+    h_cell[3*w:4*w, w:2*w] = -h_hop_s
+    h_cell[3*w:4*w, 3*w:4*w] = -h_cell_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
     g[:w, :w] = leads.square_selfenergy(w, t, e)
     g[w:, w:] = -np.conj(leads.square_selfenergy(w, t, e))
 
-    assert_almost_equal(g, modes_se(h_onslice, h_hop))
+    assert_almost_equal(g, modes_se(h_cell, h_hop))
 
 
 def test_singular_h_and_t():
@@ -235,7 +235,7 @@ def test_modes_bearded_ribbon():
                   (0, 0))] = 0.3
     sys[lat.neighbors()] = -1
     sys = sys.finalized()
-    h, t = sys.slice_hamiltonian(), sys.inter_slice_hopping()
+    h, t = sys.cell_hamiltonian(), sys.inter_cell_hopping()
     assert leads.modes(h, t).nmodes == 8  # Calculated by plotting dispersion.
 
 
diff --git a/kwant/plotter.py b/kwant/plotter.py
index 67e7b12dae468a66d569958c7e3feb872ea7d59a..233c5eab5716372cc73ab2797a2dcb0ed366c10a 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -718,14 +718,14 @@ def output_fig(fig, output_mode='auto', file=None, savefile_opts=None,
 
 # Extracting necessary data from the system.
 
-def sys_leads_sites(sys, num_lead_slices=2):
+def sys_leads_sites(sys, num_lead_cells=2):
     """Return all the sites of the system and of the leads as a list.
 
     Parameters
     ----------
     sys : kwant.builder.Builder or kwant.system.System instance
         The system, sites of which should be returned.
-    num_lead_slices : integer
+    num_lead_cells : integer
         The number of times lead sites from each lead should be returned.
         This is useful for showing several unit cells of the lead next to the
         system.
@@ -736,8 +736,8 @@ def sys_leads_sites(sys, num_lead_slices=2):
         A site is a `builder.Site` instance if the system is not finalized,
         and an integer otherwise.  For system sites `lead_number` is `None` and
         `copy_number` is `0`, for leads both are integers.
-    lead_slices : list of slices
-        `lead_slices[i]` gives the position of all the coordinates of lead
+    lead_cells : list of slices
+        `lead_cells[i]` gives the position of all the coordinates of lead
         `i` within `sites`.
 
     Notes
@@ -747,7 +747,7 @@ def sys_leads_sites(sys, num_lead_slices=2):
     unfinalized system, and sites of `system.InfiniteSystem` leads are
     returned with a finalized system.
     """
-    lead_slices = []
+    lead_cells = []
     if isinstance(sys, builder.Builder):
         sites = [(site, None, 0) for site in sys.sites()]
         for leadnr, lead in enumerate(sys.leads):
@@ -755,8 +755,8 @@ def sys_leads_sites(sys, num_lead_slices=2):
             if hasattr(lead, 'builder') and len(lead.interface):
                 sites.extend(((site, leadnr, i) for site in
                               lead.builder.sites() for i in
-                              xrange(num_lead_slices)))
-            lead_slices.append(slice(start, len(sites)))
+                              xrange(num_lead_cells)))
+            lead_cells.append(slice(start, len(sites)))
     elif isinstance(sys, system.FiniteSystem):
         sites = [(i, None, 0) for i in xrange(sys.graph.num_nodes)]
         for leadnr, lead in enumerate(sys.leads):
@@ -765,12 +765,12 @@ def sys_leads_sites(sys, num_lead_slices=2):
             if hasattr(lead, 'graph') and hasattr(lead, 'symmetry') and \
                     len(sys.lead_interfaces[leadnr]):
                 sites.extend(((site, leadnr, i) for site in
-                              xrange(lead.slice_size) for i in
-                              xrange(num_lead_slices)))
-            lead_slices.append(slice(start, len(sites)))
+                              xrange(lead.cell_size) for i in
+                              xrange(num_lead_cells)))
+            lead_cells.append(slice(start, len(sites)))
     else:
         raise TypeError('Unrecognized system type.')
-    return sites, lead_slices
+    return sites, lead_cells
 
 
 def sys_leads_pos(sys, site_lead_nr):
@@ -801,7 +801,7 @@ def sys_leads_pos(sys, site_lead_nr):
     # convert to a tuple and then to convert to numpy array ...
 
     is_builder = isinstance(sys, builder.Builder)
-    num_lead_slices = site_lead_nr[-1][2] + 1
+    num_lead_cells = site_lead_nr[-1][2] + 1
     if is_builder:
         pos = np.array(ta.array([i[0].pos for i in site_lead_nr]))
     else:
@@ -836,19 +836,19 @@ def sys_leads_pos(sys, site_lead_nr):
     vecs_doms = dict((i, get_vec_domain(i)) for i in xrange(len(sys.leads)))
     vecs_doms[None] = np.zeros((dim,)), 0
     for k, v in vecs_doms.iteritems():
-        vecs_doms[k] = [v[0] * i for i in xrange(v[1], v[1] + num_lead_slices)]
+        vecs_doms[k] = [v[0] * i for i in xrange(v[1], v[1] + num_lead_cells)]
     pos += [vecs_doms[i[1]][i[2]] for i in site_lead_nr]
     return pos
 
 
-def sys_leads_hoppings(sys, num_lead_slices=2):
+def sys_leads_hoppings(sys, num_lead_cells=2):
     """Return all the hoppings of the system and of the leads as an iterator.
 
     Parameters
     ----------
     sys : kwant.builder.Builder or kwant.system.System instance
         The system, sites of which should be returned.
-    num_lead_slices : integer
+    num_lead_cells : integer
         The number of times lead sites from each lead should be returned.
         This is useful for showing several unit cells of the lead next to the
         system.
@@ -859,8 +859,8 @@ def sys_leads_hoppings(sys, num_lead_slices=2):
         A site is a `builder.Site` instance if the system is not finalized,
         and an integer otherwise.  For system sites `lead_number` is `None` and
         `copy_number` is `0`, for leads both are integers.
-    lead_slices : list of slices
-        `lead_slices[i]` gives the position of all the coordinates of lead
+    lead_cells : list of slices
+        `lead_cells[i]` gives the position of all the coordinates of lead
         `i` within `hoppings`.
 
     Notes
@@ -871,7 +871,7 @@ def sys_leads_hoppings(sys, num_lead_slices=2):
     returned with a finalized system.
     """
     hoppings = []
-    lead_slices = []
+    lead_cells = []
     if isinstance(sys, builder.Builder):
         hoppings.extend(((hop, None, 0) for hop in sys.hoppings()))
 
@@ -893,8 +893,8 @@ def sys_leads_hoppings(sys, num_lead_slices=2):
             if hasattr(lead, 'builder') and len(lead.interface):
                 hoppings.extend(((hop, leadnr, i) for hop in
                                 lead_hoppings(lead.builder) for i in
-                                xrange(num_lead_slices)))
-            lead_slices.append(slice(start, len(hoppings)))
+                                xrange(num_lead_cells)))
+            lead_cells.append(slice(start, len(hoppings)))
     elif isinstance(sys, system.System):
         def ll_hoppings(sys):
             for i in xrange(sys.graph.num_nodes):
@@ -909,11 +909,11 @@ def sys_leads_hoppings(sys, num_lead_slices=2):
                     len(sys.lead_interfaces[leadnr]):
                 hoppings.extend(((hop, leadnr, i) for hop in
                                  ll_hoppings(lead) for i in
-                                 xrange(num_lead_slices)))
-            lead_slices.append(slice(start, len(hoppings)))
+                                 xrange(num_lead_cells)))
+            lead_cells.append(slice(start, len(hoppings)))
     else:
         raise TypeError('Unrecognized system type.')
-    return hoppings, lead_slices
+    return hoppings, lead_cells
 
 
 def sys_leads_hopping_pos(sys, hop_lead_nr):
@@ -942,7 +942,7 @@ def sys_leads_hopping_pos(sys, hop_lead_nr):
     is_builder = isinstance(sys, builder.Builder)
     if len(hop_lead_nr) == 0:
         return np.empty((0, 3)), np.empty((0, 3))
-    num_lead_slices = hop_lead_nr[-1][2] + 1
+    num_lead_cells = hop_lead_nr[-1][2] + 1
     if is_builder:
         pos = np.array(ta.array([ta.array(tuple(i[0][0].pos) +
                                           tuple(i[0][1].pos))
@@ -981,7 +981,7 @@ def sys_leads_hopping_pos(sys, hop_lead_nr):
     vecs_doms = dict((i, get_vec_domain(i)) for i in xrange(len(sys.leads)))
     vecs_doms[None] = np.zeros((dim,)), 0
     for k, v in vecs_doms.iteritems():
-        vecs_doms[k] = [v[0] * i for i in xrange(v[1], v[1] + num_lead_slices)]
+        vecs_doms[k] = [v[0] * i for i in xrange(v[1], v[1] + num_lead_cells)]
     pos += [vecs_doms[i[1]][i[2]] for i in hop_lead_nr]
     return np.copy(pos[:, : dim / 2]), np.copy(pos[:, dim / 2:])
 
@@ -998,7 +998,7 @@ _defaults = {'site_symbol': {2: 'o', 3: 'o'},
              'lead_color': {2: 'red', 3: 'red'}}
 
 
-def plot(sys, num_lead_slices=2, units='nn',
+def plot(sys, num_lead_cells=2, units='nn',
          site_symbol=None, site_size=None,
          site_color=None, site_edgecolor=None, site_lw=None,
          hop_color=None, hop_lw=None,
@@ -1014,7 +1014,7 @@ def plot(sys, num_lead_slices=2, units='nn',
     ----------
     sys : kwant.builder.Builder or kwant.system.FiniteSystem
         A system to be plotted.
-    num_lead_slices : int
+    num_lead_cells : int
         Number of lead copies to be shown with the system.
     units : 'pt', 'nn', or float
         The units in which symbol sizes, linewidths below are specified.
@@ -1080,7 +1080,7 @@ def plot(sys, num_lead_slices=2, units='nn',
     lead_site_size : number or `None`
         Relative (linear) size of the lead symbol
     lead_color : `matplotlib` color description or `None`
-        For the leads, `num_lead_slices` copies of the lead unit cell
+        For the leads, `num_lead_cells` copies of the lead unit cell
         are plotted. They are plotted in color fading from `lead_color`
         to white (alpha values in `lead_color` are supported) when moving
         from the system into the lead. Is also applied to the
@@ -1139,10 +1139,10 @@ def plot(sys, num_lead_slices=2, units='nn',
     """
 
     # Generate data.
-    sites, lead_sites_slcs = sys_leads_sites(sys, num_lead_slices)
+    sites, lead_sites_slcs = sys_leads_sites(sys, num_lead_cells)
     n_sys_sites = sum(i[1] is None for i in sites)
     sites_pos = sys_leads_pos(sys, sites)
-    hops, lead_hops_slcs = sys_leads_hoppings(sys, num_lead_slices)
+    hops, lead_hops_slcs = sys_leads_hoppings(sys, num_lead_cells)
     n_sys_hops = sum(i[1] is None for i in hops)
     end_pos, start_pos = sys_leads_hopping_pos(sys, hops)
 
@@ -1320,7 +1320,7 @@ def plot(sys, num_lead_slices=2, units='nn',
           linewidths=hop_lw, zorder=1, cmap=hop_cmap)
 
     # plot lead sites and hoppings
-    norm = matplotlib.colors.Normalize(-0.5, num_lead_slices - 0.5)
+    norm = matplotlib.colors.Normalize(-0.5, num_lead_cells - 0.5)
     lead_cmap = matplotlib.colors.LinearSegmentedColormap.from_list(None,
         [lead_color, (1, 1, 1, lead_color[3])])
 
@@ -1437,7 +1437,7 @@ def mask_interpolate(coords, values, a=None, method='nearest', oversampling=3):
 
 
 def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None,
-         a=None, method='nearest', oversampling=3, num_lead_slices=0,
+         a=None, method='nearest', oversampling=3, num_lead_cells=0,
         file=None, show=True,  dpi=None, fig_size=None):
     """Show interpolated map of a function defined for the sites of a system.
 
@@ -1472,7 +1472,7 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None,
         or "cubic"
     oversampling : integer, optional
         Number of pixels per reference length.  Defaults to 3.
-    num_lead_slices : integer, optional
+    num_lead_cells : integer, optional
         number of lead unit cells that should be plotted to indicate
         the position of leads. Defaults to 0.
     file : string or file object or `None`
@@ -1520,8 +1520,8 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None,
     image = ax.imshow(img.T, extent=(min[0], max[0], min[1], max[1]),
                       origin='lower', interpolation='none', cmap=cmap,
                       vmin=vmin, vmax=vmax)
-    if num_lead_slices:
-        plot(sys, num_lead_slices, site_symbol='no symbol', hop_lw=0,
+    if num_lead_cells:
+        plot(sys, num_lead_cells, site_symbol='no symbol', hop_lw=0,
              lead_site_symbol='s', lead_site_size=0.501,
              lead_site_lw=0, lead_hop_lw=0, lead_color='black',
              colorbar=False, ax=ax)
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index 2f3b654ddd028416027f34de2c0867f499a7da1a..3bc291d274ba43f3bb8fea36feea74c5218f9609 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -189,7 +189,7 @@ class SparseSolver(object):
                 u_in, ulinv_in = u[:, : nprop], ulinv[:, : nprop]
 
                 # Construct a matrix of 1's that translates the
-                # inter-slice hopping to a proper hopping
+                # inter-cell hopping to a proper hopping
                 # from the system to the lead.
                 iface_orbs = np.r_[tuple(slice(offsets[i], offsets[i + 1])
                                         for i in interface)]
diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py
index d6f61f7f2ac27e38e9bc477057bf6cf589d35d90..e16a2655bccddefacda82f2ce30855879f2623f7 100644
--- a/kwant/solvers/tests/_test_sparse.py
+++ b/kwant/solvers/tests/_test_sparse.py
@@ -59,11 +59,11 @@ def test_output(solve):
                         np.identity(s.shape[0]))
     assert_raises(ValueError, solve, fsys, 0, [])
     modes = solve(fsys).lead_info
-    h = fsys.leads[0].slice_hamiltonian()
-    t = fsys.leads[0].inter_slice_hopping()
+    h = fsys.leads[0].cell_hamiltonian()
+    t = fsys.leads[0].inter_cell_hopping()
     modes1 = kwant.physics.modes(h, t)
-    h = fsys.leads[1].slice_hamiltonian()
-    t = fsys.leads[1].inter_slice_hopping()
+    h = fsys.leads[1].cell_hamiltonian()
+    t = fsys.leads[1].inter_cell_hopping()
     modes2 = kwant.physics.modes(h, t)
     assert_almost_equal(modes1[0], modes[0][0])
     assert_almost_equal(modes2[1], modes[1][1])
@@ -239,7 +239,7 @@ def test_singular_graph_system(solve):
     assert_almost_equal(t_elements, t_el_should_be)
 
 
-# This test features inside the onslice Hamiltonian a hopping matrix with more
+# This test features inside the cell Hamiltonian a hopping matrix with more
 # zero eigenvalues than the lead hopping matrix. Older version of the
 # sparse solver failed here.
 def test_tricky_singular_hopping(solve):
diff --git a/kwant/system.py b/kwant/system.py
index e0843969bb32338926fee789212f415ee1f3bbaa..c6843a993c66bda409c68e20d90e0442312992a8 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -138,27 +138,27 @@ class InfiniteSystem(System):
     """
     Abstract infinite low-level system.
 
-    An infinite system consists of an infinite series of identical slices.
-    Adjacent slices are connected by identical inter-slice hoppings.
+    An infinite system consists of an infinite series of identical cells.
+    Adjacent cells are connected by identical inter-cell hoppings.
 
     Instance Variables
     ------------------
-    slice_size : integer
-        The number of sites in a single slice of the system.
+    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 slice, as well as
-    the part of the previous slice which is connected to it.  The first
-    `slice_size` sites form one complete single slice.  The remaining `N` sites
-    of the graph (`N` equals ``graph.num_nodes - slice_size``) belong to the
-    previous slice.  They are included so that hoppings between slices can be
-    represented.  The N sites of the previous slice correspond to the first `N`
-    sites of the fully included slice.  When an InfiniteSystem is used as a
-    lead, `N` acts also as the number of interface sites to which it must be
+    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 slices of an infinite system.  Each slice consists
+    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
@@ -170,7 +170,7 @@ class InfiniteSystem(System):
             |/|/|
             *-1-4
 
-        <-- order of slices
+        <-- 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
@@ -178,17 +178,17 @@ class InfiniteSystem(System):
     """
     __metaclass__ = abc.ABCMeta
 
-    def slice_hamiltonian(self, sparse=False, args=()):
-        """Hamiltonian of a single slice of the infinite system."""
-        slice_sites = xrange(self.slice_size)
-        return self.hamiltonian_submatrix(slice_sites, slice_sites,
+    def cell_hamiltonian(self, sparse=False, args=()):
+        """Hamiltonian of a single cell of the infinite system."""
+        cell_sites = xrange(self.cell_size)
+        return self.hamiltonian_submatrix(cell_sites, cell_sites,
                                           sparse=sparse, args=args)
 
-    def inter_slice_hopping(self, sparse=False, args=()):
-        """Hopping Hamiltonian between two slices of the infinite system."""
-        slice_sites = xrange(self.slice_size)
-        neighbor_sites = xrange(self.slice_size, self.graph.num_nodes)
-        return self.hamiltonian_submatrix(slice_sites, neighbor_sites,
+    def inter_cell_hopping(self, sparse=False, args=()):
+        """Hopping Hamiltonian between two cells of the infinite system."""
+        cell_sites = xrange(self.cell_size)
+        neighbor_sites = xrange(self.cell_size, self.graph.num_nodes)
+        return self.hamiltonian_submatrix(cell_sites, neighbor_sites,
                                           sparse=sparse, args=args)
 
     def modes(self, energy=0, args=()):
@@ -197,28 +197,28 @@ class InfiniteSystem(System):
         See documentation of `~kwant.physics.ModesTuple` for the return
         format details.
         """
-        ham = self.slice_hamiltonian(args=args)
+        ham = self.cell_hamiltonian(args=args)
         shape = ham.shape
         assert len(shape) == 2
         assert shape[0] == shape[1]
         # Subtract energy from the diagonal.
         ham.flat[::ham.shape[0] + 1] -= energy
-        return physics.modes(ham, self.inter_slice_hopping(args=args))
+        return physics.modes(ham, self.inter_cell_hopping(args=args))
 
     def selfenergy(self, energy=0, args=()):
         """Return self-energy of a lead.
 
         The returned matrix has the shape (s, s), where s is
         ``sum(len(self.hamiltonian(i, i))
-              for i in range(self.graph.num_nodes - self.slice_size))``.
+              for i in range(self.graph.num_nodes - self.cell_size))``.
         """
-        ham = self.slice_hamiltonian(args=args)
+        ham = self.cell_hamiltonian(args=args)
         shape = ham.shape
         assert len(shape) == 2
         assert shape[0] == shape[1]
         # Subtract energy from the diagonal.
         ham.flat[::ham.shape[0] + 1] -= energy
-        return physics.selfenergy(ham, self.inter_slice_hopping(args=args))
+        return physics.selfenergy(ham, self.inter_cell_hopping(args=args))
 
 
 class PrecalculatedLead(object):
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index f1f8a66471dfec26063bc3e977cf3b71a4230c61..447230954f225aaaa2834e7380be2b871b0f48b3 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -323,7 +323,7 @@ def test_finalization():
     assert_equal([fsys.site(i).tag for i in fsys.lead_interfaces[0]],
                  neighbors)
 
-    # Add a hopping to the lead which couples two next-nearest slices and check
+    # Add a hopping to the lead which couples two next-nearest cells and check
     # whether this leads to an error.
     a = rng.choice(lead_sites_list)
     b = rng.choice(possible_neighbors)