diff --git a/doc/source/images/ab_ring.py.diff b/doc/source/images/ab_ring.py.diff
index da410f8559cc46ad4400f0da7cf2ea930771fac6..7f09a4b46989b118ecc31b84ad89e348445fedc9 100644
--- a/doc/source/images/ab_ring.py.diff
+++ b/doc/source/images/ab_ring.py.diff
@@ -8,9 +8,9 @@
  
  
  def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
-@@ -38,12 +39,13 @@
-     for hopping in lat.nearest:
-         sys[sys.possible_hoppings(*hopping)] = -t
+@@ -37,12 +38,13 @@
+     sys[lat.shape(ring, (0, r1 + 1))] = 4 * t
+     sys[lat.nearest] = -t
  
 -    # In order to introduce a flux through the ring, we introduce a phase on
 -    # the hoppings on the line cut through one of the arms.  Since we want to
@@ -28,7 +28,7 @@
          return exp(1j * phi)
  
      def crosses_branchcut(hop):
-@@ -81,6 +83,54 @@
+@@ -82,6 +84,50 @@
      return sys
  
  
@@ -40,16 +40,14 @@
 +        rsq = x**2 + y**2
 +        return ( r1**2 < rsq < r2**2)
 +    sys[lat.shape(ring, (0, 11))] = 4 * t
-+    for hopping in lat.nearest:
-+        sys[sys.possible_hoppings(*hopping)] = -t
++    sys[lat.nearest] = -t
 +    sym_lead0 = kwant.TranslationalSymmetry((-a, 0))
 +    lead0 = kwant.Builder(sym_lead0)
 +    def lead_shape(pos):
 +        (x, y) = pos
 +        return (-1 < x < 1) and ( 0.5 * W < y < 1.5 * W )
 +    lead0[lat.shape(lead_shape, (0, W))] = 4 * t
-+    for hopping in lat.nearest:
-+        lead0[lead0.possible_hoppings(*hopping)] = -t
++    lead0[lat.nearest] = -t
 +    lead1 = lead0.reversed()
 +    sys.attach_lead(lead0)
 +    sys.attach_lead(lead1)
@@ -64,16 +62,14 @@
 +        rsq = x**2 + y**2
 +        return ( r1**2 < rsq < r2**2)
 +    sys[lat.shape(ring, (0, 11))] = 4 * t
-+    for hopping in lat.nearest:
-+        sys[sys.possible_hoppings(*hopping)] = -t
++    sys[lat.nearest] = -t
 +    sym_lead0 = kwant.TranslationalSymmetry((-a, 0))
 +    lead0 = kwant.Builder(sym_lead0)
 +    def lead_shape(pos):
 +        (x, y) = pos
 +        return (-1 < x < 1) and ( -W/2 < y < W/2  )
 +    lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
-+    for hopping in lat.nearest:
-+        lead0[lead0.possible_hoppings(*hopping)] = -t
++    lead0[lat.nearest] = -t
 +    lead1 = lead0.reversed()
 +    sys.attach_lead(lead0)
 +    sys.attach_lead(lead1, lat(0, 0))
@@ -83,7 +79,7 @@
  def plot_conductance(sys, energy, fluxes):
      # compute conductance
  
-@@ -90,18 +140,29 @@
+@@ -91,18 +137,29 @@
          smatrix = kwant.solve(sys, energy, kwargs={'phi': flux})
          data.append(smatrix.transmission(1, 0))
  
@@ -118,7 +114,7 @@
  
      # Finalize the system.
      sys = sys.finalized()
-@@ -111,6 +172,15 @@
+@@ -112,6 +169,15 @@
                                                  for i in xrange(100)])
  
  
diff --git a/doc/source/images/graphene.py.diff b/doc/source/images/graphene.py.diff
index facca27ff253f624db2b625366d3771b5a08997f..8443e7d186930bc2f4d3e52399652e4297d29cfa 100644
--- a/doc/source/images/graphene.py.diff
+++ b/doc/source/images/graphene.py.diff
@@ -8,16 +8,7 @@
  
  
  # Define the graphene lattice
-@@ -63,7 +64,7 @@
-         return (-1 < x < 1) and (-0.4 * r < y < 0.4 * r)
- 
-     lead0 = kwant.Builder(sym0)
--    lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
-+    lead0[graphene.shape(lead0_shape, (0, 0))] = - pot
-     for hopping in hoppings:
-         lead0[lead0.possible_hoppings(*hopping)] = -1
- 
-@@ -105,22 +106,40 @@
+@@ -102,22 +103,40 @@
          smatrix = kwant.solve(sys, energy)
          data.append(smatrix.transmission(0, 1))
  
