r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import itertools

2import numpy as np

3from ase.utils import pbc2pbc

4from ase.cell import Cell

7TOL = 1E-12

8MAX_IT = 100000 # in practice this is not exceeded

11class CycleChecker:

13 def __init__(self, d):

14 assert d in [2, 3]

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

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

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)

24 # flatten array for simplicity

25 H = H.ravel()

27 # check if site exists

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

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

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

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

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]

61def closest_vector(t0, u, v):

62 t = t0

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

64 rs, cs = relevant_vectors_2D(u, v)

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

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

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)

88 for it in range(MAX_IT):

89 # Sort vectors by norm

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

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

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)

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)

108 # Update basis

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

110 R = H @ B

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

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

114 return R, H

119def is_minkowski_reduced(cell, pbc=True):

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

122 Parameters:

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.

130 Returns:

132 is_reduced: bool

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

134 """

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.

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|

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

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

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

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]]]

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()

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).

203 Implements the method described in:

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 :doi:`10.1145/1597036.1597050`

210 Parameters:

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.

218 Returns:

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

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]

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

246 # undo above permutation

247 invperm = np.argsort(perm)

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

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)

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

262 elif dim == 3:

263 _, op = reduction_full(cell)

264 # maintain cell handedness

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

266 op = -op

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