1import itertools 

2import numpy as np 

3from ase.utils import pbc2pbc 

4from ase.cell import Cell 

5 

6 

7TOL = 1E-12 

8MAX_IT = 100000 # in practice this is not exceeded 

9 

10 

11class CycleChecker: 

12 

13 def __init__(self, d): 

14 assert d in [2, 3] 

15 

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

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

18 

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

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

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

22 

23 def add_site(self, H): 

24 # flatten array for simplicity 

25 H = H.ravel() 

26 

27 # check if site exists 

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

29 

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

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

32 self.visited[0] = H 

33 return found 

34 

35 

36def reduction_gauss(B, hu, hv): 

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

38 cycle_checker = CycleChecker(d=2) 

39 u = hu @ B 

40 v = hv @ B 

41 

42 for it in range(MAX_IT): 

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

44 hu, hv = hv - x * hu, hu 

45 u = hu @ B 

46 v = hv @ B 

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

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

49 return hv, hu 

50 

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

52 

53 

54def relevant_vectors_2D(u, v): 

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

56 vs = cs @ [u, v] 

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

58 return vs[indices], cs[indices] 

59 

60 

61def closest_vector(t0, u, v): 

62 t = t0 

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

64 rs, cs = relevant_vectors_2D(u, v) 

65 

66 dprev = float("inf") 

67 for it in range(MAX_IT): 

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

69 index = np.argmin(ds) 

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

71 return a 

72 

73 dprev = ds[index] 

74 r = rs[index] 

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

76 a += kopt * cs[index] 

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

78 

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

80 

81 

82def reduction_full(B): 

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

84 cycle_checker = CycleChecker(d=3) 

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

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

87 

88 for it in range(MAX_IT): 

89 # Sort vectors by norm 

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

91 

92 # Gauss-reduce smallest two vectors 

93 hw = H[2] 

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

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

96 R = H @ B 

97 

98 # Orthogonalize vectors using Gram-Schmidt 

99 u, v, _ = R 

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

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

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

103 

104 # Find closest vector to last element of R 

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

106 nb = closest_vector(pw, pu, pv) 

107 

108 # Update basis 

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

110 R = H @ B 

111 

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

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

114 return R, H 

115 

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

117 

118 

119def is_minkowski_reduced(cell, pbc=True): 

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

121 

122 Parameters: 

123 

124 cell: array 

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

126 pbc: array, optional 

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

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

129 

130 Returns: 

131 

132 is_reduced: bool 

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

134 """ 

135 

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

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

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

139 This is also good background reading for Minkowski reduction. 

140 

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

142 an already-reduced basis fulfil are: 

143 |b1| ≤ |b2| 

144 |b2| ≤ |b1 - b2| 

145 |b2| ≤ |b1 + b2| 

146 

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

148 |b1| ≤ |b2| ≤ |b3| 

149 

150 |b1 + b2| ≥ |b2| 

151 |b1 + b3| ≥ |b3| 

152 |b2 + b3| ≥ |b3| 

153 |b1 - b2| ≥ |b2| 

154 |b1 - b3| ≥ |b3| 

155 |b2 - b3| ≥ |b3| 

156 |b1 + b2 + b3| ≥ |b3| 

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

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

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

160 """ 

161 pbc = pbc2pbc(pbc) 

162 dim = pbc.sum() 

163 if dim <= 1: 

164 return True 

165 

166 if dim == 2: 

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

168 cell = cell.copy() 

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

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

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

172 

173 A = [[0, 1, 0], 

174 [1, -1, 0], 

175 [1, 1, 0]] 

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

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

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

179 else: 

180 A = [[0, 1, 0], 

181 [0, 0, 1], 

182 [1, 1, 0], 

183 [1, 0, 1], 

184 [0, 1, 1], 

185 [1, -1, 0], 

186 [1, 0, -1], 

187 [0, 1, -1], 

188 [1, 1, 1], 

189 [1, -1, 1], 

190 [1, 1, -1], 

191 [1, -1, -1]] 

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

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

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

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

196 

197 

198def minkowski_reduce(cell, pbc=True): 

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

200 has the shortest possible vector lengths and has 

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

202 

203 Implements the method described in: 

204 

205 Low-dimensional Lattice Basis Reduction Revisited 

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

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

208 https://doi.org/10.1145/1597036.1597050 

209 

210 Parameters: 

211 

212 cell: array 

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

214 pbc: array, optional 

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

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

217 

218 Returns: 

219 

220 rcell: array 

221 The reduced lattice basis. 

222 op: array 

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

224 """ 

225 cell = Cell(cell) 

226 pbc = pbc2pbc(pbc) 

227 dim = pbc.sum() 

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

229 if is_minkowski_reduced(cell, pbc): 

230 return cell, op 

231 

232 if dim == 2: 

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

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

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

236 

237 # perform gauss reduction 

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

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

240 indices = np.argsort(norms) 

241 op = op[indices] 

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

243 op[0] = hu 

244 op[1] = hv 

245 

246 # undo above permutation 

247 invperm = np.argsort(perm) 

248 op = op[invperm][:, invperm] 

249 

250 # maintain cell handedness 

251 index = np.argmin(pbc) 

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

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

254 

255 _cell = cell.copy() 

256 _cell[index] = normal 

257 _rcell = op @ cell 

258 _rcell[index] = normal 

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

260 op[index - 1] *= -1 

261 

262 elif dim == 3: 

263 _, op = reduction_full(cell) 

264 # maintain cell handedness 

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

266 op = -op 

267 

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

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

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

271 raise RuntimeError("Minkowski reduction failed") 

272 return op @ cell, op