@@ -66,7 +57,7 @@
  
  
  def main():
-@@ -133,17 +152,22 @@
+@@ -130,17 +149,22 @@
          return 0 if site.group == a else 1
  
      # Plot the closed system without leads.
diff --git a/doc/source/images/quantum_well.py.diff b/doc/source/images/quantum_well.py.diff
index d1f21475e83d17964b7f0f0a1b0c128d2852cca4..e5b45b0d3d2ab5dbfd22104dde3831b4f3a9731f 100644
--- a/doc/source/images/quantum_well.py.diff
+++ b/doc/source/images/quantum_well.py.diff
@@ -8,7 +8,7 @@
  
  
  def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
-@@ -63,19 +64,25 @@
+@@ -61,19 +62,25 @@
          smatrix = kwant.solve(sys, energy, kwargs={'pot': -welldepth})
          data.append(smatrix.transmission(1, 0))
  
diff --git a/doc/source/tutorial/ab_ring.py b/doc/source/tutorial/ab_ring.py
index 5c39ce688210b82992a76b384152ed346e890d45..e2ec01d61925ae076cfeea79d8f2a32205bb9b82 100644
--- a/doc/source/tutorial/ab_ring.py
+++ b/doc/source/tutorial/ab_ring.py
@@ -38,8 +38,7 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     # and add the corresponding lattice points using the `shape`-function
 #HIDDEN_BEGIN_lcak
     sys[lat.shape(ring, (0, r1 + 1))] = 4 * t
-    for hopping in lat.nearest:
-        sys[sys.possible_hoppings(*hopping)] = -t
+    sys[lat.nearest] = -t
 #HIDDEN_END_lcak
 
     # In order to introduce a flux through the ring, we introduce a phase on
@@ -54,13 +53,16 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     def crosses_branchcut(hop):
         ix0, iy0 = hop[0].tag
 
-        # possible_hoppings with the argument (1, 0) below
+        # builder.HoppingKind with the argument (1, 0) below
         # returns hoppings ordered as ((i+1, j), (i, j))
         return iy0 < 0 and ix0 == 1  # ix1 == 0 then implied
 
     # Modify only those hopings in x-direction that cross the branch cut
-    sys[(hop for hop in sys.possible_hoppings((1, 0), lat, lat)
-         if crosses_branchcut(hop))] = fluxphase
+    def hops_across_cut(sys):
+        for hop in kwant.builder.HoppingKind((1, 0), lat, lat)(sys):
+            if crosses_branchcut(hop):
+                yield hop
+    sys[hops_across_cut] = fluxphase
 #HIDDEN_END_lvkt
 
     #### Define the leads. ####
@@ -74,8 +76,7 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
         return (-1 < x < 1) and (-W / 2 < y < W / 2)
 
     lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
-    for hopping in lat.nearest:
-        lead0[lead0.possible_hoppings(*hopping)] = -t
+    lead0[lat.nearest] = -t
 #HIDDEN_END_qwgr
 
     # Then the lead to the right
diff --git a/doc/source/tutorial/closed_system.py b/doc/source/tutorial/closed_system.py
index d67c28992a2df393cf116a3ed7f32573399e83a4..a7bf7d2d9f48a224220884ffd19a3c10c7de29bc 100644
--- a/doc/source/tutorial/closed_system.py
+++ b/doc/source/tutorial/closed_system.py
@@ -43,9 +43,9 @@ def make_system(a=1, t=1.0, r=10):
 
     sys[lat.shape(circle, (0, 0))] = 4 * t
     # hoppings in x-direction
-    sys[sys.possible_hoppings((1, 0), lat, lat)] = hopx
+    sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
     # hoppings in y-directions
-    sys[sys.possible_hoppings((0, 1), lat, lat)] = -t
+    sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
 
     # It's a closed system for a change, so no leads
     return sys
diff --git a/doc/source/tutorial/graphene.py b/doc/source/tutorial/graphene.py
index 201a380c4bd7362cbcb27ca10a1b9c0d52fbe16f..4d12e1d76fba80ff73a27adbb965bec40ab55c49 100644
--- a/doc/source/tutorial/graphene.py
+++ b/doc/source/tutorial/graphene.py
@@ -49,13 +49,12 @@ def make_system(r=10, w=2.0, pot=0.1):
 #HIDDEN_END_shzy
 
     # specify the hoppings of the graphene lattice in the
-    # format expected by possibe_hoppings()
+    # format expected by builder.HoppingKind
 #HIDDEN_BEGIN_hsmc
     hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
 #HIDDEN_END_hsmc
 #HIDDEN_BEGIN_bfwb
