diff --git a/doc/Makefile b/doc/Makefile
index c3529688942efc29effcf7c574d861e211ff0c5c..28f86b153fab74b06e65fab437a639e557a48974 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -136,5 +136,5 @@ $(foreach name,$(SCRIPTS),$(eval $(call makedep,$(name))))
 
 # Generation of images
 .%_flag: %.py
-	cd $(dir $<) && python $(notdir $<)
+	cd $(dir $<) && python3 $(notdir $<)
 	@touch $@
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 98d2c401804651380e27101374ef1fe6a917cc47..8514d0db295d6641493ce69f077c7c25e9c03cf8 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -42,8 +42,8 @@ source_suffix = '.rst'
 master_doc = 'index'
 
 # General information about the project.
-project = u'Kwant'
-copyright = u'2011-2015, C. W. Groth (CEA), M. Wimmer, A. R. Akhmerov, X. Waintal (CEA), and others'
+project = 'Kwant'
+copyright = '2011-2015, C. W. Groth (CEA), M. Wimmer, A. R. Akhmerov, X. Waintal (CEA), and others'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -208,7 +208,7 @@ r"""\makeatletter
 # We use "et al." as it is shorter and there's not much space left horizontally.
 latex_documents = [
   ('index', 'kwant.tex', 'Kwant {0} documentation'.format(release),
-   u'C. W. Groth, M. Wimmer, A. R. Akhmerov, X. Waintal, et al.',
+   'C. W. Groth, M. Wimmer, A. R. Akhmerov, X. Waintal, et al.',
    'manual'),
 ]
 
@@ -249,9 +249,9 @@ class BoundMethodDocumenter(autodoc.FunctionDocumenter):
         # Return True iff `member` is a bound method.  Taken from
         # <http://stackoverflow.com/a/1260881>.
         return (isinstance(member, types.MethodType) and
-                member.im_self is not None and
-                not issubclass(member.im_class, type) and
-                member.im_class is not types.ClassType)
+                member.__self__ is not None and
+                not issubclass(member.__self__.__class__, type) and
+                member.__self__.__class__ is not type)
 
     def format_args(self):
         args = super(BoundMethodDocumenter, self).format_args()
diff --git a/doc/source/images/graphene.py.diff b/doc/source/images/graphene.py.diff
index 351a717cb0be71fda9a19350434553135950ed9d..d4444aff257eda32e6a21bef257f90dace0790aa 100644
--- a/doc/source/images/graphene.py.diff
+++ b/doc/source/images/graphene.py.diff
@@ -1,14 +1,13 @@
 --- original
 +++ modified
-@@ -11,6 +11,7 @@
+@@ -11,5 +11,6 @@
  #    lattice, namely graphene
  
- from __future__ import division  # so that 1/2 == 0.5, and not 0
 +import _defs
  from math import pi, sqrt, tanh
  
  import kwant
-@@ -97,22 +98,40 @@
+@@ -96,22 +97,40 @@
          smatrix = kwant.smatrix(sys, energy)
          data.append(smatrix.transmission(0, 1))
  
@@ -57,7 +56,7 @@
  
  
  def main():
-@@ -124,8 +143,11 @@
+@@ -123,8 +142,11 @@
      def family_colors(site):
          return 0 if site.family == a else 1
  
@@ -71,7 +70,7 @@
  
      # Compute some eigenvalues.
      compute_evs(sys.finalized())
-@@ -134,9 +156,11 @@
+@@ -133,9 +155,11 @@
      for lead in leads:
          sys.attach_lead(lead)
  
diff --git a/doc/source/tutorial/ab_ring.py b/doc/source/tutorial/ab_ring.py
index 5df7ecda6ef7f2331fc46b0756d7e5f5104cb126..33de87fbee803011e1222787d67577d95840e740 100644
--- a/doc/source/tutorial/ab_ring.py
+++ b/doc/source/tutorial/ab_ring.py
@@ -118,7 +118,7 @@ def main():
 
     # We should see a conductance that is periodic with the flux quantum
     plot_conductance(sys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
-                                                for i in xrange(100)])
+                                                for i in range(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/tutorial/band_structure.py b/doc/source/tutorial/band_structure.py
index ddb9976482783bee6d5d8781d4f2a5d95ab7190c..5b223d3e96f0f474e7b684f2eae779c160698388 100644
--- a/doc/source/tutorial/band_structure.py
+++ b/doc/source/tutorial/band_structure.py
@@ -24,7 +24,7 @@ def make_lead(a=1, t=1.0, W=10):
 
     # build up one unit cell of the lead, and add the hoppings
     # to the next unit cell
-    for j in xrange(W):
+    for j in range(W):
         lead[lat(0, j)] = 4 * t
 
         if j > 0:
diff --git a/doc/source/tutorial/closed_system.py b/doc/source/tutorial/closed_system.py
index 2f54bf046023d678ec9d47873c0c238bccccb077..49d2dc1d460252bde38bf2499190958947a2d01e 100644
--- a/doc/source/tutorial/closed_system.py
+++ b/doc/source/tutorial/closed_system.py
@@ -106,7 +106,7 @@ def main():
     try:
         # We should observe energy levels that flow towards Landau
         # level energies with increasing magnetic field.
-        plot_spectrum(sys, [iB * 0.002 for iB in xrange(100)])
+        plot_spectrum(sys, [iB * 0.002 for iB in range(100)])
 
         # Plot an eigenmode of a circular dot. Here we create a larger system for
         # better spatial resolution.
@@ -114,8 +114,8 @@ def main():
         plot_wave_function(sys)
     except ValueError as e:
         if e.message == "Input matrix is not real-valued.":
-            print "The calculation of eigenvalues failed because of a bug in SciPy 0.9."
-            print "Please upgrade to a newer version of SciPy."
+            print("The calculation of eigenvalues failed because of a bug in SciPy 0.9.")
+            print("Please upgrade to a newer version of SciPy.")
         else:
             raise
 
diff --git a/doc/source/tutorial/graphene.py b/doc/source/tutorial/graphene.py
index cd0ec2e6939cbd2a54bf7ff87a6e444d0ac34126..31f2c578501323a16983b698e9776a7545c30fff 100644
--- a/doc/source/tutorial/graphene.py
+++ b/doc/source/tutorial/graphene.py
@@ -10,7 +10,6 @@
 #  - Application of all the aspects of tutorials 1-3 to a more complicated
 #    lattice, namely graphene
 
-from __future__ import division  # so that 1/2 == 0.5, and not 0
 from math import pi, sqrt, tanh
 
 import kwant
@@ -102,7 +101,7 @@ def compute_evs(sys):
     sparse_mat = sys.hamiltonian_submatrix(sparse=True)
 
     evs = sla.eigs(sparse_mat, 2)[0]
-    print evs.real
+    print(evs.real)
 #HIDDEN_END_zydk
 
 
@@ -166,11 +165,11 @@ def main():
     sys = sys.finalized()
 
     # Compute the band structure of lead 0.
-    momenta = [-pi + 0.02 * pi * i for i in xrange(101)]
+    momenta = [-pi + 0.02 * pi * i for i in range(101)]
     plot_bandstructure(sys.leads[0], momenta)
 
     # Plot conductance.
-    energies = [-2 * pot + 4. / 50. * pot * i for i in xrange(51)]
+    energies = [-2 * pot + 4. / 50. * pot * i for i in range(51)]
     plot_conductance(sys, energies)
 
 
diff --git a/doc/source/tutorial/quantum_well.py b/doc/source/tutorial/quantum_well.py
index 76ab616624bba5a5cf7cf30254ae9eef6745a102..98309437c614b47418846c04dc7e8167d3a2168b 100644
--- a/doc/source/tutorial/quantum_well.py
+++ b/doc/source/tutorial/quantum_well.py
@@ -43,7 +43,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
 
     #### Define and attach the leads. ####
     lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
-    lead[(lat(0, j) for j in xrange(W))] = 4 * t
+    lead[(lat(0, j) for j in range(W))] = 4 * t
     lead[lat.neighbors()] = -t
     sys.attach_lead(lead)
     sys.attach_lead(lead.reversed())
@@ -79,7 +79,7 @@ def main():
 
     # We should see conductance steps.
     plot_conductance(sys, energy=0.2,
-                     welldepths=[0.01 * i for i in xrange(100)])
+                     welldepths=[0.01 * i for i in range(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/tutorial/quantum_wire.py b/doc/source/tutorial/quantum_wire.py
index d4d89956c27d9bb322b99660135e2bda1273986c..b32072123c2f65b492c6f5588103ad26cec3f55b 100644
--- a/doc/source/tutorial/quantum_wire.py
+++ b/doc/source/tutorial/quantum_wire.py
@@ -35,8 +35,8 @@ L = 30
 
 # Define the scattering region
 
-for i in xrange(L):
-    for j in xrange(W):
+for i in range(L):
+    for j in range(W):
         # On-site Hamiltonian
         sys[lat(i, j)] = 4 * t
 
@@ -59,7 +59,7 @@ left_lead = kwant.Builder(sym_left_lead)
 #HIDDEN_END_xcmc
 
 #HIDDEN_BEGIN_ndez
-for j in xrange(W):
+for j in range(W):
     left_lead[lat(0, j)] = 4 * t
     if j > 0:
         left_lead[lat(0, j), lat(0, j - 1)] = -t
@@ -75,7 +75,7 @@ sys.attach_lead(left_lead)
 sym_right_lead = kwant.TranslationalSymmetry((a, 0))
 right_lead = kwant.Builder(sym_right_lead)
 
-for j in xrange(W):
+for j in range(W):
     right_lead[lat(0, j)] = 4 * t
     if j > 0:
         right_lead[lat(0, j), lat(0, j - 1)] = -t
@@ -98,7 +98,7 @@ sys = sys.finalized()
 #HIDDEN_BEGIN_buzn
 energies = []
 data = []
-for ie in xrange(100):
+for ie in range(100):
     energy = ie * 0.01
 
     # compute the scattering matrix at a given energy
diff --git a/doc/source/tutorial/quantum_wire_revisited.py b/doc/source/tutorial/quantum_wire_revisited.py
index c42912ce10636e3d0f4a08d7015ef60b7df1487c..e21341576d0ed53807c46020ed763eec4cb547cc 100644
--- a/doc/source/tutorial/quantum_wire_revisited.py
+++ b/doc/source/tutorial/quantum_wire_revisited.py
@@ -39,7 +39,7 @@ def make_system(a=1, t=1.0, W=10, L=30):
     # Construct the left lead.
 #HIDDEN_BEGIN_iepx
     lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
-    lead[(lat(0, j) for j in xrange(W))] = 4 * t
+    lead[(lat(0, j) for j in range(W))] = 4 * t
     lead[lat.neighbors()] = -t
 #HIDDEN_END_iepx
 
@@ -79,7 +79,7 @@ def main():
     sys = sys.finalized()
 
     # We should see conductance steps.
-    plot_conductance(sys, energies=[0.01 * i for i in xrange(100)])
+    plot_conductance(sys, energies=[0.01 * i for i in range(100)])
 #HIDDEN_END_cjel
 
 
diff --git a/doc/source/tutorial/spin_orbit.py b/doc/source/tutorial/spin_orbit.py
index 0112d3feb9874e184303327555c4c19c0898fdcb..e9d3f7f4ff73764c3c236d58ce3f60c6b3edab75 100644
--- a/doc/source/tutorial/spin_orbit.py
+++ b/doc/source/tutorial/spin_orbit.py
@@ -55,7 +55,7 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
 
 #HIDDEN_BEGIN_yliu
-    lead[(lat(0, j) for j in xrange(W))] = 4 * t * sigma_0 + e_z * sigma_z
+    lead[(lat(0, j) for j in range(W))] = 4 * t * sigma_0 + e_z * sigma_z
     # hoppings in x-direction
     lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
         -t * sigma_0 - 1j * alpha * sigma_y
@@ -95,7 +95,7 @@ def main():
     sys = sys.finalized()
 
     # We should see non-monotonic conductance steps.
-    plot_conductance(sys, energies=[0.01 * i - 0.3 for i in xrange(100)])
+    plot_conductance(sys, energies=[0.01 * i - 0.3 for i in range(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/tutorial/superconductor_band_structure.py b/doc/source/tutorial/superconductor_band_structure.py
index 9f8da5e274daf2dbd884a69099889772ac1cefd6..6a7928363e420abbee4f24c56ac14ebeebe4dcf5 100644
--- a/doc/source/tutorial/superconductor_band_structure.py
+++ b/doc/source/tutorial/superconductor_band_structure.py
@@ -35,7 +35,7 @@ def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
 
     # build up one unit cell of the lead, and add the hoppings
     # to the next unit cell
-    for j in xrange(W):
+    for j in range(W):
         lead[lat(0, j)] = (4 * t - mu) * tau_z + Delta * tau_x
 
         if j > 0:
diff --git a/doc/source/tutorial/superconductor_transport.py b/doc/source/tutorial/superconductor_transport.py
index 505479e28a1989bd94521105d0bb90a70ce88cd6..a8184c95f3b4c5de4ce27855f7bb55fc360d994c 100644
--- a/doc/source/tutorial/superconductor_transport.py
+++ b/doc/source/tutorial/superconductor_transport.py
@@ -55,12 +55,12 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
 
     # left electron lead
     lead0 = kwant.Builder(sym_left)
-    lead0[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
+    lead0[(lat_e(0, j) for j in range(W))] = 4 * t - mu
     lead0[lat_e.neighbors()] = -t
 
     # left hole lead
     lead1 = kwant.Builder(sym_left)
-    lead1[(lat_h(0, j) for j in xrange(W))] = mu - 4 * t
+    lead1[(lat_h(0, j) for j in range(W))] = mu - 4 * t
     lead1[lat_h.neighbors()] = t
 #HIDDEN_END_ttth
 
@@ -72,7 +72,7 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     lead2 = kwant.Builder(sym_right)
     lead2 += lead0
     lead2 += lead1
-    lead2[((lat_e(0, j), lat_h(0, j)) for j in xrange(W))] = Delta
+    lead2[((lat_e(0, j), lat_h(0, j)) for j in range(W))] = Delta
 #HIDDEN_END_mhiw
 
     #### Attach the leads and return the system. ####
@@ -113,7 +113,7 @@ def main():
     # Finalize the system.
     sys = sys.finalized()
 
-    plot_conductance(sys, energies=[0.002 * i for i in xrange(100)])
+    plot_conductance(sys, energies=[0.002 * i for i in range(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/examples/3d_plot.py b/examples/3d_plot.py
index 8dcef1f40ca45a341d16afb8a208a9ed523e311a..ed8bccfeaa5bd7ea03cc6236f08a834a0c835f04 100644
--- a/examples/3d_plot.py
+++ b/examples/3d_plot.py
@@ -1,6 +1,6 @@
 from math import sqrt
 import scipy.sparse.linalg as sla
-import matplotlib.pyplot
+from matplotlib import pyplot
 import kwant
 
 sys = kwant.Builder()
@@ -9,11 +9,11 @@ lat = kwant.lattice.general([(1,0,0), (0,1,0), (0,0,1)])
 t = 1.0
 R = 10
 
-sys[(lat(x,y,z) for x in xrange(-R-1, R+1)
-     for y in xrange(-R-1, R+1) for z in xrange(R+1)
+sys[(lat(x,y,z) for x in range(-R-1, R+1)
+     for y in range(-R-1, R+1) for z in range(R+1)
      if sqrt(x**2 + y**2 + z**2) < R + 0.01)] = 4 * t
-sys[(lat(x,y,z) for x in xrange(-2*R, 2*R + 1)
-     for y in xrange(-R, R+1) for z in xrange(-R, 0))] = 4 * t
+sys[(lat(x,y,z) for x in range(-2*R, 2*R + 1)
+     for y in range(-R, R+1) for z in range(-R, 0))] = 4 * t
 sys[lat.neighbors()] = -t
 sys = sys.finalized()
 kwant.plot(sys)
diff --git a/examples/advanced_construction.py b/examples/advanced_construction.py
index dcfcfa92e6291fc577cd529732310b26d7b79c00..1d4806bacbd0eb268e0f9c7524a84b91cf98c0a1 100644
--- a/examples/advanced_construction.py
+++ b/examples/advanced_construction.py
@@ -1,6 +1,5 @@
 """An example of advanced system creation."""
 
-from __future__ import division
 from math import tanh
 from cmath import exp
 import tinyarray as ta
@@ -48,7 +47,7 @@ def make_system(R):
 
 def main():
     sys = make_system(100)
-    print kwant.smatrix(sys, 1.1, [0.1]).transmission(0, 1)
+    print(kwant.smatrix(sys, 1.1, [0.1]).transmission(0, 1))
 
 
 if __name__ == '__main__':
diff --git a/examples/logo.py b/examples/logo.py
index c3dcb1f6e59cef8b4e7427895891767f2fb56043..941c613c8c803205d139a91bb24fd6826f960e0b 100644
--- a/examples/logo.py
+++ b/examples/logo.py
@@ -1,12 +1,10 @@
 """The script generating Kwant logo. In addition to Kwant it also needs Python
