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

5import numpy as np

6from typing import List, Tuple, Dict

9class HermiteNormalForm(object):

10 """

11 Hermite Normal Form matrix.

12 """

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)

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.

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

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

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

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

58 """

59 Yield all Hermite Normal Form matrices with determinant det.

61 Parameters

62 ----------

63 det

64 Target determinant of HNFs

65 pbc

66 Periodic boundary conditions of the primitive structure

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

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)

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)

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

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

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 = []

144 for hnf in yield_hermite_normal_forms(ncells, pbc):

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

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

169class SmithNormalForm(object):

170 """

171 Smith Normal Form matrix.

172 """

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 = []

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

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

188 """

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

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)

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)

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]

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

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

237 Paramaters

238 ----------

239 hnf

240 HermiteNormalForm object

241 """

242 self.hnfs.append(hnf)

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

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

257 """

258 Switch rows in matrix.

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.

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

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

280 """

281 Switch columns in matrix.

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.

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

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.

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.

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

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.

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.

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

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.

407 Parameters

408 ----------

409 hnfs

410 List of HermiteNormalForm objects

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

425 break

426 if snf_is_new:

427 snf.set_group_order()

429 snfs.append(snf)

430 return snfs