Coverage for /builds/kinetik161/ase/ase/geometry/minkowski_reduction.py: 96.95%

1import itertools 

2 

3import numpy as np 

4 

5from ase.cell import Cell 

6from ase.utils import pbc2pbc 

7 

8TOL = 1E-12 

9MAX_IT = 100000 # in practice this is not exceeded 

10 

11 

12class CycleChecker: 

13 

14 def __init__(self, d): 

15 assert d in [2, 3] 

16 

17 # worst case is the hexagonal cell in 2D and the fcc cell in 3D 

18 n = {2: 6, 3: 12}[d] 

19 

20 # max cycle length is total number of primtive cell descriptions 

21 max_cycle_length = np.prod([n - i for i in range(d)]) * np.prod(d) 

22 self.visited = np.zeros((max_cycle_length, 3 * d), dtype=int) 

23 

24 def add_site(self, H): 

25 # flatten array for simplicity 

26 H = H.ravel() 

27 

28 # check if site exists 

29 found = (self.visited == H).all(axis=1).any() 

30 

31 # shift all visited sites down and place current site at the top 

32 self.visited = np.roll(self.visited, 1, axis=0) 

33 self.visited[0] = H 

34 return found 

35 

36 

37def reduction_gauss(B, hu, hv): 

38 """Calculate a Gauss-reduced lattice basis (2D reduction).""" 

39 cycle_checker = CycleChecker(d=2) 

40 u = hu @ B 

41 v = hv @ B 

42 

43 for it in range(MAX_IT): 

44 x = int(round(np.dot(u, v) / np.dot(u, u))) 

45 hu, hv = hv - x * hu, hu 

46 u = hu @ B 

47 v = hv @ B 

48 site = np.array([hu, hv]) 

49 if np.dot(u, u) >= np.dot(v, v) or cycle_checker.add_site(site): 

50 return hv, hu 

51 

52 raise RuntimeError(f"Gaussian basis not found after {MAX_IT} iterations") 

53 

54 

55def relevant_vectors_2D(u, v): 

56 cs = np.array([e for e in itertools.product([-1, 0, 1], repeat=2)]) 

57 vs = cs @ [u, v] 

58 indices = np.argsort(np.linalg.norm(vs, axis=1))[:7] 

59 return vs[indices], cs[indices] 

60 

61 

62def closest_vector(t0, u, v): 

63 t = t0 

64 a = np.zeros(2, dtype=int) 

65 rs, cs = relevant_vectors_2D(u, v) 

66 

67 dprev = float("inf") 

68 for it in range(MAX_IT): 

69 ds = np.linalg.norm(rs + t, axis=1) 

70 index = np.argmin(ds) 

71 if index == 0 or ds[index] >= dprev: 

72 return a 

73 

74 dprev = ds[index] 

75 r = rs[index] 

76 kopt = int(round(-np.dot(t, r) / np.dot(r, r))) 

77 a += kopt * cs[index] 

78 t = t0 + a[0] * u + a[1] * v 

79 

80 raise RuntimeError(f"Closest vector not found after {MAX_IT} iterations") 

81 

82 

83def reduction_full(B): 

84 """Calculate a Minkowski-reduced lattice basis (3D reduction).""" 

85 cycle_checker = CycleChecker(d=3) 

86 H = np.eye(3, dtype=int) 

87 norms = np.linalg.norm(B, axis=1) 

88 

89 for it in range(MAX_IT): 

90 # Sort vectors by norm 

91 H = H[np.argsort(norms, kind='merge')] 

92 

93 # Gauss-reduce smallest two vectors 

94 hw = H[2] 

95 hu, hv = reduction_gauss(B, H[0], H[1]) 

96 H = np.array([hu, hv, hw]) 

97 R = H @ B 

98 

99 # Orthogonalize vectors using Gram-Schmidt 

100 u, v, _ = R 

101 X = u / np.linalg.norm(u) 

102 Y = v - X * np.dot(v, X) 

103 Y /= np.linalg.norm(Y) 

104 

105 # Find closest vector to last element of R 

106 pu, pv, pw = R @ np.array([X, Y]).T 

107 nb = closest_vector(pw, pu, pv) 

108 

109 # Update basis 

110 H[2] = [nb[0], nb[1], 1] @ H 

111 R = H @ B 

112 

113 norms = np.linalg.norm(R, axis=1) 

114 if norms[2] >= norms[1] or cycle_checker.add_site(H): 

115 return R, H 

116 

117 raise RuntimeError(f"Reduced basis not found after {MAX_IT} iterations") 

118 

119 

120def is_minkowski_reduced(cell, pbc=True): 

121 """Tests if a cell is Minkowski-reduced. 

122 

123 Parameters: 

124 

125 cell: array 

126 The lattice basis to test (in row-vector format). 

127 pbc: array, optional 

128 The periodic boundary conditions of the cell (Default `True`). 

129 If `pbc` is provided, only periodic cell vectors are tested. 

130 

131 Returns: 

132 

133 is_reduced: bool 

134 True if cell is Minkowski-reduced, False otherwise. 

135 """ 

136 