-image library (PIL)."""
+image library Pillow."""
 
-import Image
-import ImageFont
-import ImageDraw
+from PIL import Image, ImageFont, ImageDraw
 import matplotlib
 import numpy as np
-import scipy
+import scipy.misc
 import kwant
 
 def main():
@@ -86,7 +84,7 @@ def main():
     # result is not too empty or not too dark.
     out = np.zeros(textpos.shape)
     for i, rho in enumerate(ldos**.2):
-        x1, y1 = sys.site(i).tag
+        x1, y1 = sys.sites[i].tag
         out[x1, y1] = rho
     out = normalize_data(out)
 
diff --git a/examples/square.py b/examples/square.py
index c7df5d54200b247cb7584db9a7e0a9699ebf2f58..580df113f33af5743bc32140a460aba003aa6ca8 100644
--- a/examples/square.py
+++ b/examples/square.py
@@ -2,7 +2,6 @@
 kwant.Builder.
 """
 
-from __future__ import division
 import numpy as np
 from matplotlib import pyplot
 import kwant
@@ -10,6 +9,7 @@ from kwant.physics.leads import square_selfenergy
 
 __all__ = ['System']
 
+
 class Lead(object):
     def __init__(self, width, t, potential):
         self.width = width
@@ -20,6 +20,7 @@ class Lead(object):
         return square_selfenergy(self.width, self.t,
                                  self.potential + fermi_energy)
 
+
 class System(kwant.system.FiniteSystem):
     def __init__(self, shape, hopping,
                  potential=0, lead_potentials=(0, 0),
@@ -56,7 +57,7 @@ class System(kwant.system.FiniteSystem):
             edges[:shape[across], 1] += increment[along]
             edges[shape[across]:, (0, 1)] = edges[:shape[across], (1, 0)]
             g.add_edges(edges)
-            for i in xrange(shape[along] - 2):
+            for i in range(shape[along] - 2):
                 edges += increment[along]
                 g.add_edges(edges)
         self.graph = g.compressed()
@@ -66,7 +67,7 @@ class System(kwant.system.FiniteSystem):
             # We have to use list here, as numpy.array does not understand
             # generators.
             interface = list(self.nodeid_from_pos((x, y))
-                             for y in xrange(shape[1]))
+                             for y in range(shape[1]))
             self.lead_interfaces.append(np.array(interface))
 
         self.leads = [Lead(shape[1], hopping, lead_potentials[i])
@@ -85,7 +86,7 @@ class System(kwant.system.FiniteSystem):
         return result
 
     def nodeid_from_pos(self, pos):
-        for i in xrange(2):
+        for i in range(2):
             assert int(pos[i]) == pos[i]
             assert pos[i] >= 0 and pos[i] < self.shape[i]
         return pos[0] + pos[1] * self.shape[0]
@@ -98,8 +99,8 @@ class System(kwant.system.FiniteSystem):
 
 def main():
     sys = System((10, 5), 1)
-    energies = [0.04 * i for i in xrange(100)]
-    data = [kwant.smatrix(sys, energy).transmission(1, 0)
+    energies = [0.04 * i for i in range(100)]
+    data = [kwant.greens_function(sys, energy).transmission(1, 0)
             for energy in energies]
 
     pyplot.plot(energies, data)
diff --git a/examples/tests/test_square.py b/examples/tests/test_square.py
index 4d605368caa4f43660fbd27ea94587aecf09ef3f..f8654c98a50288f49d55a6b079029a959090d5aa 100644
--- a/examples/tests/test_square.py
+++ b/examples/tests/test_square.py
@@ -4,18 +4,18 @@ from numpy.testing import assert_almost_equal
 
 def test_nodeid_to_from_pos():
     s = square.System((3, 4), 1)
-    assert_raises(StandardError, s.nodeid_from_pos, (0, -2))
-    assert_raises(StandardError, s.nodeid_from_pos, (-1, 3))
-    assert_raises(StandardError, s.nodeid_from_pos, (3, 1))
-    assert_raises(StandardError, s.pos_from_nodeid, -1)
-    assert_raises(StandardError, s.pos_from_nodeid, 12)
+    assert_raises(Exception, s.nodeid_from_pos, (0, -2))
+    assert_raises(Exception, s.nodeid_from_pos, (-1, 3))
+    assert_raises(Exception, s.nodeid_from_pos, (3, 1))
+    assert_raises(Exception, s.pos_from_nodeid, -1)
+    assert_raises(Exception, s.pos_from_nodeid, 12)
     assert_equal(s.nodeid_from_pos((0, 0)), 0)
     assert_equal(s.nodeid_from_pos(s.pos_from_nodeid(7)), 7)
     assert_equal(s.pos_from_nodeid(s.nodeid_from_pos((2, 3))), (2, 3))
 
 def test_hamiltonian():
     sys = square.System((4, 5), 1)
-    for i in xrange(sys.graph.num_nodes):
+    for i in range(sys.graph.num_nodes):
         shape = sys.hamiltonian(i, i).shape
         assert_equal(len(shape), 2)
         assert_equal(shape[0], 1)
@@ -28,7 +28,7 @@ def test_hamiltonian():
 
 def test_selfenergy():
     sys = square.System((2, 4), 1)
-    for lead in xrange(len(sys.lead_interfaces)):
+    for lead in range(len(sys.lead_interfaces)):
         se = sys.leads[lead].selfenergy(0)
         assert_equal(len(se.shape), 2)
         assert_equal(se.shape[0], se.shape[1])
diff --git a/kwant/__init__.py b/kwant/__init__.py
index c45de3fd9fac82e270e530bb39c82630e08399a4..74e80f9fb3bb1491cbee5b721c40c0a952fe5165 100644
--- a/kwant/__init__.py
+++ b/kwant/__init__.py
@@ -29,7 +29,7 @@ __all__.append('KwantDeprecationWarning')
 from ._common import version as __version__
 
 for module in ['system', 'builder', 'lattice', 'solvers', 'digest', 'rmt']:
-    exec 'from . import {0}'.format(module)
+    exec('from . import {0}'.format(module))
     __all__.append(module)
 
 # Make selected functionality available directly in the root namespace.
@@ -38,7 +38,7 @@ available = [('builder', ['Builder', 'HoppingKind']),
              ('solvers.default',
               ['smatrix', 'greens_function', 'ldos', 'wave_function'])]
 for module, names in available:
-    exec 'from .{0} import {1}'.format(module, ', '.join(names))
+    exec('from .{0} import {1}'.format(module, ', '.join(names)))
     __all__.extend(names)
 
 # Importing plotter might not work, but this does not have to be a problem --
diff --git a/kwant/_common.py b/kwant/_common.py
index e050f2edef9923a1aec16cbb46b1a2bb1fe767ff..b67d1326a53d3af9758fb752e27a12924661c9b0 100644
--- a/kwant/_common.py
+++ b/kwant/_common.py
@@ -11,7 +11,9 @@ import os
 
 __all__ = ['version', 'KwantDeprecationWarning']
 
-distr_root = os.path.dirname(os.path.dirname(os.path.realpath(__file__)))
+package_root = os.path.dirname(os.path.realpath(__file__))
+distr_root = os.path.dirname(package_root)
+version_file = '_kwant_version.py'
 
 def get_version_from_git():
     try:
@@ -22,8 +24,7 @@ def get_version_from_git():
         return
     if p.wait() != 0:
         return
-    # TODO: use os.path.samefile once we depend on Python >= 3.3.
-    if os.path.normpath(p.communicate()[0].rstrip('\n')) != distr_root:
+    if not os.path.samefile(p.communicate()[0].decode().rstrip('\n'), distr_root):
         # The top-level directory of the current Git repository is not the same
         # as the root directory of the Kwant distribution: do not extract the
         # version from Git.
@@ -42,7 +43,7 @@ def get_version_from_git():
             break
     else:
         return
-    description = p.communicate()[0].strip('v').rstrip('\n')
+    description = p.communicate()[0].decode().strip('v').rstrip('\n')
 
     release, dev, git = description.rsplit('-', 2)
     version = [release]
@@ -67,7 +68,11 @@ def get_version_from_git():
 
 
 
-from _kwant_version import version
+# populate the version_info dictionary with values stored in the version file
+version_info = {}
+with open(os.path.join(package_root, version_file), 'r') as f:
+    exec(f.read(), {}, version_info)
+version = version_info['version']
 version_is_from_git = (version == "__use_git__")
 if version_is_from_git:
     version = get_version_from_git()
diff --git a/kwant/_system.pyx b/kwant/_system.pyx
index dc3e665dd3ef9505cacb852226943ca946e82e62..beebd9465866dd61532b1f8d92ef411a9ee07555 100644
--- a/kwant/_system.pyx
+++ b/kwant/_system.pyx
@@ -11,6 +11,7 @@ import tinyarray as ta
 import numpy as np
 from scipy import sparse as sp
 from itertools import chain
+import types
 
 from .graph.core cimport CGraph, gintArraySlice
 from .graph.defs cimport gint
@@ -84,8 +85,15 @@ def make_sparse(ham, args, CGraph gr, diag,
                         rows_cols[1, k] = j + from_off[n_fs]
                         k += 1
 
-    return sp.coo_matrix((data[:k], rows_cols[:, :k]),
-                         shape=(to_off[-1], from_off[-1]))
+    # hack around a bug in Scipy + Python 3 + memoryviews
+    # see https://github.com/scipy/scipy/issues/5123 for details
+    np_data = np.asarray(data)
+    np_rows_cols = np.asarray(rows_cols)
+    np_to_off = np.asarray(to_off)
+    np_from_off = np.asarray(from_off)
+
+    return sp.coo_matrix((np_data[:k], np_rows_cols[:, :k]),
+                         shape=(np_to_off[-1], np_from_off[-1]))
 
 
 @cython.boundscheck(False)
@@ -147,8 +155,15 @@ def make_sparse_full(ham, args, CGraph gr, diag,
                         rows_cols[0, k + 1] = rows_cols[1, k] = j + from_off[fs]
                         k += 2
 
-    return sp.coo_matrix((data[:k], rows_cols[:, :k]),
-                         shape=(to_off[-1], from_off[-1]))
+    # hack around a bug in Scipy + Python 3 + memoryviews
+    # see https://github.com/scipy/scipy/issues/5123 for details
+    np_data = np.asarray(data)
+    np_rows_cols = np.asarray(rows_cols)
+    np_to_off = np.asarray(to_off)
+    np_from_off = np.asarray(from_off)
+
+    return sp.coo_matrix((np_data[:k], np_rows_cols[:, :k]),
+                         shape=(np_to_off[-1], np_from_off[-1]))
 
 
 @cython.boundscheck(False)
@@ -327,3 +342,13 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
         mat = func(ham, args, self.graph, diag, from_sites,
                    n_by_to_site, to_norb, to_off, from_norb, from_off)
     return (mat, to_norb, from_norb) if return_norb else mat
+
+# workaround for Cython functions not having __get__ and
+# Python 3 getting rid of unbound methods
+cdef class HamiltonianSubmatrix:
+
+    def __get__(self, obj, objtype):
+        if obj is None:
+            return hamiltonian_submatrix
+        else:
+            return types.MethodType(hamiltonian_submatrix, obj)
diff --git a/kwant/builder.py b/kwant/builder.py
index c594f7202322e492abfddd78f22f6d85ce09e341..cd83747df7112a8fb35169dd644db138050f5d89 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -6,15 +6,13 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
-
 __all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
            'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead']
 
 import abc
 import warnings
 import operator
-from itertools import izip, islice, chain
+from itertools import islice, chain
 import tinyarray as ta
 import numpy as np
 from . import system, graph, KwantDeprecationWarning
@@ -84,7 +82,7 @@ class Site(tuple):
         return self.family.pos(self.tag)
 
 
-class SiteFamily(object):
+class SiteFamily(object, metaclass=abc.ABCMeta):
     """Abstract base class for site families.
 
     Site families are the 'type' of `Site` objects.  Within a family, individual
