Coverage for kwant/rmt.py : 78%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Copyright 2011-2013 Kwant authors.
2#
3# This file is part of Kwant. It is subject to the license terms in the file
4# LICENSE.rst found in the top-level directory of this distribution and at
5# https://kwant-project.org/license. A list of Kwant authors can be found in
6# the file AUTHORS.rst at the top-level directory of this distribution and at
7# https://kwant-project.org/authors.
9__all__ = ['gaussian', 'circular']
11import numpy as np
12from ._common import ensure_rng
14qr = np.linalg.qr
16sym_list = 'A', 'AI', 'AII', 'AIII', 'BDI', 'CII', 'D', 'DIII', 'C', 'CI'
19def t(sym):
20 """Return the value of time-reversal symmetry squared (1, 0, or -1)"""
21 if sym not in sym_list: 21 ↛ 22line 21 didn't jump to line 22, because the condition on line 21 was never true
22 raise ValueError('Non-existent symmetry class.')
23 if sym in ('CI', 'AI', 'BDI'):
24 return 1
25 elif sym in ('CII', 'AII', 'DIII'):
26 return -1
27 else:
28 return 0
31def p(sym):
32 """Return the value of particle-hole symmetry squared (1, 0, or -1)"""
33 if sym not in sym_list: 33 ↛ 34line 33 didn't jump to line 34, because the condition on line 33 was never true
34 raise ValueError('Non-existent symmetry class.')
35 if sym in ('D', 'DIII', 'BDI'):
36 return 1
37 elif sym in ('C', 'CI', 'CII'):
38 return -1
39 else:
40 return 0
43def c(sym):
44 """Return 1 if the system has chiral symmetry, and 0 otherwise."""
45 if (t(sym) and p(sym)) or sym == 'AIII':
46 return 1
47 else:
48 return 0
51h_t_matrix = {'AI': [[1]], 'CI': [[0, 1], [1, 0]], 'BDI': [[1, 0], [0, -1]],
52 'AII': [[0, 1j], [-1j, 0]],
53 'CII': [[0, 0, 1j, 0], [0, 0, 0, 1j],
54 [-1j, 0, 0, 0], [0, -1j, 0, 0]],
55 'DIII': [[0, 1j], [-1j, 0]]}
56h_p_matrix = {'C': [[0, 1j], [-1j, 0]], 'CI': [[0, 1j], [-1j, 0]],
57 'CII': [[0, 0, 1j, 0], [0, 0, 0, -1j],
58 [-1j, 0, 0, 0], [0, 1j, 0, 0]],
59 'D': [[1]], 'DIII': [[0, 1], [1, 0]], 'BDI': [[1]]}
62def gaussian(n, sym='A', v=1., rng=None):
63 """Make a n * n random Gaussian Hamiltonian.
65 Parameters
66 ----------
67 n : int
68 Size of the Hamiltonian. It should be even for all the classes
69 except A, D, and AI, and in class CII it should be a multiple of 4.
70 sym : one of 'A', 'AI', 'AII', 'AIII', 'BDI', 'CII', 'D', 'DIII', 'C', 'CI'
71 Altland-Zirnbauer symmetry class of the Hamiltonian.
72 v : float
73 Variance every degree of freedom of the Hamiltonian. The probaility
74 distribution of the Hamiltonian is `P(H) = exp(-Tr(H^2) / 2 v^2)`.
75 rng: int or rng (optional)
76 Seed or random number generator. If no 'rng' is provided the random
77 number generator from numpy will be used.
79 Returns
80 -------
81 h : numpy.ndarray
82 A numpy array drawn from a corresponding Gaussian ensemble.
84 Notes
85 -----
86 The representations of symmetry operators are chosen according to
87 Phys. Rev. B 85, 165409.
89 Matrix indices are grouped first according to orbital number,
90 then sigma-index, then tau-index.
92 Chiral (sublattice) symmetry C always reads:
93 H = -tau_z H tau_z.
94 Time reversal symmetry T reads:
95 AI: H = H^*.
96 BDI: H = tau_z H^* tau_z.
97 CI: H = tau_x H^* tau_x.
98 AII, CII: H = sigma_y H^* sigma_y.
99 DIII: H = tau_y H^* tau_y.
100 Particle-hole symmetry reads:
101 C, CI: H = -tau_y H^* tau_y.
102 CII: H = -tau_z sigma_y H^* tau_z sigma_y.
103 D, BDI: H = -H^*.
104 DIII: H = -tau_x H^* tau_x.
106 This implementation should be sufficiently efficient for large matrices,
107 since it avoids any matrix multiplication.
108 """
109 if sym not in sym_list: 109 ↛ 110line 109 didn't jump to line 110, because the condition on line 109 was never true
110 raise ValueError('Unknown symmetry type.')
111 if (c(sym) or t(sym) == -1 or p(sym) == -1):
112 if n % 2:
113 raise ValueError('Matrix dimension should be even in'
114 ' chosen symmetry class.')
115 else:
116 tau_z = np.array((n // 2) * [1, -1])
117 idx_x = np.arange(n) + tau_z
119 if sym == 'CII' and n % 4: 119 ↛ 120line 119 didn't jump to line 120, because the condition on line 119 was never true
120 raise ValueError('Matrix dimension should be a multiple of 4 in'
121 'symmetry class CII.')
122 factor = v / np.sqrt(2)
124 # define random number generator
125 rng = ensure_rng(rng)
126 randn = rng.randn
128 # Generate a Gaussian matrix of appropriate dtype.
129 if sym == 'AI':
130 h = randn(n, n)
131 elif sym in ('D', 'BDI'):
132 h = 1j * randn(n, n)
133 else:
134 h = randn(n, n) + 1j * randn(n, n)
136 # Ensure Hermiticity.
137 h += h.T.conj()
139 # Ensure Chiral symmetry.
140 if c(sym):
141 h -= tau_z.reshape(-1, 1) * h * tau_z
142 factor *= 0.5
144 # Ensure the necessary anti-unitary symmetry.
145 if sym in ('AII', 'DIII'):
146 h += tau_z.reshape(-1, 1) * h[idx_x][:, idx_x].conj() * tau_z
147 factor /= np.sqrt(2)
148 elif sym in ('C', 'CI'):
149 h -= tau_z.reshape(-1, 1) * h[idx_x][:, idx_x].conj() * tau_z
150 factor /= np.sqrt(2)
151 elif sym == 'CII':
152 sigma_z = np.array((n // 4) * [1, 1, -1, -1])
153 idx_sigma_x = np.arange(n) + 2 * sigma_z
154 h += (sigma_z.reshape(-1, 1) * h[idx_sigma_x][:, idx_sigma_x].conj() *
155 sigma_z)
156 factor /= np.sqrt(2)
158 h *= factor
160 return h
163def circular(n, sym='A', charge=None, rng=None):
164 """Make a n * n matrix belonging to a symmetric circular ensemble.
166 Parameters
167 ----------
168 n : int
169 Size of the matrix. It should be even for the classes C, CI, CII,
170 AII, DIII (either T^2 = -1 or P^2 = -1).
171 sym : one of 'A', 'AI', 'AII', 'AIII', 'BDI', 'CII', 'D', 'DIII', 'C', 'CI'
172 Altland-Zirnbauer symmetry class of the matrix.
173 charge : int or None
174 Topological invariant of the matrix. Should be one of 1, -1 in symmetry
175 classes D and DIII, should be from 0 to n in classes AIII and BDI,
176 and should be from 0 to n / 2 in class CII. If charge is None,
177 it is drawn from a binomial distribution with p = 1 / 2.
178 rng: int or rng (optional)
179 Seed or random number generator. If no 'rng' is passed, the random
180 number generator provided by numpy will be used.
182 Returns
183 -------
184 s : numpy.ndarray
185 A numpy array drawn from a corresponding circular ensemble.
187 Notes
188 -----
189 The representations of symmetry operators are chosen according to
190 Phys. Rev. B 85, 165409, except class D.
192 Matrix indices are grouped first according to channel number,
193 then sigma-index.
195 Chiral (sublattice) symmetry C always reads:
196 s = s^+.
197 Time reversal symmetry T reads:
198 AI, BDI: r = r^T.
199 CI: r = -sigma_y r^T sigma_y.
200 AII, DIII: r = -r^T.
201 CII: r = sigma_y r^T sigma_y.
202 Particle-hole symmetry reads:
203 CI: r = -sigma_y r^* sigma_y
204 C, CII: r = sigma_y r^* sigma_y
205 D, BDI: r = r^*.
206 DIII: -r = r^*.
208 This function uses QR decomposition to probe symmetric compact groups,
209 as detailed in arXiv:math-ph/0609050. For a reason as yet unknown, scipy
210 implementation of QR decomposition also works for symplectic matrices.
211 """
212 rng = ensure_rng(rng)
213 randn = rng.randn
215 if sym not in sym_list: 215 ↛ 216line 215 didn't jump to line 216, because the condition on line 215 was never true
216 raise ValueError('Unknown symmetry type.')
217 if (t(sym) == -1 or p(sym) == -1) and n % 2:
218 raise ValueError('n must be even in chosen symmetry class.')
220 # Prepare a real, complex or symplectic Gaussian matrix
221 # Real case.
222 if p(sym) == 1:
223 h = randn(n, n)
224 # Complex case.
225 elif p(sym) == 0:
226 h = randn(n, n) + 1j * randn(n, n)
227 # Symplectic case.
228 else: # p(sym) == -1
229 h = randn(n, n) + 1j * randn(n, n)
230 tau_z = np.array((n // 2) * [1, -1])
231 idx_x = np.arange(n) + tau_z
232 h -= tau_z.reshape(-1, 1) * h[idx_x][:, idx_x].conj() * tau_z
233 h *= 1j
235 # Generate a random matrix with proper electron-hole symmetry.
236 # This matrix is a random element of groups O(N), U(N), or USp(N).
237 s, h = qr(h)
238 if p(sym) == -1 and not np.allclose(np.diag(h, 1)[::2], 0, atol=1e-8): 238 ↛ 239line 238 didn't jump to line 239, because the condition on line 238 was never true
239 raise RuntimeError('QR decomposition symmetry failure.')
240 h = np.diag(h).copy()
241 h /= np.abs(h)
242 s *= h
244 # Ensure proper topological invariant in classes D and DIII.
245 if sym in ('D', 'DIII') and charge is not None: 245 ↛ 246line 245 didn't jump to line 246, because the condition on line 245 was never true
246 if charge not in (-1, 1):
247 raise ValueError('Impossible value of topological invariant.')
248 det = np.linalg.det(s)
249 if sym == 'DIII':
250 det *= (-1) ** (n // 2)
251 if (charge > 0) != (det > 0):
252 idx = np.arange(n)
253 idx[-1] -= 1
254 idx[-2] += 1
255 s = s[idx]
257 # Add the proper time-reversal symmetry:
258 if sym in ('AI', 'CI'):
259 s = np.dot(s.T, s)
260 if sym == 'CI':
261 tau_z = np.array((n // 2) * [1, -1])
262 idx_x = np.arange(n) + tau_z
263 s = 1j * tau_z * s[:, idx_x]
264 elif sym == 'AII' or sym == 'DIII':
265 tau_z = np.array((n // 2) * [1, -1])
266 idx_x = np.arange(n) + tau_z
267 s = 1j * np.dot(s.T * tau_z, s[idx_x])
269 # Add the chiral symmetry:
270 elif sym in ('AIII', 'BDI', 'CII'):
271 if sym != 'CII':
272 if charge is None: 272 ↛ 274line 272 didn't jump to line 274, because the condition on line 272 was never false
273 diag = 2 * rng.randint(2, size=(n,)) - 1
274 elif (0 <= charge <= n) and int(charge) == charge:
275 diag = np.array(charge * [-1] + (n - charge) * [1])
276 else:
277 raise ValueError('Impossible value of topological invariant.')
278 else:
279 if charge is None: 279 ↛ 282line 279 didn't jump to line 282, because the condition on line 279 was never false
280 diag = 2 * rng.randint(2, size=(n // 2,)) - 1
281 diag = np.resize(diag, (2, n // 2)).T.flatten()
282 elif (0 <= charge <= n // 2) and int(charge) == charge:
283 charge *= 2
284 diag = np.array(charge * [-1] + (n - charge) * [1])
285 else:
286 raise ValueError('Impossible value of topological invariant.')
288 s = np.dot(diag * s.T.conj(), s)
290 return s