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

2Handling of Hermite Normal Form and Smith Normal Form matrices 

3""" 

4 

5import numpy as np 

6from typing import List, Tuple, Dict 

7 

8 

9class HermiteNormalForm(object): 

10 """ 

11 Hermite Normal Form matrix. 

12 """ 

13 

14 def __init__(self, H: np.ndarray, rotations: np.ndarray, 

15 translations: np.ndarray, basis_shifts: np.ndarray): 

16 self.H = H 

17 self.snf = SmithNormalForm(H) 

18 self.transformations = [] 

19 self.compute_transformations(rotations, translations, basis_shifts) 

20 

21 def compute_transformations(self, rotations: np.ndarray, 

22 translations: np.ndarray, 

23 basis_shifts: np.ndarray, 

24 tolerance: float = 1e-3) -> None: 

25 """ 

26 Save transformations (based on rotations) that turns the supercell 

27 into an equivalent supercell. Precompute these transformations, 

28 consisting of permutation as well as translation and basis shift, for 

29 later use. 

30 

31 Parameters 

32 ---------- 

33 rotations 

34 list of (rotational) symmetry operations 

35 translations 

36 list of translations that go together with the rotational operations 

37 basis_shifts 

38 corresponding to d_Nd in HarFor09 

39 tolerance 

40 tolerance applied when checking that matrices are all integers 

41 """ 

42 

43 for R, T, basis_shift in zip(rotations, translations, 

44 basis_shifts): 

45 check = np.dot(np.dot(np.linalg.inv(self.H), R), self.H) 

46 check = check - np.round(check) 

47 if (abs(check) < tolerance).all(): 

48 LRL = np.dot(self.snf.L, np.dot(R, np.linalg.inv(self.snf.L))) 

49 

50 # Should be an integer matrix 

51 assert (abs(LRL - np.round(LRL)) < tolerance).all() 

52 LRL = np.round(LRL).astype(np.int64) 

53 LT = np.dot(T, self.snf.L.T) 

54 self.transformations.append([LRL, LT, basis_shift]) 

55 

56 

57def yield_hermite_normal_forms(det: int, pbc: List[bool]) -> np.ndarray: 

58 """ 

59 Yield all Hermite Normal Form matrices with determinant det. 

60 

61 Parameters 

62 ---------- 

63 det 

64 Target determinant of HNFs 

65 pbc 

66 Periodic boundary conditions of the primitive structure 

67 

68 Yields 

69 ------ 

70 3x3 HNF matrix 

71 """ 

72 # 1D 

73 if sum(pbc) == 1: 

74 hnf = np.eye(3, dtype=int) 

75 for i, bc in enumerate(pbc): 75 ↛ 79line 75 didn't jump to line 79, because the loop on line 75 didn't complete

76 if bc: 

77 hnf[i, i] = det 

78 break 

79 yield hnf 

80 

81 # 2D 

82 elif sum(pbc) == 2: 

83 for a in range(1, det + 1): 

84 if det % a == 0: 

85 c = det // a 

86 for b in range(0, c): 

87 if not pbc[0]: 87 ↛ 88line 87 didn't jump to line 88, because the condition on line 87 was never true

88 hnf = [[1, 0, 0], 

89 [0, a, 0], 

90 [0, b, c]] 

91 elif not pbc[1]: 91 ↛ 92line 91 didn't jump to line 92, because the condition on line 91 was never true

92 hnf = [[a, 0, 0], 

93 [0, 1, 0], 

94 [b, 0, c]] 

95 else: 

96 hnf = [[a, 0, 0], 

97 [b, c, 0], 

98 [0, 0, 1]] 

99 yield np.array(hnf) 

100 

101 # 3D 

102 else: 

103 for a in range(1, det + 1): 

104 if det % a == 0: 

105 for c in range(1, det // a + 1): 

106 if det // a % c == 0: 

107 f = det // (a * c) 

108 for b in range(0, c): 

109 for d in range(0, f): 

110 for e in range(0, f): 

111 hnf = [[a, 0, 0], 

112 [b, c, 0], 

113 [d, e, f]] 

114 yield np.array(hnf) 

115 

116 

117def yield_reduced_hnfs(ncells: int, symmetries: Dict, 

118 pbc: List[bool], tolerance: float = 1e-3) -> HermiteNormalForm: 

119 """ 

120 For a fixed determinant N (i.e., a number of atoms N), yield all 

121 Hermite Normal Forms (HNF) that are inequivalent under symmetry 

122 operations of the parent lattice.' 

123 

124 Parameters 

125 ---------- 

126 ncells 

127 Determinant of the HNF. 

128 symmetries 

129 Symmetry operations of the parent lattice. 

130 pbc 

131 Periodic boundary conditions of the primitive structure 

132 tolerance 

133 tolerance applied when checking that matrices are all integers 

134 

135 Yields 

136 ------ 

137 HermiteNormalForm object, each one symmetrically distinct from the others 