@@ -107,7 +105,6 @@ class SiteFamily(object):
     site belonging to this family with a given tag.
 
     """
-    __metaclass__ = abc.ABCMeta
 
     def __init__(self, canonical_repr, name):
         self.canonical_repr = canonical_repr
@@ -224,7 +221,7 @@ def validate_hopping(hopping):
 
 ################ Symmetries
 
-class Symmetry(object):
+class Symmetry(object, metaclass=abc.ABCMeta):
     """Abstract base class for spatial symmetries.
 
     Many physical systems possess a discrete spatial symmetry, which results in
@@ -249,7 +246,6 @@ class Symmetry(object):
     typical example of this is when the vector defining a translational
     symmetry is not a lattice vector.
     """
-    __metaclass__ = abc.ABCMeta
 
     @abc.abstractproperty
     def num_directions(self):
@@ -430,7 +426,7 @@ class HermConjOfFunc(object):
 
 ################ Leads
 
-class Lead(object):
+class Lead(object, metaclass=abc.ABCMeta):
     """Abstract base class for leads that can be attached to a `Builder`.
 
     To attach a lead to a builder, append it to the builder's `~Builder.leads`
@@ -442,7 +438,6 @@ class Lead(object):
     interface : sequence of sites
 
     """
-    __metaclass__ = abc.ABCMeta
 
     @abc.abstractmethod
     def finalized(self):
@@ -567,7 +562,7 @@ def edges(seq):
     # izip, when given the same iterator twice, turns a sequence into a
     # sequence of pairs.
     seq_iter = iter(seq)
-    result = izip(seq_iter, seq_iter)
+    result = zip(seq_iter, seq_iter)
     next(result)                # Skip the special loop edge.
     return result
 
@@ -766,7 +761,7 @@ class Builder(object):
         result.H = self.H
         return result
 
-    def __nonzero__(self):
+    def __bool__(self):
         return bool(self.H)
 
     def expand(self, key):
@@ -974,13 +969,13 @@ class Builder(object):
         initially (but always the equivalent ones).
         """
         try:
-            return self.H.viewkeys()
+            return self.H.keys()
         except AttributeError:
             return frozenset(self.H)
 
     def site_value_pairs(self):
         """Return an iterator over all (site, value) pairs."""
-        for site, hvhv in self.H.iteritems():
+        for site, hvhv in self.H.items():
             yield site, hvhv[1]
 
     def hoppings(self):
@@ -990,7 +985,7 @@ class Builder(object):
         `Builder` symmetry, and are not necessarily the ones that were set
         initially (but always the equivalent ones).
         """
-        for tail, hvhv in self.H.iteritems():
+        for tail, hvhv in self.H.items():
             for head, value in edges(hvhv):
                 if value is Other:
                     continue
@@ -998,7 +993,7 @@ class Builder(object):
 
     def hopping_value_pairs(self):
         """Return an iterator over all (hopping, value) pairs."""
-        for tail, hvhv in self.H.iteritems():
+        for tail, hvhv in self.H.items():
             for head, value in edges(hvhv):
                 if value is Other:
                     continue
@@ -1201,7 +1196,7 @@ class Builder(object):
         #### Make graph.
         g = graph.Graph()
         g.num_nodes = len(sites)  # Some sites could not appear in any edge.
-        for tail, hvhv in self.H.iteritems():
+        for tail, hvhv in self.H.items():
             for head in islice(hvhv, 2, None, 2):
                 if tail == head:
                     continue
@@ -1225,7 +1220,7 @@ class Builder(object):
                     msg = 'When finalizing lead {0}:'.format(lead_nr)
                     warnings.warn(w.__class__(' '.join((msg,) + w.args)),
                                   stacklevel=3)
-            except ValueError, e:
+            except ValueError as e:
                 # Re-raise the exception with an additional message.
                 msg = 'Problem finalizing lead {0}:'.format(lead_nr)
                 e.args = (' '.join((msg,) + e.args),)
diff --git a/kwant/digest.py b/kwant/digest.py
index 8daac5cc6483aad90f12d9788bb75a7235670ec2..0b502d5f2992477ee12d721111a104b0fa481536 100644
--- a/kwant/digest.py
+++ b/kwant/digest.py
@@ -19,8 +19,6 @@ good enough to pass the "dieharder" battery of tests: see the function `test` of
 this module.
 """
 
-from __future__ import division
-
 from math import pi, log, sqrt, cos
 from hashlib import md5
 from struct import unpack
@@ -35,30 +33,22 @@ BPF_MASK = 2**53 - 1
 RECIP_BPF = 2**-BPF
 
 
