Coverage for icet/tools/structure_enumeration_support/normal_form_matrices.py: 98%
155 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
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]:
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_row(A.T, L.T, 0)
207 A, L = A.T, L.T
209 # Do the same thing for lower 2x2 matrix
210 while sorted(A[1, 1:])[0] != 0 or sorted(A[1:, 1])[0] != 0:
211 A, R = _gcd_reduce_row(A, R, 1)
212 A, L = _gcd_reduce_row(A.T, L.T, 1)
213 A, L = A.T, L.T
215 # If last diagonal entry is negative,
216 # make it positive
217 if A[2, 2] < 0:
218 A[2, 2] = -A[2, 2]
219 L[2] = -L[2]
221 # Check that the diagonal entry i,i divides
222 # diagonal entry i+1, i+1. Otherwise,
223 # add row i+1 to i and start over.
224 if A[2, 2] % A[1, 1] != 0:
225 A[1] = A[1] + A[2]
226 L[1] = L[1] + L[2]
227 elif A[1, 1] % A[0, 0] != 0: 227 ↛ 231line 227 didn't jump to line 231, because the condition on line 227 was never false
228 A[0] = A[0] + A[1]
229 L[0] = L[0] + L[1]
230 else:
231 break
232 assert (abs(np.dot(np.dot(L, H), R) - A) < 1e-5).all()
233 self.S_matrix = A
234 self.L = L
236 def add_hnf(self, hnf: HermiteNormalForm) -> None:
237 """Add HNF to SNF.
239 Paramaters
240 ----------
241 hnf
242 HermiteNormalForm object
243 """
244 self.hnfs.append(hnf)
246 def set_group_order(self) -> None:
247 """
248 Set group representation of an SNF matrix (the G matrix in HarFor08).
249 """
250 group_order = []
251 for i in range(self.S[0]):
252 for j in range(self.S[1]):
253 for k in range(self.S[2]):
254 group_order.append([i, j, k])
255 self.group_order = group_order
258def _gcd_reduce_row(A: np.ndarray, R: np.ndarray, i: int) -> Tuple[np.ndarray, np.ndarray]:
259 """
260 Use column operations to make A[i, i] the greatest common
261 denominator of the elements in row i and the other elements
262 zero.
264 Parameters
265 ----------
266 A
267 Matrix whose row is to be cleared.
268 R
269 Matrix that should be subject to the same operations.
270 i
271 Index of row to be treated.
273 Returns
274 -------
275 A
276 Treated matrix A.
277 R
278 Matrix that has been subject to the same operations.
279 """
280 for j in range(i, 3):
281 if A[i, j] < 0:
282 A[:, j] = -1 * A[:, j]
283 R[:, j] = -1 * R[:, j]
284 while np.sort(A[i, i:])[1 - i] > 0:
285 max_index = np.argmax(A[i, i:]) + i
286 min_index = np.argmin(A[i, i:]) + i
287 if max_index == min_index:
288 max_index += 1
289 if A[i, min_index] == 0 and i == 0:
290 if np.sort(A[i])[1] > 0: 290 ↛ 296line 290 didn't jump to line 296, because the condition on line 290 was never false
291 min_index += 1
292 min_index = min_index % 3
293 if min_index == max_index:
294 min_index += 1
295 min_index = min_index % 3
296 if A[i, min_index] == A[i, max_index]:
297 tmp = min_index
298 min_index = max_index
299 max_index = tmp
300 A[:, max_index] = A[:, max_index] - A[:, min_index]
301 R[:, max_index] = R[:, max_index] - R[:, min_index]
302 max_index = np.argmax(A[i])
303 return A, R
306def get_unique_snfs(hnfs: List[HermiteNormalForm]) -> List[SmithNormalForm]:
307 """
308 For a list of Hermite Normal Forms, obtain the set of unique Smith Normal
309 Forms.
311 Parameters
312 ----------
313 hnfs
314 List of HermiteNormalForm objects
316 Returns
317 -------
318 The unique Smith Normal Form matrices in a list
319 """
320 snfs = []
321 for hnf in hnfs:
322 # Check whether the snf is new or already encountered
323 snf = hnf.snf
324 snf_is_new = True
325 for snf_comp in snfs:
326 if snf_comp.S == snf.S:
327 snf_is_new = False
328 snf_comp.add_hnf(hnf)
329 break
330 if snf_is_new:
331 snf.set_group_order()
332 snf.add_hnf(hnf)
333 snfs.append(snf)
334 return snfs