138 """ 

139 rotations = symmetries['rotations'] 

140 translations = symmetries['translations'] 

141 basis_shifts = symmetries['basis_shifts'] 

142 hnfs = [] 

143 

144 for hnf in yield_hermite_normal_forms(ncells, pbc): 

145 

146 # Throw away HNF:s that yield equivalent supercells 

147 hnf_inv = np.linalg.inv(hnf) 

148 duplicate = False 

149 for R in rotations: 

150 HR = np.dot(hnf_inv, R) 

151 for hnf_previous in hnfs: 

152 check = np.dot(HR, hnf_previous.H) 

153 check = check - np.round(check) 

154 if (abs(check) < tolerance).all(): 

155 duplicate = True 

156 break 

157 if duplicate: 

158 break 

159 if duplicate: 

160 continue 

161 

162 # If it's not a duplicate, save the hnf 

163 # and the supercell so that it can be compared to 

164 hnf = HermiteNormalForm(hnf, rotations, translations, basis_shifts) 

165 hnfs.append(hnf) 

166 yield hnf 

167 

168 

169class SmithNormalForm(object): 

170 """ 

171 Smith Normal Form matrix. 

172 """ 

173 

174 def __init__(self, H: np.ndarray): 

175 self.compute_snf(H) 

176 self.S = tuple([self.S_matrix[i, i] for i in range(3)]) 

177 self.ncells = self.S[0] * self.S[1] * self.S[2] 

178 self.group_order = None 

179 self.hnfs = [] 

180 

181 # Help list for permuting labelings 

182 blocks = [self.ncells // self.S[0]] 

183 for i in range(1, 3): 

184 blocks.append(blocks[-1] // self.S[i]) 

185 self.blocks = blocks 

186 

187 def compute_snf(self, H: np.ndarray) -> None: 

188 """ 

189 Compute Smith Normal Form for 3x3 matrix. Note that H = L*S*R. 

190 

191 Parameters 

192 ---------- 

193 H 

194 3x3 matrix 

195 """ 

196 A = H.copy() 

197 L = np.eye(3, dtype=int) 

198 R = np.eye(3, dtype=int) 

199 while True: 

200 # Clear upper row and leftmost column in such 

201 # a way that greatest common denominator ends 

202 # up in A[0, 0], in a standard Smith Normal Form way 

203 # (Euclidean algorithm for finding greatest common divisor) 

204 while sorted(A[0])[1] != 0 or sorted(A[:, 0])[1] != 0: 

205 A, R = _gcd_reduce_row(A, R, 0) 

206 A, L = _gcd_reduce_column(A, L, 0) 

207 

208 # Do the same thing for lower 2x2 matrix 

209 while sorted(A[1, 1:])[0] != 0 or sorted(A[1:, 1])[0] != 0: 

210 A, R = _gcd_reduce_row(A, R, 1) 

211 A, L = _gcd_reduce_column(A, L, 1) 

212 

213 # If last diagonal entry is negative, 

214 # make it positive 

215 if A[2, 2] < 0: 

216 A[2, 2] = -A[2, 2] 

217 L[2] = -L[2] 

218 

219 # Check that the diagonal entry i,i divides 

220 # diagonal entry i+1, i+1. Otherwise, 

221 # add row i+1 to i and start over. 

222 if A[2, 2] % A[1, 1] != 0: 

223 A[1] = A[1] + A[2] 

224 L[1] = L[1] + L[2] 

225 elif A[1, 1] % A[0, 0] != 0: 

226 A[0] = A[0] + A[1] 

227 L[0] = L[0] + L[1] 

228 else: 

229 break 

230 assert (abs(np.dot(np.dot(L, H), R) - A) < 1e-5).all() 

231 self.S_matrix = A 

232 self.L = L 

233 

234 def add_hnf(self, hnf: HermiteNormalForm) -> None: 

235 """Add HNF to SNF. 

236 

237 Paramaters 

238 ---------- 

239 hnf 

240 HermiteNormalForm object 

241 """ 

242 self.hnfs.append(hnf) 

243 

244 def set_group_order(self) -> None: 

245 """ 

246 Set group representation of an SNF matrix (the G matrix in HarFor08). 

247 """ 

248 group_order = [] 

249 for i in range(self.S[0]): 

250 for j in range(self.S[1]): 

251 for k in range(self.S[2]): 

252 group_order.append([i, j, k]) 

253 self.group_order = group_order 

254 

255 

256def _switch_rows(A: np.ndarray, i: int, j: int) -> np.ndarray: 

257 """ 

258 Switch rows in matrix. 

259 

260 Parameters 

261 --------- 

262 A 

263 Matrix in which rows will be swapped. 

264 i 

265 Index of row 1 to be swapped. 

266 j 

267 Index of row 2 to be swapped. 

268 

269 Returns 

270 ------- 

271 Matrix with swapped rows. 

272 """ 

273 row = A[j].copy() 

274 A[j] = A[i] 

275 A[i] = row 

276 return A 

277 

278 

279def _switch_columns(A: np.ndarray, i: int, j: int): 

280 """ 

281 Switch columns in matrix. 

282 

283 Parameters 

284 --------- 

285 A : ndarray 

286 Matrix in which columns will be swapped. 