-# TODO: Remove the following workaround for Python 2.6 once we do not support it
-# anymore.
-
-if sys.version_info < (2, 7):
-    def uniform2(input, salt=''):
-        """Return two independent [0,1)-distributed numbers."""
-        try:
-            input = bytes(buffer(input)) + salt
-        except TypeError:
-            # Tinyarray does not provide the old buffer protocol, so buffer does
-            # not work.  However, bytearray does work!
-            input = bytearray(input) + salt
-        a, b = unpack('qq', md5(input).digest())
-        a &= BPF_MASK
-        b &= BPF_MASK
-        return a * RECIP_BPF, b * RECIP_BPF
-else:
-    def uniform2(input, salt=''):
-        """Return two independent [0,1)-distributed numbers."""
-        input = memoryview(input).tobytes() + salt
-        a, b = unpack('qq', md5(input).digest())
-        a &= BPF_MASK
-        b &= BPF_MASK
-        return a * RECIP_BPF, b * RECIP_BPF
+def str_to_bytes(s):
+    """Return bytes if the input is a string, else return the object as-is."""
+    try:
+        return s.encode('utf8')
+    except AttributeError:
+        return s
+
+def uniform2(input, salt=''):
+    """Return two independent [0,1)-distributed numbers."""
+    input = str_to_bytes(input)
+    salt = str_to_bytes(salt)
+    input = memoryview(input).tobytes() + salt
+    a, b = unpack('qq', md5(input).digest())
+    a &= BPF_MASK
+    b &= BPF_MASK
+    return a * RECIP_BPF, b * RECIP_BPF
 
 
 def uniform(input, salt=''):
@@ -96,8 +86,8 @@ def test(n=20000):
 
     f = tempfile.NamedTemporaryFile(delete=False)
     try:
-        for x in xrange(n):
-            for y in xrange(n):
+        for x in range(n):
+            for y in range(n):
                 a = array((x, y))
                 i = int(2**32 * uniform(a))
                 f.write(pack('I', i))
diff --git a/kwant/graph/__init__.py b/kwant/graph/__init__.py
index dc9463446aae5d3e7152d171442c0ea368838036..2905777027cb8b0672367c507000794e9dde1f1b 100644
--- a/kwant/graph/__init__.py
+++ b/kwant/graph/__init__.py
@@ -11,6 +11,6 @@
 # Merge the public interface of all submodules.
 __all__ = []
 for module in ['core', 'defs', 'slicer', 'utils']:
