diff --git a/kwant/builder.py b/kwant/builder.py
index 93fcdc236e0fa6101d259fe676ffb8bd75f01b7f..22afb119db1edf023d1c628e1a8ab9145660a2f1 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -3,7 +3,8 @@ from __future__ import division
 __all__ = ['Builder', 'Site', 'SiteGroup', 'SimpleSiteGroup', 'Symmetry',
            'Lead', 'BuilderLead', 'SelfEnergy']
 
-import abc, sys
+import abc
+import sys
 import operator
 from itertools import izip, islice, chain
 from collections import Iterable
@@ -284,6 +285,7 @@ def herm_conj(value):
 class HermConjOfFunc(object):
     """Proxy returning the hermitian conjugate of the original result."""
     __slots__ = ('function')
+
     def __init__(self, function):
         self.function = function
 
@@ -724,8 +726,8 @@ class Builder(object):
                 # These two following lines make sure we do not waste space by
                 # storing different instances of identical sites.  They also
                 # verify that sites a and b already belong to the system.
-                a = self.H[a][0]      # Might fail.
-                b = self.H[b][0]      # Might fail.
+                a = self.H[a][0]                 # Might fail.
+                b = self.H[b][0]                 # Might fail.
                 self._set_edge(a, b, value)      # Will work.
                 self._set_edge(b, a, other)      # Will work.
             else:
@@ -733,8 +735,8 @@ class Builder(object):
                 if b2 not in self.H:
                     raise KeyError()
                 assert not sym.in_fd(a2)
-                self._set_edge(a, b, value)   # Might fail.
-                self._set_edge(b2, a2, other) # Will work.
+                self._set_edge(a, b, value)      # Might fail.
+                self._set_edge(b2, a2, other)    # Will work.
         except KeyError:
             raise KeyError(hoppinglike)
 
@@ -793,7 +795,8 @@ class Builder(object):
         sites = list(site for site in self.H
                      if self._out_degree(site) < 2)
         for site in sites:
-            if site not in self.H: continue
+            if site not in self.H:
+                continue
             while site:
                 pneighbors = tuple(self._out_neighbors(site))
                 if pneighbors:
@@ -837,14 +840,16 @@ class Builder(object):
         """
         for tail, hvhv in self.H.iteritems():
             for head, value in edges(hvhv):
-                if value is other: continue
+                if value is other:
+                    continue
                 yield (tail, head)
 
     def hopping_value_pairs(self):
         """Return an iterator over all (hopping, value) pairs."""
         for tail, hvhv in self.H.iteritems():
             for head, value in edges(hvhv):
-                if value is other: continue
+                if value is other:
+                    continue
                 yield (tail, head), value
 
     def dangling(self):
@@ -1045,10 +1050,11 @@ class Builder(object):
 
         #### Make graph.
         g = graph.Graph()
-        g.num_nodes = len(sites) # Some sites could not appear in any edge.
+        g.num_nodes = len(sites)  # Some sites could not appear in any edge.
         for tail, hvhv in self.H.iteritems():
             for head in islice(hvhv, 2, None, 2):
-                if tail == head: continue
+                if tail == head:
+                    continue
                 g.add_edge(id_by_site[tail], id_by_site[head])
         g = g.compressed()
 
@@ -1091,9 +1097,9 @@ class Builder(object):
 
         #### For each site of the fundamental domain, determine whether it has
         #### neighbors or not.
-        lsites_with = []    # Fund. domain sites with neighbors in prev. dom
-        lsites_without = [] # Remaining sites of the fundamental domain
-        for tail in self.H: # Loop over all sites of the fund. domain.
+        lsites_with = []       # Fund. domain sites with neighbors in prev. dom
+        lsites_without = []    # Remaining sites of the fundamental domain
+        for tail in self.H:    # Loop over all sites of the fund. domain.
             for head in self._out_neighbors(tail):
                 fd = sym.which(head)[0]
                 if fd == 1:
diff --git a/kwant/graph/dissection.py b/kwant/graph/dissection.py
index a368d094a179eebda016f13bb9356c8dc17ded21..0dd050e9d8b5ef90cd4a541f4779a92b80081034 100644
--- a/kwant/graph/dissection.py
+++ b/kwant/graph/dissection.py
@@ -6,6 +6,7 @@ import numpy as np
 from . import core, utils, scotch
 from .defs import gint_dtype
 
+
 def edge_dissection(gr, minimum_size, is_undirected=False):
     """Returns a nested dissection tree (represented as a nested number of
     tuples) for the graph gr, based on edge separators.
