diff --git a/kwant/contrib/__init__.py b/kwant/contrib/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..c52331e99965033fa443a80c3ab4cec0427ab84c
--- /dev/null
+++ b/kwant/contrib/__init__.py
@@ -0,0 +1,10 @@
+# Copyright 2011-2013 kwant authors.
+#
+# This file is part of kwant. It is subject to the license terms in the
+# LICENSE file found in the top-level directory of this distribution and at
+# http://kwant-project.org/license. A list of kwant authors can be found in
+# the AUTHORS file at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+__all__ = ['rmt']
+from . import rmt
diff --git a/kwant/contrib/rmt.py b/kwant/contrib/rmt.py
new file mode 100644
index 0000000000000000000000000000000000000000..67b8b9c83f4fe1bde1f8478105367f2e7bdbc865
--- /dev/null
+++ b/kwant/contrib/rmt.py
@@ -0,0 +1,277 @@
+# Copyright 2011-2013 kwant authors.
+#
+# This file is part of kwant. It is subject to the license terms in the
+# LICENSE file found in the top-level directory of this distribution and at
+# http://kwant-project.org/license. A list of kwant authors can be found in
+# the AUTHORS file 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
+
+qr = np.linalg.qr
+randn = np.random.randn
+sym_list = 'A', 'AI', 'AII', 'AIII', 'BDI', 'CII', 'D', 'DIII', 'C', 'CI'
+
+
+def t(sym):
+ """Return the value of time-reversal symmetry squared (1, 0, or -1)"""
+ if sym not in sym_list:
+ raise ValueError('Non-existent symmetry class.')
+ if sym in ('CI', 'AI', 'BDI'):
+ return 1
+ elif sym in ('CII', 'AII', 'DIII'):
+ return -1
+ else:
+ return 0
+
+
+def p(sym):
+ """Return the value of particle-hole symmetry squared (1, 0, or -1)"""
+ if sym not in sym_list:
+ raise ValueError('Non-existent symmetry class.')
+ if sym in ('D', 'DIII', 'BDI'):
+ return 1
+ elif sym in ('C', 'CI', 'CII'):
+ return -1
+ else:
+ return 0
+
+
+def c(sym):
+ """Return 1 if the system has chiral symmetry, and 0 otherwise."""
+ if (t(sym) and p(sym)) or sym == 'AIII':
+ return 1
+ else:
+ return 0
+
+
+h_t_matrix = {'AI': [[1]], 'CI': [[0, 1], [1, 0]], 'BDI': [[1, 0], [0, -1]],
+ 'AII': [[0, 1j], [-1j, 0]],
+ 'CII': [[0, 0, 1j, 0], [0, 0, 0, 1j],
+ [-1j, 0, 0, 0], [0, -1j, 0, 0]],
+ 'DIII': [[0, 1j], [-1j, 0]]}
+h_p_matrix = {'C': [[0, 1j], [-1j, 0]], 'CI': [[0, 1j], [-1j, 0]],
+ 'CII': [[0, 0, 1j, 0], [0, 0, 0, -1j],
+ [-1j, 0, 0, 0], [0, 1j, 0, 0]],
+ 'D': [[1]], 'DIII': [[0, 1], [1, 0]], 'BDI': [[1]]}
+
+
+def gaussian(n, sym='A', v=1.):
+ """Make a n * n random Gaussian Hamiltonian.
+
+ Parameters
+ ----------
+ n : int
+ Size of the Hamiltonian. It should be even for all the classes
+ except A, D, and AI, and in class CII it should be a multiple of 4.
+ sym : one of 'A', 'AI', 'AII', 'AIII', 'BDI', 'CII', 'D', 'DIII', 'C', 'CI'
+ Altland-Zirnbauer symmetry class of the Hamiltonian.
+ 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)`.
+
+ Returns
+ -------
+ h : numpy.ndarray
+ A numpy array drawn from a corresponding Gaussian ensemble.
+
+ Notes
+ -----
+ The representations of symmetry operators are chosen according to
+ Phys. Rev. B 85, 165409.
+
+ Matrix indices are grouped first according to orbital number,
+ then sigma-index, then tau-index.
+
+ Chiral (sublattice) symmetry C always reads:
+ H = -tau_z H tau_z.
+ Time reversal symmetry T reads:
+ AI: H = H^*.
+ BDI: H = tau_z H^* tau_z.
+ CI: H = tau_x H^* tau_x.
+ AII, CII: H = sigma_y H^* sigma_y.
+ DIII: H = tau_y H^* tau_y.
+ Particle-hole symmetry reads:
+ C, CI: H = -tau_y H^* tau_y.
+ CII: H = -tau_z sigma_y H^* tau_z sigma_y.
+ D, BDI: H = -H^*.
+ DIII: H = -tau_x H^* tau_x.
+
+ This implementation should be sufficiently efficient for large matrices,
+ since it avoids any matrix multiplication.
+ """
+ if sym not in sym_list:
+ raise ValueError('Unknown symmetry type.')
+ if (c(sym) or t(sym) == -1 or p(sym) == -1) and n % 2:
+ raise ValueError('Matrix dimension should be even in chosen symmetry'
+ 'class.')
+ else:
+ tau_z = np.array((n // 2) * [1, -1])
+ idx_x = np.arange(n) + tau_z
+ if sym == 'CII' and n % 4:
+ raise ValueError('Matrix dimension should be a multiple of 4 in'
+ 'symmetry class CII.')
+ factor = v / np.sqrt(2)
+
+ # Generate a Gaussian matrix of appropriate dtype.
+ if sym == 'AI':
+ h = randn(n, n)
+ elif sym in ('D', 'BDI'):
+ h = 1j * randn(n, n)
+ else:
+ h = randn(n, n) + 1j * randn(n, n)
+
+ # Ensure Hermiticity.
+ h += h.T.conj()
+
+ # Ensure Chiral symmetry.
+ if c(sym):
+ h -= tau_z.reshape(-1, 1) * h * tau_z
+ factor *= 0.5
+
+ # Ensure the necessary anti-unitary symmetry.
+ if sym in ('AII', 'DIII'):
+ h += tau_z.reshape(-1, 1) * h[idx_x][:, idx_x].conj() * tau_z
+ factor /= np.sqrt(2)
+ elif sym in ('C', 'CI'):
+ h -= tau_z.reshape(-1, 1) * h[idx_x][:, idx_x].conj() * tau_z
+ factor /= np.sqrt(2)
+ elif sym == 'CII':
+ sigma_z = np.array((n // 4) * [1, 1, -1, -1])
+ idx_sigma_x = np.arange(n) + 2 * sigma_z
+ h += sigma_z.reshape(-1, 1) * h[idx_sigma_x][:, idx_sigma_x].conj() * \
+ sigma_z
+ factor /= np.sqrt(2)
+
+ h *= factor
+
+ return h
+
+
+def circular(n, sym='A', charge=None):
+ """Make a n * n matrix belonging to a symmetric circular ensemble.
+
+ Parameters
+ ----------
+ n : int
+ Size of the matrix. It should be even for the classes C, CI, CII,
+ AII, DIII (either T^2 = -1 or P^2 = -1).
+ sym : one of 'A', 'AI', 'AII', 'AIII', 'BDI', 'CII', 'D', 'DIII', 'C', 'CI'
+ Altland-Zirnbauer symmetry class of the matrix.
+ charge : int or None
+ Topological invariant of the matrix. Should be one of 1, -1 in symmetry
+ 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.
+
+ Returns
+ -------
+ s : numpy.ndarray
+ A numpy array drawn from a corresponding circular ensemble.
+
+ Notes
+ -----
+ The representations of symmetry operators are chosen according to
+ Phys. Rev. B 85, 165409, except class D.
+
+ Matrix indices are grouped first according to channel number,
+ then sigma-index.
+
+ Chiral (sublattice) symmetry C always reads:
+ s = s^+.
+ Time reversal symmetry T reads:
+ AI, BDI: r = r^T.
+ CI: r = -sigma_y r^T sigma_y.
+ AII, DIII: r = -r^T.
+ CII: r = sigma_y r^T sigma_y.
+ Particle-hole symmetry reads:
+ CI: r = -sigma_y r^* sigma_y
+ C, CII: r = sigma_y r^* sigma_y
+ D, BDI: r = r^*.
+ DIII: -r = r^*.
+
+ This function uses QR decomposition to probe symmetric compact groups,
+ as detailed in arXiv:math-ph/0609050. For a reason as yet unknown, scipy
+ implementation of QR decomposition also works for symplectic matrices.
+ """
+
+ if sym not in sym_list:
+ raise ValueError('Unknown symmetry type.')
+ if (t(sym) == -1 or p(sym) == -1) and n % 2:
+ raise ValueError('n must be even in chosen symmetry class.')
+
+ # Prepare a real, complex or symplectic Gaussian matrix
+ # Real case.
+ if p(sym) == 1:
+ h = randn(n, n)
+ # Complex case.
+ elif p(sym) == 0:
+ h = randn(n, n) + 1j * randn(n, n)
+ # Symplectic case.
+ else: # p(sym) == -1
+ h = randn(n, n) + 1j * randn(n, n)
+ tau_z = np.array((n // 2) * [1, -1])
+ idx_x = np.arange(n) + tau_z
+ h -= tau_z.reshape(-1, 1) * h[idx_x][:, idx_x].conj() * tau_z
+ h *= 1j
+
+ # Generate a random matrix with proper electron-hole symmetry.
+ # This matrix is a random element of groups O(N), U(N), or USp(N).
+ s, h = qr(h)
+ if p(sym) == -1 and not np.allclose(np.diag(h, 1)[::2], 0, atol=1e-8):
+ raise RuntimeError('QR decomposition symmetry failure.')
+ h = np.diag(h)
+ h /= np.abs(h)
+ s *= h
+
+ # Ensure proper topological invariant in classes D and DIII.
+ if sym in ('D', 'DIII') and charge is not None:
+ if charge not in (-1, 1):
+ raise ValueError('Impossible value of topological invariant.')
+ det = np.linalg.det(s)
+ if sym == 'DIII':
+ det *= (-1) ** (n // 2)
+ if (charge > 0) != (det > 0):
+ idx = np.arange(n)
+ idx[-1] -= 1
+ idx[-2] += 1
+ s = s[idx]
+
+ # Add the proper time-reversal symmetry:
+ if sym in ('AI', 'CI'):
+ s = np.dot(s.T, s)
+ if sym == 'CI':
+ tau_z = np.array((n // 2) * [1, -1])
+ idx_x = np.arange(n) + tau_z
+ s = 1j * tau_z * s[:, idx_x]
+ elif sym == 'AII' or sym == 'DIII':
+ tau_z = np.array((n // 2) * [1, -1])
+ idx_x = np.arange(n) + tau_z
+ s = 1j * np.dot(s.T * tau_z, s[idx_x])
+
+ # Add the chiral symmetry:
+ elif sym in ('AIII', 'BDI', 'CII'):
+ if sym != 'CII':
+ if charge is None:
+ diag = 2 * np.random.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 = np.resize(diag, (2, n // 2)).T.flatten()
+ elif (0 <= charge <= n // 2) and int(charge) == charge:
+ charge *= 2
+ diag = np.array(charge * [-1] + (n - charge) * [1])
+ else:
+ raise ValueError('Impossible value of topological invariant.')
+
+ s = np.dot(diag * s.T.conj(), s)
+
+ return s
diff --git a/kwant/contrib/tests/test_rmt.py b/kwant/contrib/tests/test_rmt.py
new file mode 100644
index 0000000000000000000000000000000000000000..5bd749357159ff1747cfc5adb072684cd71bbf39
--- /dev/null
+++ b/kwant/contrib/tests/test_rmt.py
@@ -0,0 +1,90 @@
+# Copyright 2011-2013 kwant authors.
+#
+# This file is part of kwant. It is subject to the license terms in the
+# LICENSE file found in the top-level directory of this distribution and at
+# http://kwant-project.org/license. A list of kwant authors can be found in
+# the AUTHORS file at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+import numpy as np
+from scipy import stats
+from nose.tools import assert_raises
+from kwant.contrib import rmt
+assert_allclose = np.testing.assert_allclose
+
+
+def test_gaussian():
+ np.random.seed(10)
+ n = 8
+ for sym in rmt.sym_list:
+ if sym not in ('A', 'D', 'AI'):
+ assert_raises(ValueError, rmt.gaussian, 5, sym)
+ h = rmt.gaussian(n, sym)
+ 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)
+ assert_allclose(h, np.dot(t_mat, np.dot(h.conj(), t_mat)),
+ err_msg='TRS broken in ' + sym)
+ if rmt.p(sym):
+ p_mat = np.array(rmt.h_p_matrix[sym])
+ p_mat = np.kron(np.identity(n / len(p_mat)), p_mat)
+ assert_allclose(h, -np.dot(p_mat, np.dot(h.conj(), p_mat)),
+ err_msg='PHS broken in ' + sym)
+ if rmt.c(sym):
+ sz = np.kron(np.identity(n / 2), np.diag([1, -1]))
+ assert_allclose(h, -np.dot(sz, np.dot(h, sz)),
+ err_msg='SLS broken in ' + sym)
+
+ for sym in rmt.sym_list:
+ matrices = np.array([rmt.gaussian(n, sym)[-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)
+ 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:
+ assert_raises(ValueError, rmt.circular, 5, sym)
+ s = rmt.circular(n, sym)
+ assert_allclose(np.dot(s, s.T.conj()), np.identity(n), atol=1e-9,
+ err_msg='Unitarity broken in ' + sym)
+ if rmt.t(sym):
+ s1 = np.copy(s.T if rmt.p(sym) != -1
+ else np.dot(sy, np.dot(s.T, sy)))
+ s1 *= rmt.t(sym) * (-1 if rmt.p(sym) == -1 else 1)
+ assert_allclose(s, s1, atol=1e-9, err_msg='TRS broken in ' + sym)
+ if rmt.p(sym):
+ s1 = np.copy(s.conj() if rmt.p(sym) != -1
+ else np.dot(sy, np.dot(s.conj(), sy)))
+ if sym in ('DIII', 'CI'):
+ s1 *= -1
+ assert_allclose(s, s1, atol=1e-9, err_msg='PHS broken in ' + sym)
+ if rmt.c(sym):
+ assert_allclose(s, s.T.conj(), atol=1e-9,
+ err_msg='SLS broken in ' + sym)
+
+ # 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)
+ for i in range(1000)])
+ assert stats.ks_2samp(sample_distr, s_sample[:, 0, 0].real)[1] > 0.1, \
+ 'Noncircular distribution in ' + sym
+ assert stats.ks_2samp(sample_distr, s_sample[:, 3, 2].real)[1] > 0.1, \
+ 'Noncircular distribution in ' + sym
+ assert stats.ks_2samp(sample_distr, s_sample[:, 1, 1].imag)[1] > 0.1, \
+ 'Noncircular distribution in ' + sym
+ 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)])
+ 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])
+ assert ks[1] > 0.1, 'Noncircular distribution in D ' + str(ks)