-    for hopping in hoppings:
-        sys[sys.possible_hoppings(*hopping)] = -1
+    sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
 #HIDDEN_END_bfwb
 
     # Modify the scattering region
@@ -75,8 +74,7 @@ def make_system(r=10, w=2.0, pot=0.1):
 
     lead0 = kwant.Builder(sym0)
     lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
-    for hopping in hoppings:
-        lead0[lead0.possible_hoppings(*hopping)] = -1
+    lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
 
     # The second lead, going to the top right
     sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))
@@ -89,8 +87,7 @@ def make_system(r=10, w=2.0, pot=0.1):
 
     lead1 = kwant.Builder(sym1)
     lead1[graphene.shape(lead1_shape, (0, 0))] = pot
-    for hopping in hoppings:
-        lead1[lead1.possible_hoppings(*hopping)] = -1
+    lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
 #HIDDEN_END_aakh
 
 #HIDDEN_BEGIN_kmmw
diff --git a/doc/source/tutorial/quantum_well.py b/doc/source/tutorial/quantum_well.py
index f9563c8f38147a704ac974cc7faa5166fcbc8d60..c4d80c26a038b25088e0891730653a816f900b5e 100644
--- a/doc/source/tutorial/quantum_well.py
+++ b/doc/source/tutorial/quantum_well.py
@@ -35,8 +35,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
         return 4 * t + potential(site, pot)
 
     sys[(lat(x, y) for x in range(L) for y in range(W))] = onsite
-    for hopping in lat.nearest:
-        sys[sys.possible_hoppings(*hopping)] = -t
+    sys[lat.nearest] = -t
 #HIDDEN_END_coid
 
     #### Define the leads. ####
@@ -45,8 +44,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     lead0 = kwant.Builder(sym_lead0)
 
     lead0[(lat(0, j) for j in xrange(W))] = 4 * t
-    for hopping in lat.nearest:
-        lead0[lead0.possible_hoppings(*hopping)] = -t
+    lead0[lat.nearest] = -t
 
     # ... then the lead to the right.  We use a method that returns a copy of
     # `lead0` with its direction reversed.
diff --git a/doc/source/tutorial/quantum_wire_revisited.py b/doc/source/tutorial/quantum_wire_revisited.py
index 920ded17d3d3963903224fcd9cc3d2fb450b8da5..9ce7abcd256165896d1a14529cac55aaf68bf860 100644
--- a/doc/source/tutorial/quantum_wire_revisited.py
+++ b/doc/source/tutorial/quantum_wire_revisited.py
@@ -4,7 +4,7 @@
 #
 # Kwant features highlighted
 # --------------------------
-#  - Using iterables and possible_hoppings() for making systems
+#  - Using iterables and builder.HoppingKind for making systems
 #  - introducing `reversed()` for the leads
 #
 # Note: Does the same as tutorial1a.py, but using other features of kwant
