Coverage for icet/tools/structure_enumeration_support/normal_form_matrices.py: 98%
154 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-09-14 04:08 +0000
1"""
2Handling of Hermite Normal Form and Smith Normal Form matrices
3"""
5import numpy as np
8class HermiteNormalForm(object):
9 """
10 Hermite Normal Form matrix.
11 """
13 def __init__(self, H: np.ndarray, rotations: np.ndarray,
14 translations: np.ndarray, basis_shifts: np.ndarray):
15 self.H = H
16 self.snf = SmithNormalForm(H)
17 self.transformations = []
18 self.compute_transformations(rotations, translations, basis_shifts)
20 def compute_transformations(self, rotations: np.ndarray,
21 translations: np.ndarray,
22 basis_shifts: np.ndarray,
23 tolerance: float = 1e-3) -> None:
24 """
25 Save transformations (based on rotations) that turns the supercell
26 into an equivalent supercell. Precompute these transformations,
27 consisting of permutation as well as translation and basis shift, for
28 later use.
30 Parameters
31 ----------
32 rotations
33 list of (rotational) symmetry operations
34 translations
35 list of translations that go together with the rotational operations
36 basis_shifts
37 corresponding to d_Nd in HarFor09
38 tolerance
39 tolerance applied when checking that matrices are all integers
40 """
42 for R, T, basis_shift in zip(rotations, translations,
43 basis_shifts):
44 check = np.dot(np.dot(np.linalg.inv(self.H), R), self.H)
45 check = check - np.round(check)
46 if (abs(check) < tolerance).all():
47 LRL = np.dot(self.snf.L, np.dot(R, np.linalg.inv(self.snf.L)))
49 # Should be an integer matrix
50 assert (abs(LRL - np.round(LRL)) < tolerance).all()
51 LRL = np.round(LRL).astype(np.int64)
52 LT = np.dot(T, self.snf.L.T)
53 self.transformations.append([LRL, LT, basis_shift])
56def yield_hermite_normal_forms(det: int, pbc: list[bool]) -> np.ndarray:
57 """
58 Yield all Hermite Normal Form matrices with determinant det.
60 Parameters
61 ----------
62 det
63 Target determinant of HNFs
64 pbc
65 Periodic boundary conditions of the primitive structure
67 Yields
68 ------
69 3x3 HNF matrix
70 """
71 # 1D
72 if sum(pbc) == 1:
73 hnf = np.eye(3, dtype=int)
74 for i, bc in enumerate(pbc): 74 ↛ 78line 74 didn't jump to line 78 because the loop on line 74 didn't complete
75 if bc:
76 hnf[i, i] = det
77 break
78 yield hnf
80 # 2D
81 elif sum(pbc) == 2:
82 for a in range(1, det + 1):
83 if det % a == 0:
84 c = det // a
85 for b in range(0, c):
86 if not pbc[0]:
87 hnf = [[1, 0, 0],
88 [0, a, 0],
89 [0, b, c]]
90 elif not pbc[1]: 90 ↛ 91line 90 didn't jump to line 91 because the condition on line 90 was never true
91 hnf = [[a, 0, 0],
92 [0, 1, 0],
93 [b, 0, c]]
94 else:
95 hnf = [[a, 0, 0],
96 [b, c, 0],
97 [0, 0, 1]]
98 yield np.array(hnf)
100 # 3D
101 else:
102 for a in range(1, det + 1):
103 if det % a == 0:
104 for c in range(1, det // a + 1):
105 if det // a % c == 0:
106 f = det // (a * c)
107 for b in range(0, c):
108 for d in range(0, f):
109 for e in range(0, f):
110 hnf = [[a, 0, 0],
111 [b, c, 0],
112 [d, e, f]]
113 yield np.array(hnf)
116def yield_reduced_hnfs(ncells: int, symmetries: dict,
117 pbc: list[bool], tolerance: float = 1e-3) -> HermiteNormalForm:
118 """
119 For a fixed determinant N (i.e., a number of atoms N), yield all
120 Hermite Normal Forms (HNF) that are inequivalent under symmetry
121 operations of the parent lattice.'
123 Parameters
124 ----------
125 ncells
126 Determinant of the HNF.
127 symmetries
128 Symmetry operations of the parent lattice.
129 pbc
130 Periodic boundary conditions of the primitive structure
131 tolerance
132 tolerance applied when checking that matrices are all integers
134 Yields
135 ------
136 HermiteNormalForm object, each one symmetrically distinct from the others
137 """
138 rotations = symmetries['rotations']
139 translations = symmetries['translations']
140 basis_shifts = symmetries['basis_shifts']
141 hnfs = []
143 for hnf in yield_hermite_normal_forms(ncells, pbc):
145 # Throw away HNF:s that yield equivalent supercells
146 hnf_inv = np.linalg.inv(hnf)
147 duplicate = False
148 for R in rotations:
149 HR = np.dot(hnf_inv, R)
150 for hnf_previous in hnfs:
151 check = np.dot(HR, hnf_previous.H)
152 check = check - np.round(check)
153 if (abs(check) < tolerance).all():
154 duplicate = True
155 break
156 if duplicate:
157 break
158 if duplicate:
159 continue
161 # If it's not a duplicate, save the hnf
162 # and the supercell so that it can be compared to
163 hnf = HermiteNormalForm(hnf, rotations, translations, basis_shifts)
164 hnfs.append(hnf)
165 yield hnf
168class SmithNormalForm(object):
169 """
170 Smith Normal Form matrix.
171 """
173 def __init__(self, H: np.ndarray):
174 self.compute_snf(H)
175 self.S = tuple([self.S_matrix[i, i] for i in range(3)])
176 self.ncells = self.S[0] * self.S[1] * self.S[2]
177 self.group_order = None
178 self.hnfs = []
180 # Help list for permuting labelings
181 blocks = [self.ncells // self.S[0]]
182 for i in range(1, 3):
183 blocks.append(blocks[-1] // self.S[i])
184 self.blocks = blocks
186 def compute_snf(self, H: np.ndarray) -> None:
187 """
188 Compute Smith Normal Form for 3x3 matrix. Note that H = L*S*R.
190 Parameters
191 ----------
192 H
193 3x3 matrix
194 """
195 A = H.copy()
196 L = np.eye(3, dtype=int)
197 R = np.eye(3, dtype=int)
198 while True:
199 # Clear upper row and leftmost column in such
200 # a way that greatest common denominator ends
201 # up in A[0, 0], in a standard Smith Normal Form way
202 # (Euclidean algorithm for finding greatest common divisor)
203 while sorted(A[0])[1] != 0 or sorted(A[:, 0])[1] != 0:
204 A, R = _gcd_reduce_row(A, R, 0)
205 A, L = _gcd_reduce_row(A.T, L.T, 0)
206 A, L = A.T, L.T
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_row(A.T, L.T, 1)
212 A, L = A.T, L.T
214 # If last diagonal entry is negative,
215 # make it positive
216 if A[2, 2] < 0:
217 A[2, 2] = -A[2, 2]
218 L[2] = -L[2]
220 # Check that the diagonal entry i,i divides
221 # diagonal entry i+1, i+1. Otherwise,
222 # add row i+1 to i and start over.
223 if A[2, 2] % A[1, 1] != 0:
224 A[1] = A[1] + A[2]
225 L[1] = L[1] + L[2]
226 elif A[1, 1] % A[0, 0] != 0: 226 ↛ 230line 226 didn't jump to line 230 because the condition on line 226 was always true
227 A[0] = A[0] + A[1]
228 L[0] = L[0] + L[1]
229 else:
230 break
231 assert (abs(np.dot(np.dot(L, H), R) - A) < 1e-5).all()
232 self.S_matrix = A
233 self.L = L
235 def add_hnf(self, hnf: HermiteNormalForm) -> None:
236 """Add HNF to SNF.
238 Paramaters
239 ----------
240 hnf
241 HermiteNormalForm object
242 """
243 self.hnfs.append(hnf)
245 def set_group_order(self) -> None:
246 """
247 Set group representation of an SNF matrix (the G matrix in HarFor08).
248 """
249 group_order = []
250 for i in range(self.S[0]):
251 for j in range(self.S[1]):
252 for k in range(self.S[2]):
253 group_order.append([i, j, k])
254 self.group_order = group_order
257def _gcd_reduce_row(A: np.ndarray, R: np.ndarray, i: int) -> tuple[np.ndarray, np.ndarray]:
258 """
259 Use column operations to make A[i, i] the greatest common
260 denominator of the elements in row i and the other elements
261 zero.
263 Parameters
264 ----------
265 A
266 Matrix whose row is to be cleared.
267 R
268 Matrix that should be subject to the same operations.
269 i
270 Index of row to be treated.
272 Returns
273 -------
274 A
275 Treated matrix A.
276 R
277 Matrix that has been subject to the same operations.
278 """
279 for j in range(i, 3):
280 if A[i, j] < 0:
281 A[:, j] = -1 * A[:, j]
282 R[:, j] = -1 * R[:, j]
283 while np.sort(A[i, i:])[1 - i] > 0:
284 max_index = np.argmax(A[i, i:]) + i
285 min_index = np.argmin(A[i, i:]) + i
286 if max_index == min_index:
287 max_index += 1
288 if A[i, min_index] == 0 and i == 0:
289 if np.sort(A[i])[1] > 0: 289 ↛ 295line 289 didn't jump to line 295 because the condition on line 289 was always true
290 min_index += 1
291 min_index = min_index % 3
292 if min_index == max_index:
293 min_index += 1
294 min_index = min_index % 3
295 if A[i, min_index] == A[i, max_index]:
296 tmp = min_index
297 min_index = max_index
298 max_index = tmp
299 A[:, max_index] = A[:, max_index] - A[:, min_index]
300 R[:, max_index] = R[:, max_index] - R[:, min_index]
301 max_index = np.argmax(A[i])
302 return A, R
305def get_unique_snfs(hnfs: list[HermiteNormalForm]) -> list[SmithNormalForm]:
306 """
307 For a list of Hermite Normal Forms, obtain the set of unique Smith Normal
308 Forms.
310 Parameters
311 ----------
312 hnfs
313 List of HermiteNormalForm objects.
315 Returns
316 -------
317 The unique Smith Normal Form matrices in a list.
318 """
319 snfs = []
320 for hnf in hnfs:
321 # Check whether the snf is new or already encountered
322 snf = hnf.snf
323 snf_is_new = True
324 for snf_comp in snfs:
325 if snf_comp.S == snf.S:
326 snf_is_new = False
327 snf_comp.add_hnf(hnf)
328 break
329 if snf_is_new:
330 snf.set_group_order()
331 snf.add_hnf(hnf)
332 snfs.append(snf)
333 return snfs