diff --git a/kwant/lattice.py b/kwant/lattice.py
index bc578039921227ce14b82af7712ba672b69e4877..070cf003a17491919f0aabd900070776420a6e9b 100644
--- a/kwant/lattice.py
+++ b/kwant/lattice.py
@@ -4,7 +4,6 @@ __all__ = ['make_lattice', 'TranslationalSymmetry',
            'PolyatomicLattice', 'MonatomicLattice']
 
 from math import sqrt
-from itertools import izip, chain
 import numpy as np
 import tinyarray as ta
 from . import builder
@@ -245,6 +244,7 @@ class TranslationalSymmetry(builder.Symmetry):
             the site group which has to be processed.  Be sure to delete the
             previously processed site groups from `site_group_data` if you want
             to modify the cache.
+
         other_vectors : list of lists of integers
             Bravais lattice vectors used to complement the periods in forming
             a basis. The fundamental domain belongs to the linear space
@@ -282,7 +282,7 @@ class TranslationalSymmetry(builder.Symmetry):
         m.T[: num_dir] = bravais_periods
         num_vec = num_dir + len(other_vectors)
         if len(other_vectors) != 0:
-            m.T[num_dir : num_vec] = other_vectors
+            m.T[num_dir:num_vec] = other_vectors
         norms = np.apply_along_axis(np.linalg.norm, 1, m)
         indices = np.argsort(norms)
         for coord in zip(indices, range(num_vec, dim)):
diff --git a/kwant/linalg/decomp_ev.py b/kwant/linalg/decomp_ev.py
index 19a5d126e11f9b5355f61ac5fd8d6c2c5e88fc7b..36ceb70612ccfc49a3fa8602b1c56d12e8e27301 100644
--- a/kwant/linalg/decomp_ev.py
+++ b/kwant/linalg/decomp_ev.py
@@ -2,6 +2,7 @@ __all__ = ['gen_eig']
 
 from . import lapack
 
+
 def gen_eig(a, b, left=False, right=True, overwrite_ab=False):
     """Compute the eigenvalues and -vectors of the matrix pencil (a,b), i.e. of
     the generalized (unsymmetric) eigenproblem a v = lambda b v where a and b
diff --git a/kwant/linalg/decomp_lu.py b/kwant/linalg/decomp_lu.py
index d7659a9d8c70a8e7d74e47caf62343dbe0505b44..32cb37ec0af1999455bdb718ad6e19d64abed08d 100644
--- a/kwant/linalg/decomp_lu.py
+++ b/kwant/linalg/decomp_lu.py
@@ -3,7 +3,8 @@ __all__ = ['lu_factor', 'lu_solve', 'rcond_from_lu']
 import numpy as np
 from . import lapack
 
-def lu_factor(a, overwrite_a = False):
+
+def lu_factor(a, overwrite_a=False):
     """Compute the LU factorization of a matrix A = P * L * U. The function
     returns a tuple (lu, p, singular), where lu contains the LU factorization
     storing the unit lower triangular matrix L in the strictly lower triangle
@@ -51,6 +52,7 @@ def lu_factor(a, overwrite_a = False):
     else:
         return lapack.cgetrf(a)
 
+
 def lu_solve((lu, ipiv, singular), b):
     """Solve a linear system of equations, a x = b, given the LU
     factorization of a
@@ -74,7 +76,7 @@ def lu_solve((lu, ipiv, singular), b):
                              "are probably unreliable")
 
     ltype, lu, b = lapack.prepare_for_lapack(False, lu, b)
-    ipiv = np.ascontiguousarray(np.asanyarray(ipiv), dtype = lapack.int_dtype)
+    ipiv = np.ascontiguousarray(np.asanyarray(ipiv), dtype=lapack.int_dtype)
 
     if b.ndim > 2:
         raise ValueError("lu_solve: b must be a vector or matrix")
@@ -91,7 +93,8 @@ def lu_solve((lu, ipiv, singular), b):
     else:
         return lapack.cgetrs(lu, ipiv, b)
 
-def rcond_from_lu((lu, ipiv, singular), norm_a, norm = "1"):
+
+def rcond_from_lu((lu, ipiv, singular), 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.
diff --git a/kwant/linalg/decomp_schur.py b/kwant/linalg/decomp_schur.py
index 96182fc8a00c35ded2e324fc5e3df2bc1788a441..01e1d3921fb6a9c7d32f751ea5afa7cace646544 100644
--- a/kwant/linalg/decomp_schur.py
+++ b/kwant/linalg/decomp_schur.py
@@ -6,6 +6,7 @@ from math import sqrt
 import numpy as np
 from . import lapack
 
+
 def schur(a, calc_q=True, calc_ev=True, overwrite_a=False):
     """Compute the Schur form of a square matrix a.
 