-    exec 'from . import {0}'.format(module)
-    exec 'from .{0} import *'.format(module)
-    exec '__all__.extend({0}.__all__)'.format(module)
+    exec('from . import {0}'.format(module))
+    exec('from .{0} import *'.format(module))
+    exec('__all__.extend({0}.__all__)'.format(module))
diff --git a/kwant/graph/tests/test_core.py b/kwant/graph/tests/test_core.py
index 967ff6d8ad95043e3515a2d8e6e340f794fa947c..14cdad0ba633b1989a1341450ee82472ba6df86e 100644
--- a/kwant/graph/tests/test_core.py
+++ b/kwant/graph/tests/test_core.py
@@ -6,8 +6,8 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from StringIO import StringIO
-from itertools import izip_longest
+from io import StringIO
+from itertools import zip_longest
 import numpy as np
 from nose.tools import assert_equal, assert_raises
 from kwant.graph.core import (Graph, NodeDoesNotExistError,
@@ -41,12 +41,12 @@ def test_num_nodes():
 def test_large():
     num_edges = 1000
     graph = Graph()
-    for i in xrange(num_edges):
+    for i in range(num_edges):
         graph.add_edge(i, i + 1)
     g = graph.compressed()
     g2 = graph.compressed(twoway=True)
     assert_equal(num_edges, g.num_nodes - 1)
-    for i in xrange(num_edges):
+    for i in range(num_edges):
         assert_equal(tuple(g.out_neighbors(i)), (i + 1,))
         assert_equal(tuple(g2.in_neighbors(i + 1)), (i,))
 
@@ -71,7 +71,7 @@ def test_small():
     check_dot(dot_expect, g)
     g = g.compressed(twoway=True)
 
-    for edge_should, edge_is in izip_longest(edges, g):
+    for edge_should, edge_is in zip_longest(edges, g):
         assert_equal(edge_should, edge_is)
 
     edge_ids = []
diff --git a/kwant/graph/tests/test_dissection.py b/kwant/graph/tests/test_dissection.py
index cb8fc8f7f1e0209a3eef171b54974fcaeab67516..fb21fe786737a2f784d6d1de09df528ce4dfc1d8 100644
--- a/kwant/graph/tests/test_dissection.py
+++ b/kwant/graph/tests/test_dissection.py
@@ -19,13 +19,13 @@ def _DISABLED_test_edge_dissection():
     size = 5
     graph = Graph()
 
-    for i in xrange(size - 1):
+    for i in range(size - 1):
         offset = i * size
-        for j in xrange(size - 1):
+        for j in range(size - 1):
             graph.add_edge(offset + j, offset + j + 1)
             graph.add_edge(offset + j + 1, offset + j)
         if i > 0:
-            for j in xrange(size):
+            for j in range(size):
                 graph.add_edge(offset + j, offset + j - size)
                 graph.add_edge(offset + j - size, offset + j)
     g = graph.compressed()
diff --git a/kwant/graph/tests/test_scotch.py b/kwant/graph/tests/test_scotch.py
index 222ff05f5f9e0b2b1e6a6eaa01c9b217801cc994..2ded31ae2f863eceee4a371b0e2db6a35eb70788 100644
--- a/kwant/graph/tests/test_scotch.py
+++ b/kwant/graph/tests/test_scotch.py
@@ -17,32 +17,32 @@ def _DISABLED_test_bisect():
     size = 5
     graph = Graph()
 
-    for i in xrange(size-1):
+    for i in range(size - 1):
         offset = i * size
-        for j in xrange(size-1):
+        for j in range(size - 1):
             graph.add_edge(offset + j, offset + j + 1)
             graph.add_edge(offset + j + 1, offset + j)
         if i > 0:
-            for j in xrange(size):
+            for j in range(size):
                 graph.add_edge(offset + j, offset + j - size)
                 graph.add_edge(offset + j - size, offset + j)
     g = graph.compressed()
 
     parts = bisect(g)
-    for i in xrange(g.num_nodes):
+    for i in range(g.num_nodes):
         assert_true(parts[i] == 0 or parts[i] == 1)
 
 def _DISABLED_test_reset():
     size = 5
     graph = Graph()
 
-    for i in xrange(size-1):
+    for i in range(size - 1):
         offset = i * size
-        for j in xrange(size-1):
+        for j in range(size - 1):
             graph.add_edge(offset + j, offset + j + 1)
             graph.add_edge(offset + j + 1, offset + j)
         if i > 0:
-            for j in xrange(size):
+            for j in range(size):
                 graph.add_edge(offset + j, offset + j - size)
                 graph.add_edge(offset + j - size, offset + j)
     g = graph.compressed()
diff --git a/kwant/graph/tests/test_slicer.py b/kwant/graph/tests/test_slicer.py
index 29e549cbdef4c04e568a6eecbea44e246a046d2d..2191310dbc759c0f64dd8f2ec2224ef5a836a9ad 100644
--- a/kwant/graph/tests/test_slicer.py
+++ b/kwant/graph/tests/test_slicer.py
@@ -11,13 +11,13 @@ from kwant.graph import slicer
 
 def assert_sanity(graph, slices):
     # Slices must comprise all of the graph.
-    slclist = [slices[j][i] for j in xrange(len(slices))
-               for i in xrange(len(slices[j]))]
+    slclist = [slices[j][i] for j in range(len(slices))
+               for i in range(len(slices[j]))]
     slclist.sort()
-    assert slclist == [i for i in xrange(graph.num_nodes)]
+    assert slclist == [i for i in range(graph.num_nodes)]
 
     # Nodes may only have neighbors in neighboring slices.
-    for j in xrange(len(slices)):
+    for j in range(len(slices)):
         for node in slices[j]:
             for neigh in graph.out_neighbors(node):
                 if j > 0 and j < len(slices) - 1:
@@ -39,8 +39,8 @@ def test_rectangle():
         sys = kwant.Builder()
         lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
         lat = kwant.lattice.square()
-        lead[(lat(0, i) for i in xrange(w))] = 0
-        sys[(lat(j, i) for j in xrange(l) for i in xrange(w))] = 0
+        lead[(lat(0, i) for i in range(w))] = 0
+        sys[(lat(j, i) for j in range(l) for i in range(w))] = 0
         for s in [lead, sys]:
             for kind in [kwant.builder.HoppingKind((1, 0), lat),
                          kwant.builder.HoppingKind((0, 1), lat)]:
@@ -57,7 +57,7 @@ def test_rectangle():
         # we know that all slices must have the same shape.
         assert len(slices) == l
 
-        for j in xrange(l):
+        for j in range(l):
             assert len(slices[j]) == w
 
         assert_sanity(fsys.graph, slices)
diff --git a/kwant/graph/tests/test_utils.py b/kwant/graph/tests/test_utils.py
index bea187e5f38344e1407dba3ea277712806eaf67e..0906765e5fd587dfe7fa2345bf06cb44c6dc8748 100644
--- a/kwant/graph/tests/test_utils.py
+++ b/kwant/graph/tests/test_utils.py
@@ -83,7 +83,7 @@ def test_induced_subgraph():
     num_nodes = 6
 
     graph = Graph()
-    for i in xrange(num_nodes-1):
+    for i in range(num_nodes - 1):
         graph.add_edge(i, i + 1)
     graph.add_edge(1, 0)
     g = graph.compressed()
diff --git a/kwant/lattice.py b/kwant/lattice.py
index 5b55b98890f984aad57944bfc3b1b8bf3d5b580a..e93228980cc7392bcd5f9842f15bae745eb584eb 100644
--- a/kwant/lattice.py
+++ b/kwant/lattice.py
@@ -6,8 +6,6 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
-
 __all__ = ['TranslationalSymmetry', 'general', 'Polyatomic', 'Monatomic',
            'chain', 'square', 'triangular', 'honeycomb', 'kagome']
 
@@ -617,7 +615,6 @@ class TranslationalSymmetry(builder.Symmetry):
 
         det_m = int(round(np.linalg.det(m)))
         if det_m == 0:
-            print m
             raise RuntimeError('Adding site family failed.')
 
         det_x_inv_m = np.array(np.round(det_m * np.linalg.inv(m)), dtype=int)
diff --git a/kwant/linalg/__init__.py b/kwant/linalg/__init__.py
index ad60f34fdde3fa5ca390b90a2ea7131204c3c161..cd46cf5aab99934e9c73282c42b1ff3700835433 100644
--- a/kwant/linalg/__init__.py
+++ b/kwant/linalg/__init__.py
@@ -11,6 +11,6 @@ from . import lapack
 
 # Merge the public interface of the other submodules.
 for module in ['decomp_lu', 'decomp_ev', 'decomp_schur']:
-    exec 'from . import {0}'.format(module)
-    exec 'from .{0} import *'.format(module)
-    exec '__all__.extend({0}.__all__)'.format(module)
+    exec('from . import {0}'.format(module))
+    exec('from .{0} import *'.format(module))
+    exec('__all__.extend({0}.__all__)'.format(module))
diff --git a/kwant/linalg/decomp_lu.py b/kwant/linalg/decomp_lu.py
index 65bc0fd723ebeb2b2695caf3b09bbae8c007a45d..79accc254feaa53206a9ea2f4a9d977a0d904138 100644
--- a/kwant/linalg/decomp_lu.py
+++ b/kwant/linalg/decomp_lu.py
@@ -61,13 +61,13 @@ def lu_factor(a, overwrite_a=False):
         return lapack.cgetrf(a)
 
 
-def lu_solve((lu, ipiv, singular), b):
+def lu_solve(matrix_factorization, b):
     """Solve a linear system of equations, a x = b, given the LU
     factorization of a
 
     Parameters
     ----------
-    (lu, piv, singular)
+    matrix_factorization
         Factorization of the coefficient matrix a, as given by lu_factor
     b : array (vector or matrix)
         Right-hand side
@@ -77,7 +77,7 @@ def lu_solve((lu, ipiv, singular), b):
     x : array (vector or matrix)
         Solution to the system
     """
-
+    (lu, ipiv, singular) = matrix_factorization
     if singular:
         raise RuntimeWarning("In lu_solve: the flag singular indicates "
                              "a singular matrix. Result of solve step "
@@ -102,7 +102,7 @@ def lu_solve((lu, ipiv, singular), b):
         return lapack.cgetrs(lu, ipiv, b)
 
 
-def rcond_from_lu((lu, ipiv, singular), norm_a, norm="1"):
+def rcond_from_lu(matrix_factorization, norm_a, norm="1"):
     """Compute the reciprocal condition number from the LU decomposition as
     returned from lu_factor(), given additionally the norm of the matrix a in
     norm_a.
@@ -112,7 +112,7 @@ def rcond_from_lu((lu, ipiv, singular), norm_a, norm="1"):
 
     Parameters
     ----------
-    (lu, piv, singular)
+    matrix_factorization
         Factorization of the matrix a, as given by lu_factor
     norm_a : float or complex
         norm of the original matrix a (type of norm is specified in norm)
@@ -126,9 +126,10 @@ def rcond_from_lu((lu, ipiv, singular), norm_a, norm="1"):
         reciprocal condition number of a with respect to the type of matrix
         norm specified in norm
     """
-
+    (lu, ipiv, singular) = matrix_factorization
     if not norm in ("1", "I"):
         raise ValueError("norm in rcond_from_lu must be either '1' or 'I'")
+    norm = norm.encode('utf8')  # lapack expects bytes
 
     ltype, lu = lapack.prepare_for_lapack(False, lu)
 
diff --git a/kwant/linalg/fortran_helpers.py b/kwant/linalg/fortran_helpers.py
index 2f82b26db9410d919981677692df5b7cfdec7b38..39ec639e606264e8dc14d85e06284c4aa51b86e6 100644
--- a/kwant/linalg/fortran_helpers.py
+++ b/kwant/linalg/fortran_helpers.py
@@ -36,7 +36,7 @@ def prepare_for_fortran(overwrite, *args):
 
     # Make sure we have NumPy arrays
     mats = [None]*len(args)
-    for i in xrange(len(args)):
+    for i in range(len(args)):
         if args[i] is not None:
             arr = np.asanyarray(args[i])
             if not np.issubdtype(arr.dtype, np.number):
diff --git a/kwant/linalg/lapack.pyx b/kwant/linalg/lapack.pyx
index 5ad31afc254a3b087307b6d1d77b724303487019..da74fbc51b1fc457a21247129c12d07abcc31012 100644
--- a/kwant/linalg/lapack.pyx
+++ b/kwant/linalg/lapack.pyx
@@ -245,7 +245,7 @@ def zgetrs(np.ndarray[np.complex128_t, ndim=2] LU,
 # Wrappers for xGECON
 
 def sgecon(np.ndarray[np.float32_t, ndim=2] LU,
-            float normA, char *norm = "1"):
+            float normA, char *norm = b"1"):
     cdef l_int N, info
     cdef float rcond
     cdef np.ndarray[np.float32_t, ndim=1] work
@@ -266,7 +266,7 @@ def sgecon(np.ndarray[np.float32_t, ndim=2] LU,
     return rcond
 
 def dgecon(np.ndarray[np.float64_t, ndim=2] LU,
-            double normA, char *norm = "1"):
+            double normA, char *norm = b"1"):
     cdef l_int N, info
     cdef double rcond
     cdef np.ndarray[np.float64_t, ndim=1] work
@@ -287,7 +287,7 @@ def dgecon(np.ndarray[np.float64_t, ndim=2] LU,
     return rcond
 
 def cgecon(np.ndarray[np.complex64_t, ndim=2] LU,
-            float normA, char *norm = "1"):
+            float normA, char *norm = b"1"):
     cdef l_int N, info
     cdef float rcond
     cdef np.ndarray[np.complex64_t, ndim=1] work
@@ -308,7 +308,7 @@ def cgecon(np.ndarray[np.complex64_t, ndim=2] LU,
     return rcond
 
 def zgecon(np.ndarray[np.complex128_t, ndim=2] LU,
-           double normA, char *norm = "1"):
+           double normA, char *norm = b"1"):
     cdef l_int N, info
     cdef double rcond
     cdef np.ndarray[np.complex128_t, ndim=1] work
diff --git a/kwant/linalg/lll.py b/kwant/linalg/lll.py
index 17b353ca11fea6f39487fc527af1ef3636d76ddc..c9f2540fc68056234bc03056ae3f89d862483590 100644
--- a/kwant/linalg/lll.py
+++ b/kwant/linalg/lll.py
@@ -6,13 +6,12 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
-
 __all__ = ['lll', 'cvp', 'voronoi']
 
 import numpy as np
 from itertools import product
 
+
 def gs_coefficient(a, b):
     """Gram-Schmidt coefficient."""
     return np.dot(a, b) / np.linalg.norm(b)**2
diff --git a/kwant/linalg/tests/_test_utils.py b/kwant/linalg/tests/_test_utils.py
index 8f76e8150f5b04ed8ae5435f6065460b4f09bd7c..5ada76379e8b6a42dc5f7ace060213384de5ed90 100644
--- a/kwant/linalg/tests/_test_utils.py
+++ b/kwant/linalg/tests/_test_utils.py
@@ -43,16 +43,16 @@ class _Random:
         mat = np.empty((n, m), dtype = dtype)
 
         if issubclass(dtype, np.complexfloating):
-            for i in xrange(n):
-                for j in xrange(m):
+            for i in range(n):
+                for j in range(m):
                     mat[i,j] = self._randf() + 1j * self._randf()
         elif issubclass(dtype, np.floating):
-            for i in xrange(n):
-                for j in xrange(m):
+            for i in range(n):
+                for j in range(m):
                     mat[i,j] = self._randf()
         elif issubclass(dtype, np.integer):
-            for i in xrange(n):
-                for j in xrange(m):
+            for i in range(n):
+                for j in range(m):
                     mat[i,j] = self._randi()
 
         return mat
@@ -61,13 +61,13 @@ class _Random:
         vec = np.empty(n, dtype = dtype)
 
         if issubclass(dtype, np.complexfloating):
-            for i in xrange(n):
+            for i in range(n):
                 vec[i] = self._randf() + 1j * self._randf()
         elif issubclass(dtype, np.floating):
-            for i in xrange(n):
+            for i in range(n):
                 vec[i] = self._randf()
         elif issubclass(dtype, np.integer):
-            for i in xrange(n):
+            for i in range(n):
                 vec[i] = self._randi()
 
         return vec
diff --git a/kwant/linalg/tests/test_linalg.py b/kwant/linalg/tests/test_linalg.py
index f059580af1545f27123aa30af48ba7fb2740a778..7b70728ced8dae5d91de80c7d32e42aea0a2ceeb 100644
--- a/kwant/linalg/tests/test_linalg.py
+++ b/kwant/linalg/tests/test_linalg.py
@@ -12,7 +12,7 @@ from kwant.linalg import (
     convert_r2c_gen_schur, order_gen_schur, evecs_from_gen_schur)
 from nose.tools import assert_equal, assert_true
 import numpy as np
-from _test_utils import _Random, assert_array_almost_equal
+from ._test_utils import _Random, assert_array_almost_equal
 
 def test_gen_eig():
     def _test_gen_eig(dtype):
diff --git a/kwant/linalg/tests/test_lll.py b/kwant/linalg/tests/test_lll.py
index e11627aee7cdd48da595de150a4fc8d3434b7d45..cb04809046d722b061bdd9435c35205f15e00b65 100644
--- a/kwant/linalg/tests/test_lll.py
+++ b/kwant/linalg/tests/test_lll.py
@@ -6,7 +6,7 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
+
 import numpy as np
 from kwant.linalg import lll
 
diff --git a/kwant/linalg/tests/test_mumps.py b/kwant/linalg/tests/test_mumps.py
index 8991a7d2c4899e2783d94b0c07a758789709c3e3..bfb5028a75e8937514127bcf221d621bd5fc7c12 100644
--- a/kwant/linalg/tests/test_mumps.py
+++ b/kwant/linalg/tests/test_mumps.py
@@ -17,7 +17,7 @@ from kwant.builder import Builder, HoppingKind
 from numpy.testing.decorators import skipif
 import numpy as np
 import scipy.sparse as sp
-from _test_utils import _Random, assert_array_almost_equal
+from ._test_utils import _Random, assert_array_almost_equal
 
 @skipif(no_mumps)
 def test_lu_with_dense():
@@ -53,7 +53,7 @@ def test_schur_complement_with_dense():
     def _test_schur_complement_with_dense(dtype):
         rand = _Random()
         a = rand.randmat(10, 10, dtype)
-        s = schur_complement(sp.coo_matrix(a), range(3))
+        s = schur_complement(sp.coo_matrix(a), list(range(3)))
         assert_array_almost_equal(dtype, np.linalg.inv(s),
                                   np.linalg.inv(a)[:3, :3])
 
diff --git a/kwant/physics/__init__.py b/kwant/physics/__init__.py
index 1c463828c1133793242c7a31af708d204ef170eb..10edc954d7d06fb5223c49b57895fa51723dd77b 100644
--- a/kwant/physics/__init__.py
+++ b/kwant/physics/__init__.py
@@ -11,6 +11,6 @@
 # Merge the public interface of all submodules.
 __all__ = []
 for module in ['leads', 'dispersion', 'noise']:
-    exec 'from . import {0}'.format(module)
-    exec 'from .{0} import *'.format(module)
-    exec '__all__.extend({0}.__all__)'.format(module)
+    exec('from . import {0}'.format(module))
+    exec('from .{0} import *'.format(module))
+    exec('__all__.extend({0}.__all__)'.format(module))
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index 586eaaed21ecd620400334a0ffe0351218a74a72..edacd0bfb58110f1fc40cb9d7a4c46ce42350094 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -6,10 +6,10 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
+
 from math import sin, cos, sqrt, pi, copysign
 from collections import namedtuple
-from itertools import izip
+
 import numpy as np
 import numpy.linalg as npl
 import scipy.linalg as la
@@ -434,7 +434,7 @@ def make_proper_modes(lmbdainv, psi, extract, tol=1e6):
     if boundaries.shape == (0,) and len(angles):
         boundaries = np.array([0, len(angles)])
 
-    for interval in izip(boundaries[:-1], boundaries[1:]):
+    for interval in zip(boundaries[:-1], boundaries[1:]):
         if interval[1] > boundaries[0] + len(angles):
             break
 
@@ -661,8 +661,8 @@ def square_selfenergy(width, hopping, fermi_energy):
     psi_p_i = np.empty((width, width))
     factor = pi / (width + 1)
     prefactor = sqrt(2 / (width + 1))
-    for p in xrange(width):
-        for i in xrange(width):
+    for p in range(width):
+        for i in range(width):
             psi_p_i[p, i] = prefactor * sin(factor * (p + 1) * (i + 1))
 
     # Precalculate the integrals of the longitudinal wave functions.
@@ -672,15 +672,15 @@ def square_selfenergy(width, hopping, fermi_energy):
         else:
             return q/2 - copysign(sqrt((q / 2) ** 2 - 1), q)
     f_p = np.empty((width,), dtype=complex)
-    for p in xrange(width):
+    for p in range(width):
         e = 2 * hopping * (1 - cos(factor * (p + 1)))
         q = (fermi_energy - e) / hopping - 2
         f_p[p] = f(q)
 
     # Put everything together into the self energy and return it.
     result = np.empty((width, width), dtype=complex)
-    for i in xrange(width):
-        for j in xrange(width):
+    for i in range(width):
+        for j in range(width):
             result[i, j] = hopping * (
                 psi_p_i[:, i] * psi_p_i[:, j].conj() * f_p).sum()
     return result
diff --git a/kwant/physics/tests/test_leads.py b/kwant/physics/tests/test_leads.py
index fb7a5b33c7d4866128d5a46a08b958ce96b3313b..15da9dc3d1e92e9e3be7c6065cd7a435b63624db 100644
--- a/kwant/physics/tests/test_leads.py
+++ b/kwant/physics/tests/test_leads.py
@@ -6,9 +6,9 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
+
 import numpy as np
-from itertools import product, izip
+from itertools import product
 from numpy.testing import assert_almost_equal
 from kwant.physics import leads
 import kwant
@@ -273,14 +273,14 @@ def test_algorithm_equivalence():
         evan_vecs.append(full_vecs[:, 2 * result.nmodes :])
 
     msg = 'Stabilization {0} failed.'
-    for vecs, algo in izip(prop_vecs, algos):
+    for vecs, algo in zip(prop_vecs, algos):
         # Propagating modes should have identical ordering, and only vary
         # By a phase
         np.testing.assert_allclose(np.abs(np.sum(vecs/prop_vecs[0],
                                                  axis=0)), vecs.shape[0],
                                    err_msg=msg.format(algo))
 
-    for vecs, algo in izip(evan_vecs, algos):
+    for vecs, algo in zip(evan_vecs, algos):
         # Evanescent modes must span the same linear space.
         assert (np.linalg.matrix_rank(np.c_[vecs, evan_vecs[0]], tol=1e-12) ==
                 vecs.shape[1]), msg.format(algo)
diff --git a/kwant/plotter.py b/kwant/plotter.py
index 42ba2bfb99fc4d3dbc8b40ae9434c1526bb4f5c9..22bd311c1b4db700653ba67c02f6f278b4834871 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -69,7 +69,7 @@ if mpl_enabled:
 # Collections that allow for symbols and linewiths to be given in data space
 # (not for general use, only implement what's needed for plotter)
 def isarray(var):
-    if hasattr(var, '__getitem__') and not isinstance(var, basestring):
+    if hasattr(var, '__getitem__') and not isinstance(var, str):
         return True
     else:
         return False
@@ -717,18 +717,18 @@ def sys_leads_sites(sys, num_lead_cells=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_cells)))
+                              range(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)]
+        sites = [(i, None, 0) for i in range(sys.graph.num_nodes)]
         for leadnr, lead in enumerate(sys.leads):
             start = len(sites)
             # We will only plot leads with a graph and with a symmetry.
             if (hasattr(lead, 'graph') and hasattr(lead, 'symmetry') and
                 len(sys.lead_interfaces[leadnr])):
                 sites.extend(((site, leadnr, i) for site in
-                              xrange(lead.cell_size) for i in
-                              xrange(num_lead_cells)))
+                              range(lead.cell_size) for i in
+                              range(num_lead_cells)))
             lead_cells.append(slice(start, len(sites)))
     else:
         raise TypeError('Unrecognized system type.')
