Hide keyboard shortcuts

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. 

8 

9__all__ = ['DiscreteSymmetry'] 

10 

11import numpy as np 

12from scipy.sparse import identity, csr_matrix, hstack 

13 

14 

15def almost_identity(mat): 

16 return np.all(abs(mat - identity(mat.shape[0])).data < 1e-10) 

17 

18 

19def _column_sum(matrix): 

20 """Sum the columns of a sparse matrix. 

21 

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]) 

29 

30 

31def cond_conj(op, conj): 

32 return op.conj() if conj else op 

33 

34 

35_conj = [True, True, False] 

36_names = ['Time reversal', 'Particle-hole', 'Chiral'] 

37_signs = [1, -1, -1] 

38 

39class DiscreteSymmetry: 

40 r"""A collection of discrete symmetries and conservation laws. 

41 

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. 

52 

53 Notes 

54 ----- 

55 

56 When computing scattering modes, the representation of the 

57 modes is chosen to reflect declared discrete symmetries and 

58 conservation laws. 

59 

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`. 

66 

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. 

81 

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.") 

105 

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.") 

112 

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.") 

118 

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)) 

129 

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 

143 

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.") 

151 

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 

170 

171 def validate(self, matrix): 

172 """Check if a matrix satisfies the discrete symmetries. 

173 

174 Parameters 

175 ---------- 

176 matrix : sparse or dense matrix 

177 If rectangular, it is padded by zeros to be square. 

178 

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 

214 

215 def __getitem__(self, item): 

216 return (self.projectors, self.time_reversal, 

217 self.particle_hole, self.chiral)[item]