287 i : int 

288 Index of column 1 to be swapped. 

289 j : int 

290 Index of column 2 to be swapped. 

291 

292 Returns 

293 ------- 

294 Matrix with swapped columns. 

295 """ 

296 col = A[:, j].copy() 

297 A[:, j] = A[:, i] 

298 A[:, i] = col 

299 return A 

300 

301 

302def _gcd_reduce_row(A: np.ndarray, R: np.ndarray, i: int) -> Tuple[np.ndarray, np.ndarray]: 

303 """ 

304 Use column operations to make A[i, i] the greatest common 

305 denominator of the elements in row i and the other elements 

306 zero. 

307 

308 Parameters 

309 ---------- 

310 A 

311 Matrix whose row is to be cleared. 

312 R 

313 Matrix that should be subject to the same operations. 

314 i 

315 Index of row to be treated. 

316 

317 Returns 

318 ------- 

319 A 

320 Treated matrix A. 

321 R 

322 Matrix that has been subject to the same operations. 

323 """ 

324 for j in range(i, 3): 

325 if A[i, j] < 0: 

326 A[:, j] = -1 * A[:, j] 

327 R[:, j] = -1 * R[:, j] 

328 while np.sort(A[i, i:])[1 - i] > 0: 

329 max_index = np.argmax(A[i, i:]) + i 

330 min_index = np.argmin(A[i, i:]) + i 

331 if max_index == min_index: 

332 max_index += 1 

333 if A[i, min_index] == 0 and i == 0: 

334 if np.sort(A[i])[1] > 0: 334 ↛ 340line 334 didn't jump to line 340, because the condition on line 334 was never false

335 min_index += 1 

336 min_index = min_index % 3 

337 if min_index == max_index: 

338 min_index += 1 

339 min_index = min_index % 3 

340 if A[i, min_index] == A[i, max_index]: 

341 tmp = min_index 

342 min_index = max_index 

343 max_index = tmp 

344 A[:, max_index] = A[:, max_index] - A[:, min_index] 

345 R[:, max_index] = R[:, max_index] - R[:, min_index] 

346 max_index = np.argmax(A[i]) 

347 A = _switch_columns(A, i, max_index) 

348 R = _switch_columns(R, i, max_index) 

349 return A, R 

350 

351 

352def _gcd_reduce_column(A: np.ndarray, L: np.ndarray, j: int) -> Tuple[np.ndarray, np.ndarray]: 

353 """ 

354 Use row operations to make A[i, i] the greatest common 

355 denominator of the elements in column i and the other elements 

356 zero. 

357 

358 Parameters 

359 ---------- 

360 A 

361 Matrix whose column is to be cleared. 

362 L 

363 Matrix that should be subject to the same operations. 

364 i 

365 Index of column to be treated. 

366 

367 Returns 

368 ------- 

369 A 

370 Treated matrix A. 

371 L 

372 Matrix that has been subject to the same operations. 

373 """ 

374 for i in range(j, 3): 

375 if A[i, j] < 0: 

376 A[i] = -1 * A[i] 

377 L[i] = -1 * L[i] 

378 while np.sort(A[j:, j])[1 - j] > 0: 

379 max_index = np.argmax(A[j:, j]) + j 

380 min_index = np.argmin(A[j:, j]) + j 

381 if max_index == min_index: 

382 max_index += 1 

383 if A[min_index, j] == 0 and j == 0: 

384 if np.sort(A[:, j])[1] > 0: 384 ↛ 390line 384 didn't jump to line 390, because the condition on line 384 was never false

385 min_index += 1 

386 min_index = min_index % 3 

387 if min_index == max_index: 

388 min_index += 1 

389 min_index = min_index % 3 

390 if A[min_index, j] == A[max_index, j]: 

391 tmp = min_index 

392 min_index = max_index 

393 max_index = tmp 

394 A[max_index] = A[max_index] - A[min_index] 

395 L[max_index] = L[max_index] - L[min_index] 

396 max_index = np.argmax(A[:, j]) 

397 A = _switch_rows(A, j, max_index) 

398 L = _switch_rows(L, j, max_index) 

399 return A, L 

400 

401 

402def get_unique_snfs(hnfs: List[HermiteNormalForm]) -> List[SmithNormalForm]: 

403 """ 

404 For a list of Hermite Normal Forms, obtain the set of unique Smith Normal 

405 Forms. 

406 

407 Parameters 

408 ---------- 

409 hnfs 

410 List of HermiteNormalForm objects 

411 

412 Returns 

413 ------- 

414 The unique Smith Normal Form matrices in a list 

415 """ 

416 snfs = [] 

417 for hnf in hnfs: 

418 # Check whether the snf is new or already encountered 

419 snf = hnf.snf 

420 snf_is_new = True 

421 for snf_comp in snfs: 

422 if snf_comp.S == snf.S: 

423 snf_is_new = False 

424 snf_comp.add_hnf(hnf) 

425 break 

426 if snf_is_new: 

427 snf.set_group_order() 

428 snf.add_hnf(hnf) 

429 snfs.append(snf) 

430 return snfs