@@ -62,10 +63,11 @@ def schur(a, calc_q=True, calc_ev=True, overwrite_a=False):
     if a.shape[0] != a.shape[1]:
         raise ValueError("Expect square matrix")
 
-    gees = getattr(lapack, ltype+"gees")
+    gees = getattr(lapack, ltype + "gees")
 
     return gees(a, calc_q, calc_ev)
 
+
 def convert_r2c_schur(t, q):
     """Convert a real Schur form (with possibly 2x2 blocks on the diagonal)
     into a complex Schur form that is completely triangular.
@@ -120,10 +122,10 @@ def convert_r2c_schur(t, q):
         y = c
         norm = sqrt(-b * c + c * c)
 
-        U = np.array([[x/norm, -y/norm],[y/norm, -x/norm]])
+        U = np.array([[x / norm, -y / norm], [y / norm, -x / norm]])
 
         t2[i, i] = a + x
-        t2[i+1, i] = 0.0
+        t2[i+1, i] = 0
         t2[i, i+1] = -b - c
         t2[i+1, i+1] = a - x
 
@@ -184,7 +186,7 @@ def order_schur(select, t, q, calc_ev=True, overwrite_tq=False):
 
     ltype, t, q = lapack.prepare_for_lapack(overwrite_tq, t, q)
 
-    trsen = getattr(lapack, ltype+"trsen")
+    trsen = getattr(lapack, ltype + "trsen")
 
     # Figure out if select is a function or array.
     isfun = isarray = True
@@ -200,22 +202,22 @@ def order_schur(select, t, q, calc_ev=True, overwrite_tq=False):
     if not (isarray or isfun):
         raise ValueError("select must be either a function or an array")
     elif isarray:
-        select = np.array(select, dtype = lapack.logical_dtype,
-                          order = 'F')
+        select = np.array(select, dtype=lapack.logical_dtype, order='F')
     else:
         select = np.array(np.vectorize(select)(np.arange(t.shape[0])),
-                          dtype= lapack.logical_dtype, order = 'F')
+                          dtype=lapack.logical_dtype, order='F')
 
     # Now check if the reordering can actually be done as desired,
     # if we have a real Schur form (i.e. if the 2x2 blocks would be
     # separated). If this is the case, convert to complex Schur form first.
     for i in np.diagonal(t, -1).nonzero()[0]:
         if bool(select[i]) != bool(select[i+1]):
-            t, q =convert_r2c_schur(t, q)
+            t, q = convert_r2c_schur(t, q)
             return order_schur(select, t, q, calc_ev, True)
 
     return trsen(select, t, q, calc_ev)
 
+
 def evecs_from_schur(t, q, select=None, left=False, right=True,
                      overwrite_tq=False):
     """Compute eigenvectors from Schur form.
@@ -263,7 +265,7 @@ def evecs_from_schur(t, q, select=None, left=False, right=True,
         or t.shape[0] != q.shape[0]):
         raise ValueError("Invalid Schur decomposition as input")
 
-    trevc = getattr(lapack, ltype+"trevc")
+    trevc = getattr(lapack, ltype + "trevc")
 
     # check if select is a function or an array
     if select is not None:
@@ -282,16 +284,17 @@ def evecs_from_schur(t, q, select=None, left=False, right=True,
             raise ValueError("select must be either a function, "
                              "an array or None")
         elif isarray:
-            selectarr = np.array(select, dtype = lapack.logical_dtype,
-                                 order = 'F')
+            selectarr = np.array(select, dtype=lapack.logical_dtype,
+                                 order='F')
         else:
             selectarr = np.array(np.vectorize(select)(np.arange(t.shape[0])),
-                                 dtype= lapack.logical_dtype, order = 'F')
+                                 dtype=lapack.logical_dtype, order='F')
     else:
         selectarr = None
 
     return trevc(t, q, selectarr, left, right)
 
