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-2016 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__ = ['DiscreteSymmetry']
11import numpy as np
12from scipy.sparse import identity, csr_matrix, hstack
15def almost_identity(mat):
16 return np.all(abs(mat - identity(mat.shape[0])).data < 1e-10)
19def _column_sum(matrix):
20 """Sum the columns of a sparse matrix.
22 This is fully analogous to ``matrix.sum(0)``, and uses an implementation
23 similar to that in scipy v1.1.0, however it avoids using ``numpy.matrix``
24 interface and therefore does not raise a ``PendingDeprecationWarning``.
25 This should be removed once we depend on scipy v1.2.0, where the warning is
26 silenced.
27 """
28 return matrix.T @ np.ones(matrix.shape[0])
31def cond_conj(op, conj):
32 return op.conj() if conj else op
35_conj = [True, True, False]
36_names = ['Time reversal', 'Particle-hole', 'Chiral']
37_signs = [1, -1, -1]
39class DiscreteSymmetry:
40 r"""A collection of discrete symmetries and conservation laws.
42 Parameters
43 ----------
44 projectors : iterable of rectangular or square sparse matrices
45 Projectors that block-diagonalize the Hamiltonian.
46 time_reversal : square sparse matrix
47 The unitary part of the time-reversal symmetry operator.
48 particle_hole : square sparse matrix
49 The unitary part of the particle-hole symmetry operator.
50 chiral : square sparse matrix
51 The chiral symmetry operator.
53 Notes
54 -----
56 When computing scattering modes, the representation of the
57 modes is chosen to reflect declared discrete symmetries and
58 conservation laws.
60 `projectors` block diagonalize the Hamiltonian, and modes are computed
61 separately in each block. The ordering of blocks is the same as of
62 `projectors`. If `conservation_law` is declared in
63 `~kwant.builder.Builder`, `projectors` is computed as the projectors
64 onto its orthogonal eigensubspaces. The projectors are stored in the
65 order of ascending eigenvalues of `conservation_law`.
67 Symmetrization using discrete symmetries varies depending on whether
68 a conservation law/projectors are declared. Consider the case with no
69 conservation law declared. With `time_reversal` declared, the outgoing
70 modes are chosen as the time-reversed partners of the incoming modes,
71 i.e. :math:`\psi_{out}(-k) = T \psi_{in}(k)` with k the momentum.
72 `chiral` also relates incoming and outgoing modes, such that
73 :math:`\psi_{out}(k) = C \psi_{in}(k)`. `particle_hole` gives symmetric
74 incoming and outgoing modes separately, such that
75 :math:`\psi_{in/out}(-k) = P \psi_{in/out}(k)`, except when k=-k, at
76 k = 0 or :math:`\pi`. In this case, each mode is chosen as an eigenstate
77 :math:`P \psi = \psi` if :math:`P^2=1`. If :math:`P^2=-1`, we
78 symmetrize the modes by generating pairs of orthogonal modes
79 :math:`\psi` and :math:`P\psi`. Because `chiral` and `particle_hole`
80 flip the sign of energy, they only apply at zero energy.
82 Discrete symmetries can be combined with a conservation law if they
83 leave each block invariant or transform it to another block. With S
84 a discrete symmetry and :math:`P_i` and :math:`P_j` projectors onto
85 blocks i and j of the Hamiltonian, :math:`S_{ji} = P_j^+ S P_i` is the
86 symmetry projection that maps from block i to block j.
87 :math:`S_{ji} = P_j^+ S P_i` must for each j be nonzero for exactly
88 one i. If S leaves block i invariant, the modes within block i are
89 symmetrized using the nonzero projection :math:`S_{ii}`, like in the
90 case without a conservation law. If S transforms between blocks i and
91 j, the modes of the block with the larger index are obtained by
92 transforming the modes of the block with the lower index. Thus, with
93 :math:`\psi_i` and :math:`\psi_j` the modes of blocks i and j, we have
94 :math:`\psi_j = S_{ji} \psi_i`.
95 """
96 def __init__(self, projectors=None, time_reversal=None, particle_hole=None,
97 chiral=None):
98 # Normalize format
99 if projectors is not None:
100 try:
101 projectors = [proj.tocsr() for proj in projectors]
102 except AttributeError:
103 raise TypeError("projectors must be a sequence of "
104 "sparse matrices.")
106 symms = (time_reversal, particle_hole, chiral)
107 try:
108 symms = [symm.tocsr() if symm is not None else None
109 for symm in symms]
110 except AttributeError:
111 raise TypeError("Symmetries must be sparse matrices.")
113 if None not in (time_reversal, particle_hole, chiral):
114 product = chiral.dot(particle_hole).dot(time_reversal.conj())
115 if not almost_identity(product):
116 raise ValueError("Product of all symmetries "
117 "should be equal to identity.")
119 # Auto-compute a missing symmetry.
120 if None not in (time_reversal, particle_hole) and chiral is None:
121 chiral = particle_hole.dot(time_reversal.conj())
122 elif None not in (particle_hole, chiral) and time_reversal is None:
123 # Product of particle hole with itself to get the sign right.
124 time_reversal = (particle_hole.dot(particle_hole.conj()).dot(
125 particle_hole).dot(chiral.conj()))
126 elif None not in (chiral, time_reversal) and particle_hole is None:
127 particle_hole = (time_reversal.dot(time_reversal.conj()).dot(
128 chiral).dot(time_reversal))
130 symms = [time_reversal, particle_hole, chiral]
131 for i, symm, conj, name in zip(range(3), symms, _conj, _names):
132 if symm is None:
133 continue
134 if not almost_identity(symm.T.conj().dot(symm)):
135 raise ValueError("{} symmetry should be "
136 "{}unitary.".format(name, "anti" * conj))
137 symm_squared = symm.dot(cond_conj(symm, conj))
138 if not (almost_identity(symm_squared) or
139 (conj and almost_identity(-symm_squared))):
140 raise ValueError("{} symmetry should square to "
141 "{}1.".format(name, "±" * conj))
142 symms[i] = symm
144 if projectors is not None:
145 projectors = [
146 p[:, _column_sum(abs(p)) > 1e-10] for p in projectors
147 ]
148 if not almost_identity(sum(projector.dot(projector.conj().T) for
149 projector in projectors)):
150 raise ValueError("Projectors should stack to a unitary matrix.")
152 # Check that the symmetries and conservation law are in canonical
153 # form. In the projector basis, each discrete symmetry should have
154 # at most one nonzero block per row.
155 for (symm, name, conj) in zip(symms, _names, _conj):
156 if symm is None:
157 continue
158 for p_i in projectors:
159 p_hc = p_i.T.conj()
160 nonzero_blocks = sum(
161 np.any(abs(p_hc.dot(symm).dot(
162 (cond_conj(p_j, conj))).data > 1e-7))
163 for p_j in projectors)
164 if nonzero_blocks > 1:
165 raise ValueError(name + ' symmetry is not in canonical '
166 'form: see DiscreteSymmetry '
167 'documentation.')
168 self.projectors = projectors
169 self.time_reversal, self.particle_hole, self.chiral = symms
171 def validate(self, matrix):
172 """Check if a matrix satisfies the discrete symmetries.
174 Parameters
175 ----------
176 matrix : sparse or dense matrix
177 If rectangular, it is padded by zeros to be square.
179 Returns
180 -------
181 broken_symmetries : list
182 List of strings, the names of symmetries broken by the
183 matrix: any combination of "Conservation law", "Time reversal",
184 "Particle-hole", "Chiral". If no symmetries are broken, returns
185 an empty list.
186 """
187 # Extra transposes are to enforse sparse dot product in case matrix is
188 # dense.
189 n, m = matrix.shape
190 if n != m:
191 if isinstance(matrix, np.ndarray):
192 new_matrix = np.empty((n, n), dtype=matrix.dtype)
193 new_matrix[:, :m] = matrix
194 new_matrix[:, m:] = 0
195 matrix = new_matrix
196 else:
197 matrix = hstack([matrix, csr_matrix((n, n-m))])
198 broken_symmetries = []
199 if self.projectors is not None:
200 for proj in self.projectors:
201 full = proj.dot(proj.T.conj())
202 commutator = full.dot(matrix) - (full.T.dot(matrix.T)).T
203 if np.linalg.norm(commutator.data) > 1e-8:
204 broken_symmetries.append('Conservation law')
205 break
206 for symm, conj, sign, name in zip(self[1:], _conj, _signs, _names):
207 if symm is None:
208 continue
209 commutator = symm.T.conj().dot((symm.T.dot(matrix.T)).T)
210 commutator = commutator - sign * cond_conj(matrix, conj)
211 if np.linalg.norm(commutator.data) > 1e-8:
212 broken_symmetries.append(name)
213 return broken_symmetries
215 def __getitem__(self, item):
216 return (self.projectors, self.time_reversal,
217 self.particle_hole, self.chiral)[item]