@@ -796,10 +796,10 @@ def sys_leads_pos(sys, site_lead_nr):
         # Conversion to numpy array here useful for efficiency
         vec = np.array(sym.periods)[0]
         return vec, dom
-    vecs_doms = dict((i, get_vec_domain(i)) for i in xrange(len(sys.leads)))
+    vecs_doms = dict((i, get_vec_domain(i)) for i in range(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_cells)]
+    for k, v in vecs_doms.items():
+        vecs_doms[k] = [v[0] * i for i in range(v[1], v[1] + num_lead_cells)]
     pos += [vecs_doms[i[1]][i[2]] for i in site_lead_nr]
     return pos
 
@@ -856,11 +856,11 @@ def sys_leads_hoppings(sys, num_lead_cells=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_cells)))
+                                 range(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):
+            for i in range(sys.graph.num_nodes):
                 for j in sys.graph.out_neighbors(i):
                     if i < j:
                         yield i, j
@@ -872,7 +872,7 @@ def sys_leads_hoppings(sys, num_lead_cells=2):
             if (hasattr(lead, 'graph') and hasattr(lead, 'symmetry') and
                 len(sys.lead_interfaces[leadnr])):
                 hoppings.extend(((hop, leadnr, i) for hop in ll_hoppings(lead)
-                                 for i in xrange(num_lead_cells)))
+                                 for i in range(num_lead_cells)))
             lead_cells.append(slice(start, len(hoppings)))
     else:
         raise TypeError('Unrecognized system type.')
@@ -942,12 +942,12 @@ def sys_leads_hopping_pos(sys, hop_lead_nr):
         vec = np.array(sym.periods)[0]
         return np.r_[vec, vec], dom
 
