Coverage for ase / geometry / minkowski_reduction.py: 96.95%
131 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3import itertools
5import numpy as np
7from ase.cell import Cell
8from ase.utils import pbc2pbc
10TOL = 1E-12
11MAX_IT = 100000 # in practice this is not exceeded
14class CycleChecker:
16 def __init__(self, d):
17 assert d in [2, 3]
19 # worst case is the hexagonal cell in 2D and the fcc cell in 3D
20 n = {2: 6, 3: 12}[d]
22 # max cycle length is total number of primtive cell descriptions
23 max_cycle_length = np.prod([n - i for i in range(d)]) * np.prod(d)
24 self.visited = np.zeros((max_cycle_length, 3 * d), dtype=int)
26 def add_site(self, H):
27 # flatten array for simplicity
28 H = H.ravel()
30 # check if site exists
31 found = (self.visited == H).all(axis=1).any()
33 # shift all visited sites down and place current site at the top
34 self.visited = np.roll(self.visited, 1, axis=0)
35 self.visited[0] = H
36 return found
39def reduction_gauss(B, hu, hv):
40 """Calculate a Gauss-reduced lattice basis (2D reduction)."""
41 cycle_checker = CycleChecker(d=2)
42 u = hu @ B
43 v = hv @ B
45 for _ in range(MAX_IT):
46 x = int(round(np.dot(u, v) / np.dot(u, u)))
47 hu, hv = hv - x * hu, hu
48 u = hu @ B
49 v = hv @ B
50 site = np.array([hu, hv])
51 if np.dot(u, u) >= np.dot(v, v) or cycle_checker.add_site(site):
52 return hv, hu
54 raise RuntimeError(f"Gaussian basis not found after {MAX_IT} iterations")
57def relevant_vectors_2D(u, v):
58 cs = np.array([e for e in itertools.product([-1, 0, 1], repeat=2)])
59 vs = cs @ [u, v]
60 indices = np.argsort(np.linalg.norm(vs, axis=1))[:7]
61 return vs[indices], cs[indices]
64def closest_vector(t0, u, v):
65 t = t0
66 a = np.zeros(2, dtype=int)
67 rs, cs = relevant_vectors_2D(u, v)
69 dprev = float("inf")
70 for _ in range(MAX_IT):
71 ds = np.linalg.norm(rs + t, axis=1)
72 index = np.argmin(ds)
73 if index == 0 or ds[index] >= dprev:
74 return a
76 dprev = ds[index]
77 r = rs[index]
78 kopt = int(round(-np.dot(t, r) / np.dot(r, r)))
79 a += kopt * cs[index]
80 t = t0 + a[0] * u + a[1] * v
82 raise RuntimeError(f"Closest vector not found after {MAX_IT} iterations")
85def reduction_full(B):
86 """Calculate a Minkowski-reduced lattice basis (3D reduction)."""
87 cycle_checker = CycleChecker(d=3)
88 H = np.eye(3, dtype=int)
89 norms = np.linalg.norm(B, axis=1)
91 for _ in range(MAX_IT):
92 # Sort vectors by norm
93 H = H[np.argsort(norms, kind='merge')]
95 # Gauss-reduce smallest two vectors
96 hw = H[2]
97 hu, hv = reduction_gauss(B, H[0], H[1])
98 H = np.array([hu, hv, hw])
99 R = H @ B
101 # Orthogonalize vectors using Gram-Schmidt
102 u, v, _ = R
103 X = u / np.linalg.norm(u)
104 Y = v - X * np.dot(v, X)
105 Y /= np.linalg.norm(Y)
107 # Find closest vector to last element of R
108 pu, pv, pw = R @ np.array([X, Y]).T
109 nb = closest_vector(pw, pu, pv)
111 # Update basis
112 H[2] = [nb[0], nb[1], 1] @ H
113 R = H @ B
115 norms = np.linalg.norm(R, axis=1)
116 if norms[2] >= norms[1] or cycle_checker.add_site(H):
117 return R, H
119 raise RuntimeError(f"Reduced basis not found after {MAX_IT} iterations")
122def is_minkowski_reduced(cell, pbc=True):
123 """Tests if a cell is Minkowski-reduced.
125 Parameters
126 ----------
128 cell: array
129 The lattice basis to test (in row-vector format).
130 pbc: array, optional
131 The periodic boundary conditions of the cell (Default `True`).
132 If `pbc` is provided, only periodic cell vectors are tested.
134 Returns
135 -------
137 is_reduced: bool
138 True if cell is Minkowski-reduced, False otherwise.
139 """
141 """These conditions are due to Minkowski, but a nice description in English
142 can be found in the thesis of Carine Jaber: "Algorithmic approaches to
143 Siegel's fundamental domain", https://www.theses.fr/2017UBFCK006.pdf
144 This is also good background reading for Minkowski reduction.
146 0D and 1D cells are trivially reduced. For 2D cells, the conditions which
147 an already-reduced basis fulfil are:
148 |b1| ≤ |b2|
149 |b2| ≤ |b1 - b2|
150 |b2| ≤ |b1 + b2|
152 For 3D cells, the conditions which an already-reduced basis fulfil are:
153 |b1| ≤ |b2| ≤ |b3|
155 |b1 + b2| ≥ |b2|
156 |b1 + b3| ≥ |b3|
157 |b2 + b3| ≥ |b3|
158 |b1 - b2| ≥ |b2|
159 |b1 - b3| ≥ |b3|
160 |b2 - b3| ≥ |b3|
161 |b1 + b2 + b3| ≥ |b3|
162 |b1 - b2 + b3| ≥ |b3|
163 |b1 + b2 - b3| ≥ |b3|
164 |b1 - b2 - b3| ≥ |b3|
165 """
166 pbc = pbc2pbc(pbc)
167 dim = pbc.sum()
168 if dim <= 1:
169 return True
171 if dim == 2:
172 # reorder cell vectors to [shortest, longest, aperiodic]
173 cell = cell.copy()
174 cell[np.argmin(pbc)] = 0
175 norms = np.linalg.norm(cell, axis=1)
176 cell = cell[np.argsort(norms)[[1, 2, 0]]]
178 A = [[0, 1, 0],
179 [1, -1, 0],
180 [1, 1, 0]]
181 lhs = np.linalg.norm(A @ cell, axis=1)
182 norms = np.linalg.norm(cell, axis=1)
183 rhs = norms[[0, 1, 1]]
184 else:
185 A = [[0, 1, 0],
186 [0, 0, 1],
187 [1, 1, 0],
188 [1, 0, 1],
189 [0, 1, 1],
190 [1, -1, 0],
191 [1, 0, -1],
192 [0, 1, -1],
193 [1, 1, 1],
194 [1, -1, 1],
195 [1, 1, -1],
196 [1, -1, -1]]
197 lhs = np.linalg.norm(A @ cell, axis=1)
198 norms = np.linalg.norm(cell, axis=1)
199 rhs = norms[[0, 1, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2]]
200 return (lhs >= rhs - TOL).all()
203def minkowski_reduce(cell, pbc=True):
204 """Calculate a Minkowski-reduced lattice basis. The reduced basis
205 has the shortest possible vector lengths and has
206 norm(a) <= norm(b) <= norm(c).
208 Implements the method described in:
210 Low-dimensional Lattice Basis Reduction Revisited
211 Nguyen, Phong Q. and Stehlé, Damien,
212 ACM Trans. Algorithms 5(4) 46:1--46:48, 2009
213 :doi:`10.1145/1597036.1597050`
215 Parameters
216 ----------
218 cell: array
219 The lattice basis to reduce (in row-vector format).
220 pbc: array, optional
221 The periodic boundary conditions of the cell (Default `True`).
222 If `pbc` is provided, only periodic cell vectors are reduced.
224 Returns
225 -------
227 rcell: array
228 The reduced lattice basis.
229 op: array
230 The unimodular matrix transformation (rcell = op @ cell).
231 """
232 cell = Cell(cell)
233 pbc = pbc2pbc(pbc)
234 dim = pbc.sum()
235 op = np.eye(3, dtype=int)
236 if is_minkowski_reduced(cell, pbc):
237 return cell, op
239 if dim == 2:
240 # permute cell so that first two vectors are the periodic ones
241 perm = np.argsort(pbc, kind='merge')[::-1] # stable sort
242 pcell = cell[perm][:, perm]
244 # perform gauss reduction
245 norms = np.linalg.norm(pcell, axis=1)
246 norms[2] = float("inf")
247 indices = np.argsort(norms)
248 op = op[indices]
249 hu, hv = reduction_gauss(pcell, op[0], op[1])
250 op[0] = hu
251 op[1] = hv
253 # undo above permutation
254 invperm = np.argsort(perm)
255 op = op[invperm][:, invperm]
257 # maintain cell handedness
258 index = np.argmin(pbc)
259 normal = np.cross(cell[index - 2], cell[index - 1])
260 normal /= np.linalg.norm(normal)
262 _cell = cell.copy()
263 _cell[index] = normal
264 _rcell = op @ cell
265 _rcell[index] = normal
266 if _cell.handedness != Cell(_rcell).handedness:
267 op[index - 1] *= -1
269 elif dim == 3:
270 _, op = reduction_full(cell)
271 # maintain cell handedness
272 if cell.handedness != Cell(op @ cell).handedness:
273 op = -op
275 norms1 = np.sort(np.linalg.norm(cell, axis=1))
276 norms2 = np.sort(np.linalg.norm(op @ cell, axis=1))
277 if (norms2 > norms1 + TOL).any():
278 raise RuntimeError("Minkowski reduction failed")
279 return op @ cell, op