diff --git a/kwant/_common.py b/kwant/_common.py
index 23e5084d53a5c765e6189cf964a1610b728b4bdf..d49177b833adbc7e64b1874f08ad53550970aaab 100644
--- a/kwant/_common.py
+++ b/kwant/_common.py
@@ -8,6 +8,8 @@
 
 import subprocess
 import os
+import numpy as np
+import numbers
 
 __all__ = ['version', 'KwantDeprecationWarning', 'UserCodeError']
 
@@ -110,3 +112,21 @@ def ensure_isinstance(obj, typ, msg=None):
     if msg is None:
         msg = "Expecting an instance of {}.".format(typ.__name__)
     raise TypeError(msg)
+
+def ensure_rng(rng=None):
+    """Turn rng into a random number generator instance
+
+    If rng is None, return the RandomState instance used by np.random.
+    If rng is an integer, return a new RandomState instance seeded with rng.
+    If rng is already a RandomState instance, return it.
+    Otherwise raise ValueError.
+    """
+    if rng is None:
+        return np.random.mtrand._rand
+    if isinstance(rng, numbers.Integral):
+        return np.random.RandomState(rng)
+    if all(hasattr(rng, attr) for attr in ('random_sample', 'randn',
+                                           'randint', 'choice')):
+        return rng
+    raise ValueError("Expecting a seed or an object that offers the "
+                     "numpy.random.RandomState interface.")
diff --git a/kwant/linalg/tests/test_lll.py b/kwant/linalg/tests/test_lll.py
index cb04809046d722b061bdd9435c35205f15e00b65..c4d96f3091b73351faf7ecb810c28dd13a625b27 100644
--- a/kwant/linalg/tests/test_lll.py
+++ b/kwant/linalg/tests/test_lll.py
@@ -9,24 +9,25 @@
 
 import numpy as np
 from kwant.linalg import lll
+from kwant._common import ensure_rng
 
 def test_lll():
-    np.random.seed(1)
+    rng = ensure_rng(1)
     for i in range(50):
-        x = np.random.randint(4) + 1
-        mat = np.random.randn(x, x + np.random.randint(2))
-        c = 1.34 + .5 * np.random.rand()
+        x = rng.randint(4) + 1
+        mat = rng.randn(x, x + rng.randint(2))
+        c = 1.34 + .5 * rng.random_sample()
         reduced_mat, coefs = lll.lll(mat)
         assert lll.is_c_reduced(reduced_mat, c)
         assert np.allclose(np.dot(mat.T, coefs), reduced_mat.T)
 
 
 def test_cvp():
-    np.random.seed(0)
+    rng = ensure_rng(0)
     for i in range(10):
-        mat = np.random.randn(4, 4)
+        mat = rng.randn(4, 4)
         mat = lll.lll(mat)[0]
         for j in range(4):
-            point = 50 * np.random.randn(4)
+            point = 50 * rng.randn(4)
             assert np.array_equal(lll.cvp(point, mat, 10)[:3],
                                   lll.cvp(point, mat, 3))
diff --git a/kwant/physics/tests/test_leads.py b/kwant/physics/tests/test_leads.py
index 9b8527e066804d3b50397e5d9d5fb6026f88ce1f..4f8d8f4ad7b56b7e88825ce5d7ee433f7f11ef30 100644
--- a/kwant/physics/tests/test_leads.py
+++ b/kwant/physics/tests/test_leads.py
@@ -11,6 +11,7 @@ import numpy as np
 from numpy.testing import assert_almost_equal
 import scipy.linalg as la
 from kwant.physics import leads
+from kwant._common import ensure_rng
 import kwant
 
 modes_se = leads.selfenergy