+
 def gen_schur(a, b, calc_q=True, calc_z=True, calc_ev=True,
               overwrite_ab=False):
     """Compute the generalized Schur form of a matrix pencil (a, b).
@@ -365,10 +368,11 @@ def gen_schur(a, b, calc_q=True, calc_z=True, calc_ev=True,
     if a.shape[0] != b.shape[0] or a.shape[0] != b.shape[1]:
         raise ValueError("Shape of b is incompatible to matrix a")
 
-    gges = getattr(lapack, ltype+"gges")
+    gges = getattr(lapack, ltype + "gges")
 
     return gges(a, b, calc_q, calc_z, calc_ev)
 
+
 def order_gen_schur(select, s, t, q=None, z=None, calc_ev=True,
                     overwrite_stqz=False):
     """Reorder the generalized Schur form.
@@ -442,7 +446,7 @@ def order_gen_schur(select, s, t, q=None, z=None, calc_ev=True,
                             s.shape[0] != z.shape[0]))):
         raise ValueError("Invalid Schur decomposition as input")
 
-    tgsen = getattr(lapack, ltype+"tgsen")
+    tgsen = getattr(lapack, ltype + "tgsen")
 
     # Figure out if select is a function or array.
     isfun = isarray = True
@@ -458,11 +462,10 @@ def order_gen_schur(select, s, t, q=None, z=None, calc_ev=True,
     if not (isarray or isfun):
         raise ValueError("select must be either a function or an array")
     elif isarray:
-        select = np.array(select, dtype = lapack.logical_dtype,
-                          order = 'F')
+        select = np.array(select, dtype=lapack.logical_dtype, order='F')
     else:
         select = np.array(np.vectorize(select)(np.arange(t.shape[0])),
-                          dtype= lapack.logical_dtype, order = 'F')
+                          dtype=lapack.logical_dtype, order='F')
 
     # Now check if the reordering can actually be done as desired, if we have a
     # real Schur form (i.e. if the 2x2 blocks would be separated). If this is