137 """These conditions are due to Minkowski, but a nice description in English 

138 can be found in the thesis of Carine Jaber: "Algorithmic approaches to 

139 Siegel's fundamental domain", https://www.theses.fr/2017UBFCK006.pdf 

140 This is also good background reading for Minkowski reduction. 

141 

142 0D and 1D cells are trivially reduced. For 2D cells, the conditions which 

143 an already-reduced basis fulfil are: 

144 |b1| ≤ |b2| 

145 |b2| ≤ |b1 - b2| 

146 |b2| ≤ |b1 + b2| 

147 

148 For 3D cells, the conditions which an already-reduced basis fulfil are: 

149 |b1| ≤ |b2| ≤ |b3| 

150 

151 |b1 + b2| ≥ |b2| 

152 |b1 + b3| ≥ |b3| 

153 |b2 + b3| ≥ |b3| 

154 |b1 - b2| ≥ |b2| 

155 |b1 - b3| ≥ |b3| 

156 |b2 - b3| ≥ |b3| 

157 |b1 + b2 + b3| ≥ |b3| 

158 |b1 - b2 + b3| ≥ |b3| 

159 |b1 + b2 - b3| ≥ |b3| 

160 |b1 - b2 - b3| ≥ |b3| 

161 """ 

162 pbc = pbc2pbc(pbc) 

163 dim = pbc.sum() 

164 if dim <= 1: 

165 return True 

166 

167 if dim == 2: 

168 # reorder cell vectors to [shortest, longest, aperiodic] 

169 cell = cell.copy() 

170 cell[np.argmin(pbc)] = 0 

171 norms = np.linalg.norm(cell, axis=1) 

172 cell = cell[np.argsort(norms)[[1, 2, 0]]] 

173 

174 A = [[0, 1, 0], 

175 [1, -1, 0], 

176 [1, 1, 0]] 

177 lhs = np.linalg.norm(A @ cell, axis=1) 

178 norms = np.linalg.norm(cell, axis=1) 

179 rhs = norms[[0, 1, 1]] 

180 else: 

181 A = [[0, 1, 0], 

182 [0, 0, 1], 

183 [1, 1, 0], 

184 [1, 0, 1], 

185 [0, 1, 1], 

186 [1, -1, 0], 

187 [1, 0, -1], 

188 [0, 1, -1], 

189 [1, 1, 1], 

190 [1, -1, 1], 

191 [1, 1, -1], 

192 [1, -1, -1]] 

193 lhs = np.linalg.norm(A @ cell, axis=1) 

194 norms = np.linalg.norm(cell, axis=1) 

195 rhs = norms[[0, 1, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2]] 

196 return (lhs >= rhs - TOL).all() 

197 

198 

199def minkowski_reduce(cell, pbc=True): 

200 """Calculate a Minkowski-reduced lattice basis. The reduced basis 

201 has the shortest possible vector lengths and has 

202 norm(a) <= norm(b) <= norm(c). 

203 

204 Implements the method described in: 

205 

206 Low-dimensional Lattice Basis Reduction Revisited 

207 Nguyen, Phong Q. and Stehlé, Damien, 

208 ACM Trans. Algorithms 5(4) 46:1--46:48, 2009 

209 :doi:`10.1145/1597036.1597050` 

210 

211 Parameters: 

212 

213 cell: array 

214 The lattice basis to reduce (in row-vector format). 

215 pbc: array, optional 

216 The periodic boundary conditions of the cell (Default `True`). 

217 If `pbc` is provided, only periodic cell vectors are reduced. 

218 

219 Returns: 

220 

221 rcell: array 

222 The reduced lattice basis. 

223 op: array 

224 The unimodular matrix transformation (rcell = op @ cell). 

225 """ 

226 cell = Cell(cell) 

227 pbc = pbc2pbc(pbc) 

228 dim = pbc.sum() 

229 op = np.eye(3, dtype=int) 

230 if is_minkowski_reduced(cell, pbc): 

231 return cell, op 

232 

233 if dim == 2: 

234 # permute cell so that first two vectors are the periodic ones 

235 perm = np.argsort(pbc, kind='merge')[::-1] # stable sort 

236 pcell = cell[perm][:, perm] 

237 

238 # perform gauss reduction 

239 norms = np.linalg.norm(pcell, axis=1) 

240 norms[2] = float("inf") 

241 indices = np.argsort(norms) 

242 op = op[indices] 

243 hu, hv = reduction_gauss(pcell, op[0], op[1]) 

244 op[0] = hu 

245 op[1] = hv 

246 

247 # undo above permutation 

248 invperm = np.argsort(perm) 

249 op = op[invperm][:, invperm] 

250 

251 # maintain cell handedness 

252 index = np.argmin(pbc) 

253 normal = np.cross(cell[index - 2], cell[index - 1]) 

254 normal /= np.linalg.norm(normal) 

255 

256 _cell = cell.copy() 

257 _cell[index] = normal 

258 _rcell = op @ cell 

259 _rcell[index] = normal 

260 if _cell.handedness != Cell(_rcell).handedness: 

261 op[index - 1] *= -1 

262 

263 elif dim == 3: 

264 _, op = reduction_full(cell) 

265 # maintain cell handedness 

266 if cell.handedness != Cell(op @ cell).handedness: 

267 op = -op 

268 

269 norms1 = np.sort(np.linalg.norm(cell, axis=1)) 

270 norms2 = np.sort(np.linalg.norm(op @ cell, axis=1)) 

271 if (norms2 > norms1 + TOL).any(): 

272 raise RuntimeError("Minkowski reduction failed") 

273 return op @ cell, op