@@ -294,14 +295,14 @@ def check_equivalence(h, t, n, sym='', particle_hole=None, chiral=None, time_rev
 def test_symm_algorithm_equivalence():
     """Test different stabilization methods in the computation of modes,
     in the presence and/or absence of the discrete symmetries."""
-    np.random.seed(400)
+    rng = ensure_rng(400)
     for n in (12, 20, 40):
         for sym in kwant.rmt.sym_list:
             # Random onsite and hopping matrices in symmetry class
-            h_cell = kwant.rmt.gaussian(n, sym)
+            h_cell = kwant.rmt.gaussian(n, sym, rng=rng)
             # Hopping is an offdiagonal block of a Hamiltonian. We rescale it
             # to ensure that there are modes at the Fermi level.
-            h_hop = 10 * kwant.rmt.gaussian(2*n, sym)[:n, n:]
+            h_hop = 10 * kwant.rmt.gaussian(2*n, sym, rng=rng)[:n, n:]
 
             if kwant.rmt.p(sym):
                 p_mat = np.array(kwant.rmt.h_p_matrix[sym])
@@ -430,10 +431,10 @@ def test_PHS_TRIM_degenerate_ordering():
     sy = np.array([[0,-1j],[1j,0]])
     sz = np.array([[1,0],[0,-1]])
     ### P squares to 1 ###
-    np.random.seed(42)
+    rng = ensure_rng(42)
     dims = (4, 10, 20)
     ts = (1.0, 1.7, 13.8)
-    rand_hop = 1j*(0.1+np.random.rand())
+    rand_hop = 1j*(0.1+rng.random_sample())
     hop = la.block_diag(*[t*rand_hop*np.eye(dim) for t, dim in zip(ts, dims)])
 
     # Particle-hole operator
@@ -447,11 +448,11 @@ def test_PHS_TRIM_degenerate_ordering():
     ###########
 
     ### P squares to -1 ###
-    np.random.seed(1337)
+    ensure_rng(1337)
     dims = (1, 4, 16)
     ts = (1.0, 17.2, 13.4)
 
-    hop_mat = np.kron(sz, 1j*(0.1+np.random.rand())*np.eye(2))
+    hop_mat = np.kron(sz, 1j*(0.1+rng.random_sample())*np.eye(2))
     blocks = []
     for t, dim in zip(ts, dims):
         blocks += dim*[t*hop_mat]
@@ -486,14 +487,14 @@ def test_PHS_TRIM_degenerate_ordering():
 
 
 def test_modes_symmetries():
-    np.random.seed(10)
+    rng = ensure_rng(10)
     for n in (4, 8, 40, 60):
         for sym in kwant.rmt.sym_list:
             # Random onsite and hopping matrices in symmetry class
-            h_cell = kwant.rmt.gaussian(n, sym)
+            h_cell = kwant.rmt.gaussian(n, sym, rng=rng)
             # Hopping is an offdiagonal block of a Hamiltonian. We rescale it
             # to ensure that there are modes at the Fermi level.
-            h_hop = 10 * kwant.rmt.gaussian(2*n, sym)[:n, n:]
+            h_hop = 10 * kwant.rmt.gaussian(2*n, sym, rng=rng)[:n, n:]
 
             if kwant.rmt.p(sym):
                 p_mat = np.array(kwant.rmt.h_p_matrix[sym])
@@ -557,7 +558,7 @@ def test_modes_symmetries():
 
 def test_PHS_TRIM():
     """Test the function that makes particle-hole symmetric modes at a TRIM. """
-    np.random.seed(10)
+    rng = ensure_rng(10)
     for n in (4, 8, 16, 40, 60):
         for sym in kwant.rmt.sym_list:
             if kwant.rmt.p(sym):
@@ -569,13 +570,13 @@ def test_PHS_TRIM():
                     for nmodes in (1, 3, n//4, n//2, n):
                         # Random matrix of 'modes.' Take part of a unitary matrix to
                         # ensure that the modes form a basis.
-                        modes = np.random.rand(n, n) + 1j*np.random.rand(n, n)
+                        modes = rng.random_sample((n, n)) + 1j*rng.random_sample((n, n))
                         modes = la.expm(1j*(modes + modes.T.conj()))[:n, :nmodes]
                         # Ensure modes are particle-hole symmetric and normalized
                         modes = modes + p_mat.dot(modes.conj())
                         modes = np.array([col/np.linalg.norm(col) for col in modes.T]).T
                         # Mix the modes with a random unitary transformation
-                        U = np.random.rand(nmodes, nmodes) + 1j*np.random.rand(nmodes, nmodes)
+                        U = rng.random_sample((nmodes, nmodes)) + 1j*rng.random_sample((nmodes, nmodes))
                         U = la.expm(1j*(U + U.T.conj()))
                         modes = modes.dot(U)
                         # Make the modes PHS symmetric using the method for a TRIM.
@@ -589,13 +590,13 @@ def test_PHS_TRIM():
                     for nmodes in (2, 4, n//2, n):
                         # Random matrix of 'modes.' Take part of a unitary matrix to
                         # ensure that the modes form a basis.
-                        modes = np.random.rand(n, n) + 1j*np.random.rand(n, n)
+                        modes = rng.random_sample((n, n)) + 1j*rng.random_sample((n, n))
                         modes = la.expm(1j*(modes + modes.T.conj()))[:n, :nmodes]
                         # Ensure modes are particle-hole symmetric and orthonormal.
                         modes[:, nmodes//2:] = p_mat.dot(modes[:, :nmodes//2].conj())
                         modes = la.qr(modes, mode='economic')[0]
                         # Mix the modes with a random unitary transformation
-                        U = np.random.rand(nmodes, nmodes) + 1j*np.random.rand(nmodes, nmodes)
+                        U = rng.random_sample((nmodes, nmodes)) + 1j*rng.random_sample((nmodes, nmodes))
                         U = la.expm(1j*(U + U.T.conj()))
                         modes = modes.dot(U)
                         # Make the modes PHS symmetric using the method for a TRIM.
diff --git a/kwant/physics/tests/test_noise.py b/kwant/physics/tests/test_noise.py
index 23dbbc1523b42352ee2e92fd274af5d5af3749b1..144a298f7173d3ea13b5a00b00e8868353cc1a10 100644
--- a/kwant/physics/tests/test_noise.py
+++ b/kwant/physics/tests/test_noise.py
@@ -11,18 +11,19 @@ from pytest import raises
 from numpy.testing import assert_almost_equal
 import kwant
 from kwant.physics import two_terminal_shotnoise
+from kwant._common import ensure_rng
 
 n = 5
 chain = kwant.lattice.chain()
 
 def twoterminal_system():
-    np.random.seed(11)
+    rng = ensure_rng(11)
     system = kwant.Builder()
     lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
-    h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+    h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
     h += h.conjugate().transpose()
     h *= 0.8
-    t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+    t = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
     lead[chain(0)] = h
     system[chain(0)] = h * 1.2
     lead[chain(0), chain(1)] = t
diff --git a/kwant/plotter.py b/kwant/plotter.py
index 7213caa555feac7f1572655c9305d4e4a6b506f2..ae1e70e499f5f5c80169a52146035f499b8d7349 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -16,7 +16,6 @@ system in two or three dimensions.
 
 from collections import defaultdict
 import warnings
-import random
 import numpy as np
 import tinyarray as ta
 from scipy import spatial, interpolate
@@ -43,7 +42,7 @@ except ImportError:
                   "functions will work.", RuntimeWarning)
     mpl_enabled = False
 
-from . import system, builder, physics
+from . import system, builder, physics, _common
 
 __all__ = ['plot', 'map', 'bands', 'sys_leads_sites', 'sys_leads_hoppings',
            'sys_leads_pos', 'sys_leads_hopping_pos', 'mask_interpolate']
@@ -80,9 +79,10 @@ def nparray_if_array(var):
     return np.asarray(var) if isarray(var) else var
 
 
-def _sample_array(array, n_samples):
+def _sample_array(array, n_samples, rng=None):
+    rng = _common.ensure_rng(rng)
     la = len(array)
-    return array[random.sample(range(la), min(n_samples, la))]
+    return array[rng.choice(range(la), min(n_samples, la))]
 
 
 if mpl_enabled:
diff --git a/kwant/rmt.py b/kwant/rmt.py
index 6442fc0a90e220277086f35a9a07576f327e4299..58c06b573adc4196f959d603d038c8f98a717235 100644
--- a/kwant/rmt.py
+++ b/kwant/rmt.py
@@ -9,9 +9,10 @@
 __all__ = ['gaussian', 'circular']
 
 import numpy as np
+from ._common import ensure_rng
 
 qr = np.linalg.qr
-randn = np.random.randn
+
 sym_list = 'A', 'AI', 'AII', 'AIII', 'BDI', 'CII', 'D', 'DIII', 'C', 'CI'
 
 
@@ -58,7 +59,7 @@ h_p_matrix = {'C': [[0, 1j], [-1j, 0]], 'CI': [[0, 1j], [-1j, 0]],
               'D': [[1]], 'DIII': [[0, 1], [1, 0]], 'BDI': [[1]]}
 
 
-def gaussian(n, sym='A', v=1.):
+def gaussian(n, sym='A', v=1., rng=None):
     """Make a n * n random Gaussian Hamiltonian.
 
     Parameters
@@ -71,6 +72,9 @@ def gaussian(n, sym='A', v=1.):
     v : float
         Variance every degree of freedom of the Hamiltonian. The probaility
         distribution of the Hamiltonian is `P(H) = exp(-Tr(H^2) / 2 v^2)`.
+    rng: int or rng (optional)
+        Seed or random number generator. If no 'rng' is provided the random
+        number generator from numpy will be used.
 
     Returns
     -------
@@ -117,6 +121,10 @@ def gaussian(n, sym='A', v=1.):
                          'symmetry class CII.')
     factor = v / np.sqrt(2)
 
+    # define random number generator
+    rng = ensure_rng(rng)
+    randn = rng.randn
+
     # Generate a Gaussian matrix of appropriate dtype.
     if sym == 'AI':
         h = randn(n, n)
@@ -152,7 +160,7 @@ def gaussian(n, sym='A', v=1.):
     return h
 
 
-def circular(n, sym='A', charge=None):
+def circular(n, sym='A', charge=None, rng=None):
     """Make a n * n matrix belonging to a symmetric circular ensemble.
 
     Parameters
@@ -167,6 +175,9 @@ def circular(n, sym='A', charge=None):
         classes D and DIII, should be from 0 to n in classes AIII and BDI,
         and should be from 0 to n / 2 in class CII. If charge is None,
         it is drawn from a binomial distribution with p = 1 / 2.
+    rng: int or rng (optional)
+        Seed or random number generator. If no 'rng' is passed, the random
+        number generator provided by numpy will be used.
 
     Returns
     -------
@@ -198,6 +209,8 @@ def circular(n, sym='A', charge=None):
     as detailed in arXiv:math-ph/0609050. For a reason as yet unknown, scipy
     implementation of QR decomposition also works for symplectic matrices.
     """
+    rng = ensure_rng(rng)
+    randn = rng.randn
 
     if sym not in sym_list:
         raise ValueError('Unknown symmetry type.')
@@ -257,14 +270,14 @@ def circular(n, sym='A', charge=None):
     elif sym in ('AIII', 'BDI', 'CII'):
         if sym != 'CII':
             if charge is None:
-                diag = 2 * np.random.randint(2, size=(n,)) - 1
+                diag = 2 * rng.randint(2, size=(n,)) - 1
             elif (0 <= charge <= n) and int(charge) == charge:
                 diag = np.array(charge * [-1] + (n - charge) * [1])
             else:
                 raise ValueError('Impossible value of topological invariant.')
         else:
             if charge is None:
-                diag = 2 * np.random.randint(2, size=(n // 2,)) - 1
+                diag = 2 * rng.randint(2, size=(n // 2,)) - 1
                 diag = np.resize(diag, (2, n // 2)).T.flatten()
             elif (0 <= charge <= n // 2) and int(charge) == charge:
                 charge *= 2
diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py
index 99e2bd729efb32ec6656c641d0350d4dbe7e89a2..1d5fd44c7d3681214b7c214399a372fd3b1f134a 100644
--- a/kwant/solvers/tests/_test_sparse.py
+++ b/kwant/solvers/tests/_test_sparse.py
@@ -12,6 +12,7 @@ import numpy as np
 from pytest import raises
 from numpy.testing import assert_almost_equal
 import kwant
+from kwant._common import ensure_rng
 
 n = 5
 chain = kwant.lattice.chain()
@@ -39,19 +40,20 @@ def assert_modes_equal(modes1, modes2):
 # and that solving for a subblock of a scattering matrix is the same as taking
 # a subblock of the full scattering matrix.
 def test_output(smatrix):
-    np.random.seed(3)
+    rng = ensure_rng(3)
     system = kwant.Builder()
     left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
     right_lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
     for b, site in [(system, chain(0)), (system, chain(1)),
                     (left_lead, chain(0)), (right_lead, chain(0))]:
-        h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+        h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
         h += h.conjugate().transpose()
         b[site] = h
     for b, hopp in [(system, (chain(0), chain(1))),
                     (left_lead, (chain(0), chain(1))),
                     (right_lead, (chain(0), chain(1)))]:
-        b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+        b[hopp] = (10 * rng.random_sample((n, n)) +
+                   1j * rng.random_sample((n, n)))
     system.attach_lead(left_lead)
     system.attach_lead(right_lead)
     fsyst = system.finalized()
@@ -81,18 +83,19 @@ def test_output(smatrix):
 
 # Test that a system with one lead has unitary scattering matrix.
 def test_one_lead(smatrix):
-    np.random.seed(3)
+    rng = ensure_rng(3)
     system = kwant.Builder()
     lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
     for b, site in [(system, chain(0)), (system, chain(1)),
                     (system, chain(2)), (lead, chain(0))]:
-        h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+        h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
         h += h.conjugate().transpose()
         b[site] = h
     for b, hopp in [(system, (chain(0), chain(1))),
                     (system, (chain(1), chain(2))),
                     (lead, (chain(0), chain(1)))]:
-        b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+        b[hopp] = (10 * rng.random_sample((n, n)) +
+                   1j * rng.random_sample((n, n)))
     system.attach_lead(lead)
     fsyst = system.finalized()
 
@@ -161,13 +164,13 @@ def test_two_equal_leads(smatrix):
         t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
         assert_almost_equal(t_elements, t_el_should_be)
         assert_almost_equal(sol.transmission(1,0), n_modes)
-    np.random.seed(11)
+    rng = ensure_rng(11)
     system = kwant.Builder()
     lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
-    h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+    h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
     h += h.conjugate().transpose()
     h *= 0.8
-    t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+    t = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
     lead[chain(0)] = system[chain(0)] = h
     lead[chain(0), chain(1)] = t
     system.attach_lead(lead)
@@ -191,14 +194,14 @@ def test_two_equal_leads(smatrix):
 
 # Test a more complicated graph with non-singular hopping.
 def test_graph_system(smatrix):
-    np.random.seed(11)
+    rng = ensure_rng(11)
     system = kwant.Builder()
     lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
-    h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+    h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
     h += h.conjugate().transpose()
     h *= 0.8
-    t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
-    t1 = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+    t = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
+    t1 = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
     lead[sq(0, 0)] = system[[sq(0, 0), sq(1, 0)]] = h
     lead[sq(0, 1)] = system[[sq(0, 1), sq(1, 1)]] = 4 * h
     for builder in [system, lead]:
@@ -224,15 +227,15 @@ def test_graph_system(smatrix):
 
 # Test a system with singular hopping.
 def test_singular_graph_system(smatrix):
-    np.random.seed(11)
+    rng = ensure_rng(11)
 
     system = kwant.Builder()
     lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
-    h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+    h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
     h += h.conjugate().transpose()
     h *= 0.8
-    t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
-    t1 = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+    t = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
+    t1 = 4 * rng.random_sample((n, n)) + 4j * rng.random_sample((n, n))
     lead[sq(0, 0)] = system[[sq(0, 0), sq(1, 0)]] = h
     lead[sq(0, 1)] = system[[sq(0, 1), sq(1, 1)]] = 4 * h
     for builder in [system, lead]:
@@ -345,19 +348,20 @@ def test_many_leads(*factories):
 # Test equivalence between self-energy and scattering matrix representations.
 # Also check that transmission works.
 def test_selfenergy(greens_function, smatrix):
-    np.random.seed(4)
+    rng = ensure_rng(4)
     system = kwant.Builder()
     left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
     right_lead = kwant.Builder(kwant.TranslationalSymmetry((1,)))
     for b, site in [(system, chain(0)), (system, chain(1)),
                  (left_lead, chain(0)), (right_lead, chain(0))]:
-        h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+        h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
         h += h.conjugate().transpose()
         b[site] = h
     for b, hopp in [(system, (chain(0), chain(1))),
                     (left_lead, (chain(0), chain(1))),
                     (right_lead, (chain(0), chain(1)))]:
-        b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+        b[hopp] = (10 * rng.random_sample((n, n)) +
+                   1j * rng.random_sample((n, n)))
     system.attach_lead(left_lead)
     system.attach_lead(right_lead)
     fsyst = system.finalized()
@@ -390,17 +394,18 @@ def test_selfenergy(greens_function, smatrix):
 
 
 def test_selfenergy_reflection(greens_function, smatrix):
-    np.random.seed(4)
+    rng = ensure_rng(4)
     system = kwant.Builder()
     left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
     for b, site in [(system, chain(0)), (system, chain(1)),
                  (left_lead, chain(0))]:
-        h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+        h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
         h += h.conjugate().transpose()
         b[site] = h
     for b, hopp in [(system, (chain(0), chain(1))),
                     (left_lead, (chain(0), chain(1)))]:
-        b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+        b[hopp] = (10 * rng.random_sample((n, n)) +
+                   1j * rng.random_sample((n, n)))
     system.attach_lead(left_lead)
     fsyst = system.finalized()
 
@@ -457,7 +462,7 @@ def test_wavefunc_ldos_consistency(wave_function, ldos):
     L = 2
     W = 3
 
-    np.random.seed(31)
+    rng = ensure_rng(31)
     syst = kwant.Builder()
     left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
     top_lead = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
@@ -466,12 +471,13 @@ def test_wavefunc_ldos_consistency(wave_function, ldos):
                      (left_lead, [square(0, y) for y in range(W)]),
                      (top_lead, [square(x, 0) for x in range(L)])]:
         for site in sites:
-            h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+            h = rng.random_sample((n, n)) + 1j * rng.random_sample((n, n))
             h += h.conjugate().transpose()
             b[site] = h
         for hopping_kind in square.neighbors():
             for hop in hopping_kind(b):
-                b[hop] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+                b[hop] = (10 * rng.random_sample((n, n)) +
+                          1j * rng.random_sample((n, n)))
     syst.attach_lead(left_lead)
     syst.attach_lead(top_lead)
     syst = syst.finalized()
diff --git a/kwant/tests/test_lattice.py b/kwant/tests/test_lattice.py
index 45337c1074eb0b929789ce603f1f17a4af5290aa..b15ea16ddba997ec93ad74acb5d672457a74e06d 100644
--- a/kwant/tests/test_lattice.py
+++ b/kwant/tests/test_lattice.py
@@ -11,18 +11,19 @@ import numpy as np
 import tinyarray as ta
 from pytest import raises
 from kwant import lattice, builder
+from kwant._common import ensure_rng
 
 
 def test_closest():
-    np.random.seed(4)
+    rng = ensure_rng(4)
     lat = lattice.general(((1, 0), (0.5, sqrt(3)/2)))
     for i in range(50):
-        point = 20 * np.random.rand(2)
+        point = 20 * rng.random_sample(2)
         closest = lat(*lat.closest(point)).pos
         assert np.linalg.norm(point - closest) <= 1 / sqrt(3)
-    lat = lattice.general(np.random.randn(3, 3))
+    lat = lattice.general(rng.randn(3, 3))
     for i in range(50):
-        tag = np.random.randint(10, size=(3,))
+        tag = rng.randint(10, size=(3,))
         assert lat.closest(lat(*tag).pos) == tag
 
 
@@ -82,11 +83,11 @@ def test_shape():
 
 
 def test_wire():
-    np.random.seed(5)
-    vecs = np.random.randn(3, 3)
+    rng = ensure_rng(5)
+    vecs = rng.randn(3, 3)
     vecs[0] = [1, 0, 0]
-    center = np.random.randn(3)
-    lat = lattice.general(vecs, np.random.randn(4, 3))
+    center = rng.randn(3)
+    lat = lattice.general(vecs, rng.randn(4, 3))
     syst = builder.Builder(lattice.TranslationalSymmetry((2, 0, 0)))
     def wire_shape(pos):
         pos = np.array(pos)
@@ -156,13 +157,13 @@ def test_translational_symmetry():
     # Test add_site_family on random lattices and symmetries by ensuring that
     # it's possible to add site groups that are compatible with a randomly
     # generated symmetry with proper vectors.
-    np.random.seed(30)
-    vec = np.random.randn(3, 5)
+    rng = ensure_rng(30)
+    vec = rng.randn(3, 5)
     lat = lattice.general(vec)
     total = 0
     for k in range(1, 4):
         for i in range(10):
-            sym_vec = np.random.randint(-10, 10, size=(k, 3))
+            sym_vec = rng.randint(-10, 10, size=(k, 3))
             if np.linalg.matrix_rank(sym_vec) < k:
                 continue
             total += 1
@@ -173,12 +174,12 @@ def test_translational_symmetry():
 
 
 def test_translational_symmetry_reversed():
-    np.random.seed(30)
+    rng = ensure_rng(30)
     lat = lattice.general(np.identity(3))
     sites = [lat(i, j, k) for i in range(-2, 6) for j in range(-2, 6)
                           for k in range(-2, 6)]
     for i in range(4):
-            periods = np.random.randint(-5, 5, (3, 3))
+            periods = rng.randint(-5, 5, (3, 3))
             try:
                 sym = lattice.TranslationalSymmetry(*periods)
                 rsym = sym.reversed()
diff --git a/kwant/tests/test_rmt.py b/kwant/tests/test_rmt.py
index 3cad147fd15fc28078af1dc59e2fe683b37312ab..b7e3e64fdfacb5c57c1c2951ecfb7697f80df949 100644
--- a/kwant/tests/test_rmt.py
+++ b/kwant/tests/test_rmt.py
@@ -10,17 +10,18 @@ import numpy as np
 from scipy import stats
 from pytest import raises
 from kwant import rmt
+from kwant._common import ensure_rng
 assert_allclose = np.testing.assert_allclose
 
 
 def test_gaussian_symmetries():
-    np.random.seed(10)
+    rng = ensure_rng(10)
     for n in (5, 8, 100, 200):
         for sym in rmt.sym_list:
             if sym not in ('A', 'D', 'AI') and n % 2:
                 raises(ValueError, rmt.gaussian, 5, sym)
                 continue
-            h = rmt.gaussian(n, sym)
+            h = rmt.gaussian(n, sym, rng=rng)
             if rmt.t(sym):
                 t_mat = np.array(rmt.h_t_matrix[sym])
                 t_mat = np.kron(np.identity(n // len(t_mat)), t_mat)
@@ -38,23 +39,23 @@ def test_gaussian_symmetries():
 
 
 def test_gaussian_distributions():
-    np.random.seed(1)
+    rng = ensure_rng(1)
     n = 8
     for sym in rmt.sym_list:
-        matrices = np.array([rmt.gaussian(n, sym)[-1, 0] for i in range(3000)])
+        matrices = np.array([rmt.gaussian(n, sym, rng=rng)[-1, 0] for i in range(3000)])
         matrices = matrices.imag if sym in ('D', 'BDI') else matrices.real
         ks = stats.kstest(matrices, 'norm')
         assert (ks[1] > 0.1), (sym, ks)
 
 
 def test_circular():
-    np.random.seed(10)
+    rng = ensure_rng(10)
     n = 6
     sy = np.kron(np.identity(n // 2), [[0, 1j], [-1j, 0]])
     for sym in rmt.sym_list:
         if rmt.t(sym) == -1 or rmt.p(sym) == -1:
             raises(ValueError, rmt.circular, 5, sym)
-        s = rmt.circular(n, sym)
+        s = rmt.circular(n, sym, rng=rng)
         assert_allclose(np.dot(s, s.T.conj()), np.identity(n), atol=1e-9,
                         err_msg='Unitarity broken in ' + sym)
         if rmt.t(sym):
@@ -75,8 +76,8 @@ def test_circular():
     # Check for distribution properties if the ensemble is a symmetric group.
     f = lambda x: x[0] / np.linalg.norm(x)
     for sym in ('A', 'C'):
-        sample_distr = np.apply_along_axis(f, 0, np.random.randn(2 * n, 1000))
-        s_sample = np.array([rmt.circular(n, sym)
+        sample_distr = np.apply_along_axis(f, 0, rng.randn(2 * n, 1000))
+        s_sample = np.array([rmt.circular(n, sym, rng=rng)
                              for i in range(1000)])
         assert stats.ks_2samp(sample_distr, s_sample[:, 0, 0].real)[1] > 0.1, \
                 'Noncircular distribution in ' + sym
@@ -87,8 +88,8 @@ def test_circular():
         assert stats.ks_2samp(sample_distr, s_sample[:, 2, 3].imag)[1] > 0.1, \
                 'Noncircular distribution in ' + sym
 
-    sample_distr = np.apply_along_axis(f, 0, np.random.randn(n, 500))
-    s_sample = np.array([rmt.circular(n, 'D') for i in range(500)])
+    sample_distr = np.apply_along_axis(f, 0, rng.randn(n, 500))
+    s_sample = np.array([rmt.circular(n, 'D', rng=rng) for i in range(500)])
     ks = stats.ks_2samp(sample_distr, s_sample[:, 0, 0])
     assert ks[1] > 0.1, 'Noncircular distribution in D ' + str(ks)
     ks = stats.ks_2samp(sample_distr, s_sample[:, 3, 2])
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index f06da88a2b9f8296a866bec9ce772146effb2833..1f2ae57c52cc53e727facf3ee741e901ce44a99e 100644
--- a/kwant/tests/test_system.py
+++ b/kwant/tests/test_system.py
@@ -12,6 +12,7 @@ from pytest import raises
 import numpy as np
 from scipy import sparse
 import kwant
+from kwant._common import ensure_rng
 
 def test_hamiltonian_submatrix():
     syst = kwant.Builder()
@@ -59,10 +60,10 @@ def test_hamiltonian_submatrix():
     np.testing.assert_array_equal(mat_sp, mat_dense)
 
     # Test precalculation of modes.
-    np.random.seed(5)
+    rng = ensure_rng(5)
     lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
     lead[chain(0)] = np.zeros((2, 2))
-    lead[chain(0), chain(1)] = np.random.randn(2, 2)
+    lead[chain(0), chain(1)] = rng.randn(2, 2)
     syst.attach_lead(lead)
     syst2 = syst.finalized()
     smatrix = kwant.smatrix(syst2, .1).data