@@ -477,12 +480,13 @@ def order_gen_schur(select, s, t, q=None, z=None, calc_ev=True,
             elif z is not None:
                 s, t, z = convert_r2c_gen_schur(s, t, q=None, z=z)
             else:
-                s,t = convert_r2c_gen_schur(s, t)
+                s, t = convert_r2c_gen_schur(s, t)
 
             return order_gen_schur(select, s, t, q, z, calc_ev, True)
 
     return tgsen(select, s, t, q, z, calc_ev)
 
+
 def convert_r2c_gen_schur(s, t, q=None, z=None):
     """Convert a real generallzed Schur form (with possibly 2x2 blocks on the
     diagonal) into a complex Schur form that is completely triangular.  If the
@@ -563,8 +567,8 @@ def convert_r2c_gen_schur(s, t, q=None, z=None):
         # form. If necessary, order_gen_schur is used to ensure the desired
         # order of eigenvalues.
 
-        sb, tb, qb, zb, alphab, betab = gen_schur(s2[i:i+2,i:i+2],
-                                                  t2[i:i+2,i:i+2])
+        sb, tb, qb, zb, alphab, betab = gen_schur(s2[i:i+2, i:i+2],
+                                                  t2[i:i+2, i:i+2])
 
         # Ensure order of eigenvalues. (betab is positive)
         if alphab[0].imag < alphab[1].imag:
@@ -593,6 +597,7 @@ def convert_r2c_gen_schur(s, t, q=None, z=None):
     else:
         return s2, t2
 
+
 def evecs_from_gen_schur(s, t, q=None, z=None, select=None,
                          left=False, right=True, overwrite_qz=False):
     """Compute eigenvectors from Schur form.
@@ -664,7 +669,7 @@ def evecs_from_gen_schur(s, t, q=None, z=None, select=None,
     if right and z is None:
         raise ValueError("Matrix z must be provided for right eigenvectors")
 
-    tgevc = getattr(lapack, ltype+"tgevc")
+    tgevc = getattr(lapack, ltype + "tgevc")
 
     # Check if select is a function or an array.
     if select is not None:
@@ -683,11 +688,11 @@ def evecs_from_gen_schur(s, t, q=None, z=None, select=None,
             raise ValueError("select must be either a function, "
                              "an array or None")
         elif isarray:
-            selectarr = np.array(select, dtype = lapack.logical_dtype,
-                                 order = 'F')
+            selectarr = np.array(select, dtype=lapack.logical_dtype,
+                                 order='F')
         else:
             selectarr = np.array(np.vectorize(select)(np.arange(t.shape[0])),
-                                 dtype= lapack.logical_dtype, order = 'F')
+                                 dtype=lapack.logical_dtype, order='F')
     else:
         selectarr = None
 
diff --git a/kwant/physics/selfenergy.py b/kwant/physics/selfenergy.py
index eae503cfa203a885fa04cd686dc2aa0af6578300..0211fcebb02c00f64c45b51bf3181ad8ba197804 100644
--- a/kwant/physics/selfenergy.py
+++ b/kwant/physics/selfenergy.py
@@ -10,6 +10,7 @@ dot = np.dot
 
 __all__ = ['self_energy', 'modes', 'Modes']
 
+
 def setup_linsys(h_onslice, h_hop, tol=1e6):
     """
     Make an eigenvalue problem for eigenvectors of translation operator.
@@ -59,10 +60,10 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
 
         A = np.empty((2*n, 2*n), dtype=np.common_type(h_onslice, h_hop))
 
-        A[0: n, 0: n] = kla.lu_solve(sol, -h_onslice)
-        A[0: n, n: 2*n] = kla.lu_solve(sol, -h_hop.T.conj())
-        A[n: 2*n, 0: n] = np.identity(n)
-        A[n: 2*n, n: 2*n] = 0
+        A[0:n, 0:n] = kla.lu_solve(sol, -h_onslice)
+        A[0:n, n:2*n] = kla.lu_solve(sol, -h_hop.T.conj())
+        A[n:2*n, 0:n] = np.identity(n)
+        A[n:2*n, n:2*n] = 0
 
         return A
     else:
@@ -71,8 +72,8 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
 
         # Recast the svd of h_hop = u s v^dagger such that
         # u, v are matrices with shape n x n_nonsing.
-        u = u[:, : n_nonsing]
-        s = s[: n_nonsing]
+        u = u[:, :n_nonsing]
+        s = s[:n_nonsing]
         # pad v with zeros if necessary
         v = np.zeros((n, n_nonsing), dtype=vh.dtype)
         v[:vh.shape[1], :] = vh[:n_nonsing, :].T.conj()
@@ -113,7 +114,7 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
             max_h = np.amax(np.abs(h_onslice))
             max_temp = np.amax(np.abs(temp))
 
-            gamma = max_h/max_temp * 1j
+            gamma = max_h / max_temp * 1j
 
             h = h_onslice + gamma * temp
 
@@ -132,9 +133,9 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
         def extract_wf(psi, lmbdainv):
             return kla.lu_solve(sol,
                                 gamma * dot(v, psi[: n_nonsing]) +
-                                gamma * dot(u, psi[n_nonsing :] * lmbdainv) -
+                                gamma * dot(u, psi[n_nonsing:] * lmbdainv) -
                                 dot(u * s, psi[: n_nonsing] * lmbdainv) -
-                                dot(v * s, psi[n_nonsing :]))
+                                dot(v * s, psi[n_nonsing:]))
 
         # Project a full wave function back.
 
@@ -147,26 +148,27 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
         A = np.empty((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
         B = np.empty((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
 
-        A[: n_nonsing, : n_nonsing] = -np.eye(n_nonsing)
+        A[:n_nonsing, :n_nonsing] = -np.eye(n_nonsing)
 
         B[n_nonsing: 2 * n_nonsing,
           n_nonsing: 2 * n_nonsing] = np.eye(n_nonsing)
 
         temp = kla.lu_solve(sol, v)
         temp2 = dot(u.T.conj(), temp)
-        A[n_nonsing: 2 * n_nonsing, : n_nonsing] = gamma * temp2
-        A[n_nonsing: 2 * n_nonsing, n_nonsing: 2 * n_nonsing] = - temp2 * s
+        A[n_nonsing : 2 * n_nonsing, :n_nonsing] = gamma * temp2
+        A[n_nonsing : 2 * n_nonsing, n_nonsing: 2 * n_nonsing] = - temp2 * s
         temp2 = dot(v.T.conj(), temp)
-        A[: n_nonsing, : n_nonsing] += gamma * temp2
-        A[: n_nonsing, n_nonsing: 2 * n_nonsing] = - temp2 * s
+        A[:n_nonsing, :n_nonsing] += gamma * temp2
+        A[:n_nonsing, n_nonsing : 2 * n_nonsing] = - temp2 * s
 
         temp = kla.lu_solve(sol, u)
         temp2 = dot(u.T.conj(), temp)
-        B[n_nonsing:2*n_nonsing, :n_nonsing] = temp2 * s
-        B[n_nonsing:2*n_nonsing, n_nonsing:2*n_nonsing] -= gamma * temp2
+        B[n_nonsing : 2 * n_nonsing, :n_nonsing] = temp2 * s
+        B[n_nonsing : 2 * n_nonsing, n_nonsing : 2 * n_nonsing] -= \
+            gamma * temp2
         temp2 = dot(v.T.conj(), temp)
         B[:n_nonsing, :n_nonsing] = temp2 * s
-        B[:n_nonsing, n_nonsing:2*n_nonsing] = - gamma * temp2
+        B[:n_nonsing, n_nonsing : 2 * n_nonsing] = - gamma * temp2
 
         # Solving a generalized eigenproblem is about twice as expensive
         # as solving a regular eigenvalue problem.
@@ -182,10 +184,10 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
         # I put a more stringent condition here - errors can accumulate
         # from here to the eigenvalue calculation later.
         if rcond > eps * tol**2:
-            return (kla.lu_solve(lu_b, A), (u, s, v[:m,:]),
+            return (kla.lu_solve(lu_b, A), (u, s, v[:m, :]),
                     (extract_wf, project_wf))
         else:
-            return (A, B, (u, s, v[: m]), (extract_wf, project_wf))
+            return (A, B, (u, s, v[:m]), (extract_wf, project_wf))
 
 
 def split_degenerate(evs, tol=1e6):
@@ -291,7 +293,7 @@ def make_proper_modes(lmbdainv, psi, h_hop, extract=None,
             if extract is not None:
                 full_psi = extract(psi[:, indx], lmbdainv[indx])
             else:
-                full_psi = psi[: n, indx]
+                full_psi = psi[:n, indx]
 
             # Finding the true modes is done in two steps:
 
@@ -313,8 +315,8 @@ def make_proper_modes(lmbdainv, psi, h_hop, extract=None,
             if project:
                 psi[:, indx] = project(full_psi, lmbdainv[indx])
             else:
-                psi[: n, indx] = full_psi * lmbdainv[indx]
-                psi[n: 2*n, indx] = full_psi
+                psi[:n, indx] = full_psi * lmbdainv[indx]
+                psi[n:2*n, indx] = full_psi
 
             # 2. Moving infinitesimally away from the degeneracy
             # point, the modes should diagonalize the velocity
@@ -328,11 +330,11 @@ def make_proper_modes(lmbdainv, psi, h_hop, extract=None,
             # as we are happy with any superposition in this case.
 
             if h_hop.ndim == 2:
-                vel_op = -1j * dot(psi[n :, indx].T.conj(),
-                                     dot(h_hop, psi[: m, indx]))
+                vel_op = -1j * dot(psi[n:, indx].T.conj(),
+                                     dot(h_hop, psi[:m, indx]))
             else:
-                vel_op = -1j * dot(psi[n :, indx].T.conj() * h_hop,
-                                     psi[: n, indx])
+                vel_op = -1j * dot(psi[n:, indx].T.conj() * h_hop,
+                                     psi[:n, indx])
 
             vel_op = vel_op + vel_op.T.conj()
 
@@ -367,11 +369,11 @@ def make_proper_modes(lmbdainv, psi, h_hop, extract=None,
             k = indx[0]
 
             if h_hop.ndim == 2:
-                v[k] = 2 * dot(dot(psi[n: 2*n, k: k + 1].T.conj(), h_hop),
-                               psi[: m, k: k+1]).imag
+                v[k] = 2 * dot(dot(psi[n:2*n, k:k+1].T.conj(), h_hop),
+                               psi[:m, k:k+1]).imag
             else:
-                v[k] = 2 * dot(psi[n: 2*n, k: k + 1].T.conj() * h_hop,
-                               psi[0: n, k: k + 1]).imag
+                v[k] = 2 * dot(psi[n:2*n, k:k+1].T.conj() * h_hop,
+                               psi[0:n, k:k+1]).imag
 
             if v[k] > vel_eps:
                 rightselect[k] = True
@@ -459,7 +461,7 @@ def unified_eigenproblem(h_onslice, h_hop, tol):
                           eps * tol * np.abs(beta))
 
             warning_settings = np.seterr(divide='ignore', invalid='ignore')
-            ev = alpha/beta
+            ev = alpha / beta
             np.seterr(**warning_settings)
             # Note: the division is OK here, as we later only access
             #       eigenvalues close to the unit circle
@@ -571,17 +573,17 @@ def self_energy(h_onslice, h_hop, tol=1e6):
 
         if nprop > 0:
             nmodes = np.sum(rselect)
-            vecs[:, : nmodes] = prop_vecs[:, rselect]
+            vecs[:, :nmodes] = prop_vecs[:, rselect]
         else:
             nmodes = 0
 
         vecs[:, nmodes:] = vec_gen(select)
 
         if v is not None:
-            return dot(v * w, dot(vecs[n :], dot(npl.inv(vecs[: n]),
-                                                 v.T.conj())))
+            return dot(v * w, dot(vecs[n:], dot(npl.inv(vecs[:n]),
+                                                v.T.conj())))
         else:
-            return dot(h_hop.T.conj(), dot(vecs[n :], npl.inv(vecs[: n])))
+            return dot(h_hop.T.conj(), dot(vecs[n:], npl.inv(vecs[:n])))
     else:
         # Reorder all the right-going eigenmodes to the top left part of
         # the Schur decomposition.
@@ -592,13 +594,14 @@ def self_energy(h_onslice, h_hop, tol=1e6):
         z = ord_schur(select)
 
         if v is not None:
-            return dot(v * w, dot(z[n :, : n], dot(npl.inv(z[: n, : n]),
+            return dot(v * w, dot(z[n:, :n], dot(npl.inv(z[:n, :n]),
                                                    v.T.conj())))
         else:
-            return dot(h_hop.T.conj(), dot(z[n:, : n], npl.inv(z[: n, : n])))
+            return dot(h_hop.T.conj(), dot(z[n:, :n], npl.inv(z[:n, :n])))
 
 Modes = namedtuple('Modes', ['vecs', 'vecslmbdainv', 'nmodes', 'svd'])
 
+
 def modes(h_onslice, h_hop, tol=1e6):
     """
     Compute the eigendecomposition of a translation operator of a lead.
@@ -666,9 +669,9 @@ def modes(h_onslice, h_hop, tol=1e6):
 
     nprop = np.sum(propselect)
     nevan = n - nprop // 2
-    evanselect_bool = np.zeros((2 * n), dtype='bool')
+    evanselect_bool = np.zeros((2*n), dtype='bool')
     evanselect_bool[evanselect] = True
-    evan_vecs = ord_schur(evanselect)[:, : nevan]
+    evan_vecs = ord_schur(evanselect)[:, :nevan]
 
     if nprop > 0:
         # Compute the propagating eigenvectors.
@@ -697,12 +700,12 @@ def modes(h_onslice, h_hop, tol=1e6):
         nmodes = np.sum(rprop)
         vecs = np.c_[prop_vecs[n:, lprop], prop_vecs[n:, rprop],
                      evan_vecs[n:]]
-        vecslmbdainv = np.c_[prop_vecs[: n, lprop], prop_vecs[: n, rprop],
-                             evan_vecs[: n]]
+        vecslmbdainv = np.c_[prop_vecs[:n, lprop], prop_vecs[:n, rprop],
+                             evan_vecs[:n]]
 
     else:
         vecs = evan_vecs[n:]
-        vecslmbdainv = evan_vecs[: n]
+        vecslmbdainv = evan_vecs[:n]
         nmodes = 0
 
     svd = None if s is None else (u, s, v)
@@ -739,9 +742,9 @@ def square_self_energy(width, hopping, potential, fermi_energy):
     # Precalculate the integrals of the longitudinal wave functions.
     def f(q):
         if abs(q) <= 2:
-            return q/2 - 1j * sqrt(1 - (q/2)**2)
+            return q/2 - 1j * sqrt(1 - (q / 2) ** 2)
         else:
-            return q/2 - copysign(sqrt((q/2)**2 - 1), q)
+            return q/2 - copysign(sqrt((q / 2) ** 2 - 1), q)
     f_p = np.empty((width,), dtype=complex)
     for p in xrange(width):
         e = 2 * hopping * (1 - cos(factor * (p + 1)))
diff --git a/kwant/run.py b/kwant/run.py
index 0eb1ebc3a96730d112184a80cb053472e26636af..0ac97131de83b2e033d5916a8c37d1e61adb141d 100644
--- a/kwant/run.py
+++ b/kwant/run.py
@@ -1,18 +1,22 @@
 """Support for running functions from the command line"""
 
 from __future__ import division
-import os, struct, sys
-import numpy, scipy
+import os
+import struct
+import sys
+import numpy
+import scipy
 from .version import version
 __all__ = ['randomize', 'exec_argv']
 
 numpy_version = numpy.version.version
 if not numpy.version.release:
-    numpy_version +=  '-non-release'
+    numpy_version += '-non-release'
 
 scipy_version = scipy.version.version
 if not scipy.version.release:
-    scipy_version +=  '-non-release'
+    scipy_version += '-non-release'
+
 
 def randomize():
     """Seed numpy's random generator according to RNG_SEED environment
@@ -36,13 +40,14 @@ def randomize():
 
     # numpy.random.seed only uses the lower 32 bits of the seed argument, so we
     # split our 64 bit seed into two 32 bit numbers.
-    assert seed >= 0 and seed < 1<<64
-    seed_lo = int((seed & 0xffffffff) - (1<<31))
-    seed_hi = int((seed >> 32) - (1<<31))
+    assert seed >= 0 and seed < 1 << 64
+    seed_lo = int((seed & 0xffffffff) - (1 << 31))
+    seed_hi = int((seed >> 32) - (1 << 31))
     numpy.random.seed((seed_lo, seed_hi))
 
     return seed
 
+
 def exec_argv(vars):
     """Execute command line arguments as python statements.
 
diff --git a/kwant/system.py b/kwant/system.py
index 0ae6e90c6fcff1d760a90bad4f09c9a608dda73e..53fc4aa37c78cb91e8a4d3f1e1e273efd66a3bdb 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -3,11 +3,14 @@
 from __future__ import division
 __all__ = ['System', 'FiniteSystem', 'InfiniteSystem']
 
-import abc, math, types
+import abc
+import math
+import types
 import numpy as np
 from kwant import physics
 import _system
 
+
 class System(object):
     """Abstract general low-level system.
 
@@ -56,7 +59,8 @@ class FiniteSystem(System):
     Instance Variables
     ------------------
     leads : sequence of lead objects
-        Each lead object has to provide at least a method ``self_energy(energy)``.
+        Each lead object has to provide at least a method
+        ``self_energy(energy)``.
     lead_neighbor_seqs : sequence of sequences of integers
         Each sub-sequence contains the indices of the system sites to which the
         lead is connected.
@@ -127,7 +131,6 @@ class InfiniteSystem(System):
 
     def inter_slice_hopping(self):
         """Hopping Hamiltonian between two slices of the infinite system."""
-        slice_size = self.slice_size
         slice_sites = xrange(self.slice_size)
         neighbor_sites = xrange(self.slice_size, self.graph.num_nodes)
         return self.hamiltonian_submatrix(slice_sites, neighbor_sites)[0]
@@ -172,7 +175,7 @@ class InfiniteSystem(System):
         hop = self.inter_slice_hopping()
         result.hop = np.empty(result.ham.shape, dtype=complex)
         result.hop[:, : hop.shape[1]] = hop
-        result.hop[:, hop.shape[1] :] = 0
+        result.hop[:, hop.shape[1]:] = 0
         return result
 
 
diff --git a/kwant/version.py b/kwant/version.py
index 92ff8d09e0967aab3fedd043a51ec9d3fcb1891b..7b5efeeddefa3787feb3ec8016aab7fcf1b1949b 100644
--- a/kwant/version.py
+++ b/kwant/version.py
@@ -1,7 +1,9 @@
-import subprocess, os
+import subprocess
+import os
 
 __all__ = ['version']
 
+
 # When changing this function, remember to also change its twin in ../setup.py.
 def get_version_from_git():
     kwant_dir = os.path.dirname(os.path.abspath(__file__))