@@ -29,8 +29,7 @@ def make_system(a=1, t=1.0, W=10, L=30):
     sys[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
 #HIDDEN_END_vvjt
 #HIDDEN_BEGIN_nooi
-    for hopping in lat.nearest:
-        sys[sys.possible_hoppings(*hopping)] = -t
+    sys[lat.nearest] = -t
 #HIDDEN_END_nooi
 
     #### Define the leads. ####
@@ -41,8 +40,7 @@ def make_system(a=1, t=1.0, W=10, L=30):
     lead0 = kwant.Builder(sym_lead0)
 
     lead0[(lat(0, j) for j in xrange(W))] = 4 * t
-    for hopping in lat.nearest:
-        lead0[lead0.possible_hoppings(*hopping)] = -t
+    lead0[lat.nearest] = -t
 #HIDDEN_END_iepx
 
     # ... then the lead to the right.  We use a method that returns a copy of
diff --git a/doc/source/tutorial/spin_orbit.py b/doc/source/tutorial/spin_orbit.py
index 55cce77809326bd24a78bfdef1eb744735635b2f..1d14331af4ab491a580b2b270e3e1b4d28cff0d2 100644
--- a/doc/source/tutorial/spin_orbit.py
+++ b/doc/source/tutorial/spin_orbit.py
@@ -41,10 +41,10 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     sys[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t * sigma_0 + \
         e_z * sigma_z
     # hoppings in x-direction
-    sys[sys.possible_hoppings((1, 0), lat, lat)] = -t * sigma_0 - \
+    sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = -t * sigma_0 - \
         1j * alpha * sigma_y
     # hoppings in y-directions
-    sys[sys.possible_hoppings((0, 1), lat, lat)] = -t * sigma_0 + \
+    sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t * sigma_0 + \
         1j * alpha * sigma_x
 #HIDDEN_END_uxrm
 
@@ -56,10 +56,10 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
 #HIDDEN_BEGIN_yliu
     lead0[(lat(0, j) for j in xrange(W))] = 4 * t * sigma_0 + e_z * sigma_z
     # hoppings in x-direction
-    lead0[lead0.possible_hoppings((1, 0), lat, lat)] = -t * sigma_0 - \
+    lead0[kwant.builder.HoppingKind((1, 0), lat, lat)] = -t * sigma_0 - \
         1j * alpha * sigma_y
     # hoppings in y-directions
-    lead0[lead0.possible_hoppings((0, 1), lat, lat)] = -t * sigma_0 + \
+    lead0[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t * sigma_0 + \
         1j * alpha * sigma_x
 #HIDDEN_END_yliu
 
diff --git a/doc/source/tutorial/superconductor_transport.py b/doc/source/tutorial/superconductor_transport.py
index ba7743f945751731b10d47069a0c7677d38af24e..55d2e6d0c3ae34a151b2d670175318074838b82f 100644
--- a/doc/source/tutorial/superconductor_transport.py
+++ b/doc/source/tutorial/superconductor_transport.py
@@ -36,10 +36,10 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
          for y in range(W))] = mu - 4 * t - barrier
 
     # hoppings in x and y-directions, for both electrons and holes
-    sys[sys.possible_hoppings((1, 0), lat_e, lat_e)] = -t
-    sys[sys.possible_hoppings((0, 1), lat_e, lat_e)] = -t
-    sys[sys.possible_hoppings((1, 0), lat_h, lat_h)] = t
-    sys[sys.possible_hoppings((0, 1), lat_h, lat_h)] = t
+    sys[kwant.builder.HoppingKind((1, 0), lat_e, lat_e)] = -t
+    sys[kwant.builder.HoppingKind((0, 1), lat_e, lat_e)] = -t
+    sys[kwant.builder.HoppingKind((1, 0), lat_h, lat_h)] = t
+    sys[kwant.builder.HoppingKind((0, 1), lat_h, lat_h)] = t
 
     # Superconducting order parameter enters as hopping between
     # electrons and holes
@@ -56,15 +56,15 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     lead0 = kwant.Builder(sym_left)
     lead0[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
     # hoppings in x and y-direction
-    lead0[lead0.possible_hoppings((1, 0), lat_e, lat_e)] = -t
-    lead0[lead0.possible_hoppings((0, 1), lat_e, lat_e)] = -t
+    lead0[kwant.builder.HoppingKind((1, 0), lat_e, lat_e)] = -t
+    lead0[kwant.builder.HoppingKind((0, 1), lat_e, lat_e)] = -t
 
     # left hole lead
     lead1 = kwant.Builder(sym_left)
     lead1[(lat_h(0, j) for j in xrange(W))] = mu - 4 * t
     # hoppings in x and y-direction
-    lead1[lead1.possible_hoppings((1, 0), lat_h, lat_h)] = t
-    lead1[lead1.possible_hoppings((0, 1), lat_h, lat_h)] = t
+    lead1[kwant.builder.HoppingKind((1, 0), lat_h, lat_h)] = t
+    lead1[kwant.builder.HoppingKind((0, 1), lat_h, lat_h)] = t
 #HIDDEN_END_ttth
 
     # Then the lead to the right
@@ -77,10 +77,10 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     lead2[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
     lead2[(lat_h(0, j) for j in xrange(W))] = mu - 4 * t
     # hoppings in x and y-direction
-    lead2[lead2.possible_hoppings((1, 0), lat_e, lat_e)] = -t
-    lead2[lead2.possible_hoppings((0, 1), lat_e, lat_e)] = -t
-    lead2[lead2.possible_hoppings((1, 0), lat_h, lat_h)] = t
-    lead2[lead2.possible_hoppings((0, 1), lat_h, lat_h)] = t
+    lead2[kwant.builder.HoppingKind((1, 0), lat_e, lat_e)] = -t
+    lead2[kwant.builder.HoppingKind((0, 1), lat_e, lat_e)] = -t
+    lead2[kwant.builder.HoppingKind((1, 0), lat_h, lat_h)] = t
+    lead2[kwant.builder.HoppingKind((0, 1), lat_h, lat_h)] = t
     lead2[((lat_e(0, j), lat_h(0, j)) for j in xrange(W))] = Delta
 #HIDDEN_END_mhiw