-    vecs_doms = dict((i, get_vec_domain(i)) for i in xrange(len(sys.leads)))
+    vecs_doms = dict((i, get_vec_domain(i)) for i in range(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_cells)]
+    for k, v in vecs_doms.items():
+        vecs_doms[k] = [v[0] * i for i in range(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:])
+    return np.copy(pos[:, : dim // 2]), np.copy(pos[:, dim // 2:])
 
 
 # Useful plot functions (to be extended).
@@ -1136,7 +1136,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
         if name in ('site_size', 'site_lw') and isinstance(value, tuple):
             raise TypeError('{0} may not be a tuple, use list or '
                             'array instead.'.format(name))
-        if isinstance(value, (basestring, tuple)):
+        if isinstance(value, (str, tuple)):
             return
         try:
             if len(value) != n_sys_sites:
@@ -1226,7 +1226,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
         for i, symbol in enumerate(site_symbol):
             symbol_dict[symbol].append(i)
         symbol_slcs = []
-        for symbol, indx in symbol_dict.iteritems():
+        for symbol, indx in symbol_dict.items():
             symbol_slcs.append((symbol, np.array(indx)))
         fancy_indexing = True
     else:
@@ -1285,7 +1285,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
                        else defaults['hop_lw'][dim])
 
     hop_cmap = None
-    if not isinstance(cmap, basestring):
+    if not isinstance(cmap, str):
         try:
             cmap, hop_cmap = cmap
         except TypeError:
diff --git a/kwant/rmt.py b/kwant/rmt.py
index cff30d4fc381884c6092f410d42bc723e5ace17a..6442fc0a90e220277086f35a9a07576f327e4299 100644
--- a/kwant/rmt.py
+++ b/kwant/rmt.py
@@ -6,8 +6,6 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
-
 __all__ = ['gaussian', 'circular']
 
 import numpy as np
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index d74a9473a225c6d2b41ab2c6d5cf481047931c22..0d59ce553d97f0620b9e2425565e7e91c7acf1b1 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -15,6 +15,7 @@ import numpy as np
 import scipy.sparse as sp
 from .._common import ensure_isinstance
 from .. import system
+from functools import reduce
 
 # Currently, scipy.sparse does not support matrices with one dimension being
 # zero: http://projects.scipy.org/scipy/ticket/1602 We use NumPy dense matrices
@@ -28,7 +29,7 @@ from .. import system
 LinearSys = namedtuple('LinearSys', ['lhs', 'rhs', 'indices', 'num_orb'])
 
 
-class SparseSolver(object):
+class SparseSolver(object, metaclass=abc.ABCMeta):
     """Solver class for computing physical quantities based on solving
     a liner system of equations.
 
@@ -51,7 +52,6 @@ class SparseSolver(object):
       too avoid excessive memory usage, but for some solvers not too small for
       performance reasons.
     """
-    __metaclass__ = abc.ABCMeta
 
     @abc.abstractmethod
     def _factorized(self, a):
@@ -333,11 +333,11 @@ class SparseSolver(object):
 
         n = len(sys.lead_interfaces)
         if in_leads is None:
-            in_leads = range(n)
+            in_leads = list(range(n))
         else:
             in_leads = list(in_leads)
         if out_leads is None:
-            out_leads = range(n)
+            out_leads = list(range(n))
         else:
             out_leads = list(out_leads)
         if (np.any(np.diff(in_leads) <= 0) or np.any(np.diff(out_leads) <= 0)):
@@ -420,11 +420,11 @@ class SparseSolver(object):
 
         n = len(sys.lead_interfaces)
         if in_leads is None:
-            in_leads = range(n)
+            in_leads = list(range(n))
         else:
             in_leads = list(in_leads)
         if out_leads is None:
-            out_leads = range(n)
+            out_leads = list(range(n))
         else:
             out_leads = list(out_leads)
         if (np.any(np.diff(in_leads) <= 0) or np.any(np.diff(out_leads) <= 0)):
@@ -489,7 +489,7 @@ class SparseSolver(object):
                 raise NotImplementedError("ldos for leads with only "
                                           "self-energy is not implemented yet.")
 
-        linsys = self._make_linear_sys(sys, xrange(len(sys.leads)), energy,
+        linsys = self._make_linear_sys(sys, range(len(sys.leads)), energy,
                                        args, check_hermiticity)[0]
 
         ldos = np.zeros(linsys.num_orb, float)
@@ -503,7 +503,7 @@ class SparseSolver(object):
         # See comment about zero-shaped sparse matrices at the top of common.py.
         rhs = sp.bmat([[i for i in linsys.rhs if i.shape[1]]],
                       format=self.rhsformat)
-        for j in xrange(0, rhs.shape[1], self.nrhs):
+        for j in range(0, rhs.shape[1], self.nrhs):
             jend = min(j + self.nrhs, rhs.shape[1])
             psi = self._solve_linear_sys(factored, rhs[:, j:jend],
                                          slice(linsys.num_orb))
@@ -555,7 +555,7 @@ class WaveFunction(object):
                 msg = ('Wave functions for leads with only self-energy'
                        ' are not available yet.')
                 raise NotImplementedError(msg)
-        linsys = solver._make_linear_sys(sys, xrange(len(sys.leads)), energy,
+        linsys = solver._make_linear_sys(sys, range(len(sys.leads)), energy,
                                          args, check_hermiticity)[0]
         self.solve = solver._solve_linear_sys
         self.rhs = linsys.rhs
@@ -568,11 +568,10 @@ class WaveFunction(object):
         return result.transpose()
 
 
-class BlockResult(object):
+class BlockResult(object, metaclass=abc.ABCMeta):
     """
     ABC for a linear system solution with variable grouping.
     """
-    __metaclass__ = abc.ABCMeta
 
     def __init__(self, data, lead_info, out_leads, in_leads, sizes,
                  current_conserving=False):
@@ -676,7 +675,7 @@ class BlockResult(object):
         Section 2.4 of the book by S. Datta.
         """
         n = len(self.lead_info)
-        rn = xrange(n)
+        rn = range(n)
         result = np.array([[-self.transmission(i, j) if i != j else 0
                             for j in rn] for i in rn])
         # Set the diagonal elements such that the sums of columns are zero.
diff --git a/kwant/solvers/mumps.py b/kwant/solvers/mumps.py
index fc4b0e8a14f12b96dc320c6a81094e28c81891c1..f0b9ac0c5286e462dff87ca48c15e52efd6bfdcf 100644
--- a/kwant/solvers/mumps.py
+++ b/kwant/solvers/mumps.py
@@ -84,14 +84,14 @@ class Solver(common.SparseSolver):
             self.nrhs = nrhs
 
         if ordering is not None:
-            if ordering not in mumps.orderings.keys() + ['kwant_decides']:
-                raise ValueError("Invalid ordering: " + ordering)
             if ordering == 'kwant_decides':
                 # Choose what is considered to be the best ordering.
                 sorted_orderings = [order
                                     for order in ['metis', 'scotch', 'auto']
                                     if order in mumps.possible_orderings()]
                 ordering = sorted_orderings[0]
+            elif ordering not in mumps.orderings:
+                raise ValueError("Invalid ordering: " + ordering)
             self.ordering = ordering
 
         if sparse_rhs is not None:
@@ -111,7 +111,7 @@ class Solver(common.SparseSolver):
         solve = factorized_a.solve
         sols = []
 
-        for j in xrange(0, b.shape[1], self.nrhs):
+        for j in range(0, b.shape[1], self.nrhs):
             tmprhs = b[:, j:min(j + self.nrhs, b.shape[1])]
 
             if not self.sparse_rhs:
diff --git a/kwant/solvers/sparse.py b/kwant/solvers/sparse.py
index 044ebc269bab15677e4d89d786469ef5e7a7f3e0..69f9a3a83534eb37885c6488d543c35edcaa55fc 100644
--- a/kwant/solvers/sparse.py
+++ b/kwant/solvers/sparse.py
@@ -110,7 +110,7 @@ class Solver(common.SparseSolver):
 
         sols = []
         vec = np.empty(b.shape[0], complex)
-        for j in xrange(b.shape[1]):
+        for j in range(b.shape[1]):
             vec[:] = b[:, j].todense().flatten()
             sols.append(factorized_a(vec)[kept_vars])
 
diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py
index 5eaa20b1d4af9b2e131efe5e29d5eca2fc5514f6..3c26d263be9fae52e85fba704318e91331bf5130 100644
--- a/kwant/solvers/tests/_test_sparse.py
+++ b/kwant/solvers/tests/_test_sparse.py
@@ -6,7 +6,7 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
+
 from math import cos, sin
 import numpy as np
 from nose.tools import assert_raises
@@ -263,18 +263,18 @@ def test_tricky_singular_hopping(smatrix):
     lead = kwant.Builder(kwant.TranslationalSymmetry((4, 0)))
 
     interface = []
-    for i in xrange(n):
+    for i in range(n):
         site = sq(-1, i)
         interface.append(site)
         system[site] = 0
-        for j in xrange(4):
+        for j in range(4):
             lead[sq(j, i)] = 0
-    for i in xrange(n-1):
+    for i in range(n-1):
         system[sq(-1, i), sq(-1, i+1)] = -1
-        for j in xrange(4):
+        for j in range(4):
             lead[sq(j, i), sq(j, i+1)] = -1
-    for i in xrange(n):
-        for j in xrange(4):
+    for i in range(n):
+        for j in range(4):
             lead[sq(j, i), sq(j+1, i)] = -1
     del lead[sq(1, 0), sq(2, 0)]
 
@@ -300,10 +300,10 @@ def test_many_leads(*factories):
 
     # Build a square system with four leads and a hole in the center.
     sys = kwant.Builder()
-    sys[(sq(x, y) for x in xrange(-4, 4) for y in xrange(-4, 4)
+    sys[(sq(x, y) for x in range(-4, 4) for y in range(-4, 4)
          if x**2 + y**2 >= 2)] = 3
     sys[sq.neighbors()] = phase
-    for r in [xrange(-4, -1), xrange(4)]:
+    for r in [range(-4, -1), range(4)]:
         lead = kwant.Builder(kwant.TranslationalSymmetry([-1, 0]))
         lead[(sq(0, y) for y in r)] = 4
         lead[sq.neighbors()] = phase
@@ -311,7 +311,7 @@ def test_many_leads(*factories):
         sys.attach_lead(lead.reversed())
     sys = sys.finalized()
 
-    r4 = range(4)
+    r4 = list(range(4))
     br = factories[0](sys, E, out_leads=r4, in_leads=r4)
     trans = np.array([[br._transmission(i, j) for j in r4] for i in r4])
     cmat = br.conductance_matrix()
@@ -480,7 +480,7 @@ def test_wavefunc_ldos_consistency(wave_function, ldos):
         for energy in [0, 1000]:
             wf = wave_function(sys, energy)
             ldos2 = np.zeros(wf.num_orb, float)
-            for lead in xrange(len(sys.leads)):
+            for lead in range(len(sys.leads)):
                 temp = abs(wf(lead))
                 temp **= 2
                 ldos2 += temp.sum(axis=0)
diff --git a/kwant/solvers/tests/test_mumps.py b/kwant/solvers/tests/test_mumps.py
index a185a04c64f63b862ddb7de9827177fd3221a2d8..e07aa120162d69304f405eaeef9443ff60a050fd 100644
--- a/kwant/solvers/tests/test_mumps.py
+++ b/kwant/solvers/tests/test_mumps.py
@@ -11,7 +11,7 @@ from numpy.testing.decorators import skipif
 try:
     from kwant.solvers.mumps import (
         smatrix, greens_function, ldos, wave_function, options, reset_options)
-    import _test_sparse
+    from . import _test_sparse
     no_mumps = False
 except ImportError:
     no_mumps = True
diff --git a/kwant/solvers/tests/test_sparse.py b/kwant/solvers/tests/test_sparse.py
index c9920f9d8acd42f9bbd162788ea9cef9a30e429a..f93d7771f9c88f3d287539dec47f647e4f32f7e7 100644
--- a/kwant/solvers/tests/test_sparse.py
+++ b/kwant/solvers/tests/test_sparse.py
@@ -9,7 +9,7 @@
 from nose.plugins.skip import SkipTest
 from  kwant.solvers.sparse import smatrix, greens_function, ldos, wave_function
 import kwant.solvers.sparse
-import _test_sparse
+from . import _test_sparse
 
 def test_output():
     _test_sparse.test_output(smatrix)
diff --git a/kwant/system.py b/kwant/system.py
index a323eebe4f36882c746d336e6735f02e1f11cf0c..ca541b4ad9c832d3971b75fe61afe707b7e27230 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -8,7 +8,6 @@
 
 """Low-level interface of systems"""
 
-from __future__ import division
 __all__ = ['System', 'FiniteSystem', 'InfiniteSystem']
 
 import abc
@@ -17,7 +16,7 @@ from copy import copy
 from . import _system
 
 
-class System(object):
+class System(object, metaclass=abc.ABCMeta):
     """Abstract general low-level system.
 
     Attributes
@@ -33,7 +32,6 @@ class System(object):
     Optionally, a class derived from `System` can provide a method `pos` which
     is assumed to return the real-space position of a site given its index.
     """
-    __metaclass__ = abc.ABCMeta
 
     @abc.abstractmethod
     def hamiltonian(self, i, j, *args):
@@ -49,11 +47,10 @@ class System(object):
         pass
 
 # Add a C-implemented function as an unbound method to class System.
-System.hamiltonian_submatrix = types.MethodType(
-    _system.hamiltonian_submatrix, None, System)
+System.hamiltonian_submatrix = _system.HamiltonianSubmatrix()
 
 
-class FiniteSystem(System):
+class FiniteSystem(System, metaclass=abc.ABCMeta):
     """Abstract finite low-level system, possibly with leads.
 
     Attributes
@@ -80,7 +77,6 @@ class FiniteSystem(System):
     the first ``len(lead_interfaces[n])`` sites of the InfiniteSystem.
 
     """
-    __metaclass__ = abc.ABCMeta
 
     def precalculate(self, energy=0, args=(), leads=None,
                      what='modes'):
@@ -123,7 +119,7 @@ class FiniteSystem(System):
 
         result = copy(self)
         if leads is None:
-            leads = range(len(self.leads))
+            leads = list(range(len(self.leads)))
         new_leads = []
         for nr, lead in enumerate(self.leads):
             if nr not in leads:
@@ -141,7 +137,7 @@ class FiniteSystem(System):
         result.leads = new_leads
         return result
 
-class InfiniteSystem(System):
+class InfiniteSystem(System, metaclass=abc.ABCMeta):
     """
     Abstract infinite low-level system.
 
@@ -183,18 +179,17 @@ class InfiniteSystem(System):
     infinite system.  The other scheme has the numbers of site 0 and 1
     exchanged, as well as of site 3 and 4.
     """
-    __metaclass__ = abc.ABCMeta
 
     def cell_hamiltonian(self, args=(), sparse=False):
         """Hamiltonian of a single cell of the infinite system."""
-        cell_sites = xrange(self.cell_size)
+        cell_sites = range(self.cell_size)
         return self.hamiltonian_submatrix(args, cell_sites, cell_sites,
                                           sparse=sparse)
 
     def inter_cell_hopping(self, args=(), sparse=False):
         """Hopping Hamiltonian between two cells of the infinite system."""
-        cell_sites = xrange(self.cell_size)
-        interface_sites = xrange(self.cell_size, self.graph.num_nodes)
+        cell_sites = range(self.cell_size)
+        interface_sites = range(self.cell_size, self.graph.num_nodes)
         return self.hamiltonian_submatrix(args, cell_sites, interface_sites,
                                           sparse=sparse)
 
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index d60d0884e4ac06c54cbaef9b68ed6331b6327a86..11b95c4b3be0d1e0572ffee7f65d5eb8892a60f2 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -6,7 +6,6 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
 import warnings
 from random import Random
 from nose.tools import assert_raises
@@ -84,7 +83,7 @@ def test_bad_keys():
                     try:
                         assert_raises(error, func, key)
                     except AssertionError:
-                        print func, error, key
+                        print(func, error, key)
                         raise
 
 
@@ -259,14 +258,14 @@ def random_hopping_integral(rng):
 
 def check_onsite(fsys, sites, subset=False, check_values=True):
     freq = {}
-    for node in xrange(fsys.graph.num_nodes):
+    for node in range(fsys.graph.num_nodes):
         site = fsys.sites[node].tag
         freq[site] = freq.get(site, 0) + 1
         if check_values and site in sites:
             assert fsys.onsite_hamiltonians[node] is sites[site]
     if not subset:
         # Check that all sites of `fsys` are in `sites`.
-        for site in freq.iterkeys():
+        for site in freq.keys():
             assert site in sites
     # Check that all sites of `sites` are in `fsys`.
     for site in sites:
@@ -302,7 +301,7 @@ def test_finalization():
 
     def set_hops(dest, sites):
         while len(dest) < n_hops:
-            a, b = rng.sample(sites, 2)
+            a, b = rng.sample(list(sites), 2)
             if (a, b) not in dest and (b, a) not in dest:
                 dest[a, b] = random_hopping_integral(rng)
 
@@ -327,7 +326,7 @@ def test_finalization():
     set_hops(lead_hops, lead_sites)
     lead_sites_list = list(lead_sites)
     neighbors = set()
-    for i in xrange(n_hops):
+    for i in range(n_hops):
         while True:
             a = rng.choice(lead_sites_list)
             b = rng.choice(possible_neighbors)
@@ -343,9 +342,9 @@ def test_finalization():
     # Build scattering region from blueprint and test it.
     sys = builder.Builder()
     fam = kwant.lattice.general(ta.identity(2))
-    for site, value in sr_sites.iteritems():
+    for site, value in sr_sites.items():
         sys[fam(*site)] = value
-    for hop, value in sr_hops.iteritems():
+    for hop, value in sr_hops.items():
         sys[fam(*hop[0]), fam(*hop[1])] = value
     fsys = sys.finalized()
     check_onsite(fsys, sr_sites)
@@ -353,7 +352,7 @@ def test_finalization():
 
     # Build lead from blueprint and test it.
     lead = builder.Builder(kwant.TranslationalSymmetry((size, 0)))
-    for site, value in lead_sites.iteritems():
+    for site, value in lead_sites.items():
         shift = rng.randrange(-5, 6) * size
         site = site[0] + shift, site[1]
         lead[fam(*site)] = value
@@ -363,7 +362,7 @@ def test_finalization():
         assert_equal(len(w), 1)
         assert issubclass(w[0].category, RuntimeWarning)
         assert "disconnected" in str(w[0].message)
-    for (a, b), value in lead_hops.iteritems():
+    for (a, b), value in lead_hops.items():
         shift = rng.randrange(-5, 6) * size
         a = a[0] + shift, a[1]
         b = b[0] + shift, b[1]
@@ -645,21 +644,21 @@ def test_ModesLead_and_SelfEnergyLead():
     energies = [0.9, 1.7]
 
     sys = builder.Builder()
-    for x in xrange(L):
-        for y in xrange(L):
+    for x in range(L):
+        for y in range(L):
             sys[lat(x, y)] = 4 * t + rng.random() - 0.5
     sys[hoppings] = -t
 
     # Attach a lead from the left.
     lead = builder.Builder(VerySimpleSymmetry(-1))
-    for y in xrange(L):
+    for y in range(L):
         lead[lat(0, y)] = 4 * t
     lead[hoppings] = -t
     sys.attach_lead(lead)
 
     # Make the right lead and attach it.
     lead = builder.Builder(VerySimpleSymmetry(1))
-    for y in xrange(L):
+    for y in range(L):
         lead[lat(0, y)] = 4 * t
     lead[hoppings] = -t
     sys.attach_lead(lead)
@@ -669,7 +668,7 @@ def test_ModesLead_and_SelfEnergyLead():
 
     # Replace lead with it's finalized copy.
     lead = fsys.leads[1]
-    interface = [lat(L-1, lead.sites[i].tag[1]) for i in xrange(L)]
+    interface = [lat(L-1, lead.sites[i].tag[1]) for i in range(L)]
 
     # Re-attach right lead as ModesLead.
     sys.leads[1] = builder.ModesLead(lead.modes, interface)
diff --git a/kwant/tests/test_lattice.py b/kwant/tests/test_lattice.py
index 5b389d370d80289c724ec83c8d1b02389f5e579a..f2169863a75a48aa91fc6c7c2d3229433fec08ba 100644
--- a/kwant/tests/test_lattice.py
+++ b/kwant/tests/test_lattice.py
@@ -6,7 +6,6 @@
 # the file AUTHORS.rst at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-from __future__ import division
 from math import sqrt
 import numpy as np
 import tinyarray as ta
@@ -63,15 +62,15 @@ def test_shape():
     sites = list(lat.shape(in_circle, (0, 0))())
     sites_alt = list()
     sl0, sl1 = lat.sublattices
-    for x in xrange(-2, 3):
-        for y in xrange(-2, 3):
+    for x in range(-2, 3):
+        for y in range(-2, 3):
             tag = (x, y)
             for site in (sl0(*tag), sl1(*tag)):
                 if in_circle(site.pos):
                     sites_alt.append(site)
     assert len(sites) == len(sites_alt)
     assert_equal(set(sites), set(sites_alt))
-    assert_raises(ValueError, lat.shape(in_circle, (10, 10))().next)
+    assert_raises(ValueError, lat.shape(in_circle, (10, 10))().__next__)
 
     # Check if narrow ribbons work.
     for period in (0, 1), (1, 0), (1, -1):
diff --git a/kwant/tests/test_plotter.py b/kwant/tests/test_plotter.py
index d3074b0cb99e46690ed15f84cb18ba1e5ef51a05..e34d1aeb14c7d647817850eb30e408ca94e0e806 100644
--- a/kwant/tests/test_plotter.py
+++ b/kwant/tests/test_plotter.py
@@ -31,7 +31,7 @@ def test_importable_without_matplotlib():
 
     with warnings.catch_warnings(record=True) as w:
         warnings.simplefilter("always")
-        exec code               # Trigger the warning.
+        exec(code)               # Trigger the warning.
         nose.tools.assert_equal(len(w), 1)
         assert issubclass(w[0].category, RuntimeWarning)
         assert "only iterator-providing functions" in str(w[0].message)
@@ -103,7 +103,7 @@ def test_plot():
             for sys in (sys2d, sys3d):
                 fig = plot(sys, site_color=color, cmap='binary', file=out)
                 if (color != 'k' and
-                    isinstance(color(iter(sys2d.sites()).next()), float)):
+                    isinstance(color(next(iter(sys2d.sites()))), float)):
                     assert fig.axes[0].collections[0].get_array() is not None
                 assert len(fig.axes[0].collections) == (8 if sys is sys2d else
                                                         6)
@@ -114,7 +114,7 @@ def test_plot():
             for sys in (sys2d, sys3d):
                 fig = plot(sys2d, hop_color=color, cmap='binary', file=out,
                            fig_size=(2, 10), dpi=30)
-                if color != 'k' and isinstance(color(iter(sys2d.sites()).next(),
+                if color != 'k' and isinstance(color(next(iter(sys2d.sites())),
                                                           None), float):
                     assert fig.axes[0].collections[1].get_array() is not None
 
@@ -150,10 +150,10 @@ def test_map():
                                  pos_transform=bad_transform, file=out)
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
-            plotter.map(sys.finalized(), xrange(len(sys.sites())),
+            plotter.map(sys.finalized(), range(len(sys.sites())),
                               file=out)
         nose.tools.assert_raises(ValueError, plotter.map, sys,
-                                 xrange(len(sys.sites())), file=out)
+                                 range(len(sys.sites())), file=out)
 
 
 def test_mask_interpolate():
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index 29af26c2eb6e6043141536fe3605777af84a5a7f..3eb1e15b6100364bd27ff5d6712a35bb9dd010a4 100644
--- a/kwant/tests/test_system.py
+++ b/kwant/tests/test_system.py
@@ -14,9 +14,9 @@ import kwant
 def test_hamiltonian_submatrix():
     sys = kwant.Builder()
     chain = kwant.lattice.chain()
-    for i in xrange(3):
+    for i in range(3):
         sys[chain(i)] = 0.5 * i
-    for i in xrange(2):
+    for i in range(2):
         sys[chain(i), chain(i + 1)] = 1j * (i + 1)
 
     sys2 = sys.finalized()
@@ -87,8 +87,8 @@ def test_hamiltonian_submatrix():
         return p1 - p2
 
     sys = kwant.Builder()
-    sys[(chain(i) for i in xrange(3))] = onsite
-    sys[((chain(i), chain(i + 1)) for i in xrange(2))] = hopping
+    sys[(chain(i) for i in range(3))] = onsite
+    sys[((chain(i), chain(i + 1)) for i in range(2))] = hopping
     sys2 = sys.finalized()
     mat = sys2.hamiltonian_submatrix((2, 1))
     mat_should_be = [[3, 1, 0], [1, 4, 1], [0, 1, 5]]
diff --git a/setup.py b/setup.py
index 3b884b97e239eedb01950de25217d13d3deb39f4..2d14a51675cd70641d6fe87a7ef2e806a6282b02 100755
--- a/setup.py
+++ b/setup.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 
 # Copyright 2011-2015 Kwant authors.
 #
@@ -16,7 +16,7 @@ import os
 import glob
 import imp
 import subprocess
-import ConfigParser
+import configparser
 from setuptools import setup, find_packages, Extension, Command
 from sysconfig import get_platform
 from distutils.errors import DistutilsError, DistutilsModuleError, \
@@ -42,16 +42,23 @@ CYTHON_OPTION = '--cython'
 TUT_DIR = 'tutorial'
 TUT_GLOB = 'doc/source/tutorial/*.py'
 TUT_HIDDEN_PREFIX = '#HIDDEN'
+CLASSIFIERS = """\
+    Development Status :: 5 - Production/Stable
+    Intended Audience :: Science/Research
+    Intended Audience :: Developers
+    Programming Language :: Python :: 3 :: Only
+    Topic :: Software Development
+    Topic :: Scientific/Engineering
+    Operating System :: POSIX
+    Operating System :: Unix
+    Operating System :: MacOS :: MacOS X
+    Operating System :: Microsoft :: Windows"""
+
 
 # Let Kwant itself determine its own version.  We cannot simply import kwant, as
 # it is not built yet.
 _dont_write_bytecode_saved = sys.dont_write_bytecode
 sys.dont_write_bytecode = True
-try:
-    imp.load_source(STATIC_VERSION_PATH[-1].split('.')[0],
-                    os.path.join(*STATIC_VERSION_PATH))
-except IOError:
-    pass
 _common = imp.load_source('_common', 'kwant/_common.py')
 sys.dont_write_bytecode = _dont_write_bytecode_saved
 
@@ -166,7 +173,7 @@ def git_lsfiles():
 
     if p.wait() != 0:
         return
-    return p.communicate()[0].split('\n')[:-1]
+    return p.communicate()[0].decode().split('\n')[:-1]
 
 
 # Make the command "sdist" depend on "build".  This verifies that the
@@ -263,7 +270,7 @@ def search_mumps():
     except OSError:
         pass
     else:
-        p.communicate(input='int main() {}\n')
+        p.communicate(input=b'int main() {}\n')
         if p.wait() == 0:
             return {'libraries': libs}
     return {}
@@ -297,7 +304,7 @@ def extensions():
                       'kwant/graph/c_slicer/slicer.h']})]
 
     #### Add components of Kwant with external compile-time dependencies.
-    config = ConfigParser.ConfigParser()
+    config = configparser.ConfigParser()
     try:
         with open(CONFIG_FILE) as f:
             config.readfp(f)
@@ -334,7 +341,7 @@ def extensions():
         if kwrds:
             build_summary.append('Auto-configured MUMPS')
     if kwrds:
-        for name, value in lapack.iteritems():
+        for name, value in lapack.items():
             kwrds.setdefault(name, []).extend(value)
         kwrds.setdefault('depends', []).extend(
             [CONFIG_FILE, 'kwant/linalg/cmumps.pxd'])
@@ -467,7 +474,8 @@ def main():
           include_dirs=include_dirs,
           setup_requires=['numpy > 1.6.1', 'nose >= 1.0'],
           install_requires=['numpy > 1.6.1', 'scipy >= 0.9', 'tinyarray'],
-          extras_require={'plotting': 'matplotlib >= 1.2'})
+          extras_require={'plotting': 'matplotlib >= 1.2'},
+          classifiers=[c.strip() for c in CLASSIFIERS.split('\n')])
 
 if __name__ == '__main__':
     main()