r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np

4def cellvector_products(cell):

6 g0 = np.empty(6, dtype=float)

7 g0[0] = cell[0] @ cell[0]

8 g0[1] = cell[1] @ cell[1]

9 g0[2] = cell[2] @ cell[2]

10 g0[3] = 2 * (cell[1] @ cell[2])

11 g0[4] = 2 * (cell[2] @ cell[0])

12 g0[5] = 2 * (cell[0] @ cell[1])

13 return g0

17 # Add "infinitely long" lattice vectors for non-periodic directions,

18 # perpendicular to the periodic ones.

19 maxlen = max(cell.lengths())

21 cell = cell.complete()

22 cell[~mask] *= 2 * maxlen

23 return cell

26def niggli_reduce_cell(cell, epsfactor=None):

27 from ase.cell import Cell

29 cell = Cell.new(cell)

30 npbc = cell.rank

32 if epsfactor is None:

33 epsfactor = 1e-5

35 vol_normalization_exponent = 1 if npbc == 0 else 1 / npbc

36 vol_normalization = cell.complete().volume**vol_normalization_exponent

37 eps = epsfactor * vol_normalization

39 g0 = cellvector_products(cell)

40 g, C = _niggli_reduce(g0, eps)

42 abc = np.sqrt(g[:3])

43 # Prevent division by zero e.g. for cell==zeros((3, 3)):

44 abcprod = max(abc.prod(), 1e-100)

45 cosangles = abc * g[3:] / (2 * abcprod)

46 angles = 180 * np.arccos(cosangles) / np.pi

48 # Non-periodic directions have artificial infinitely long lattice vectors.

49 # We re-zero their lengths before returning:

50 abc[npbc:] = 0.0

52 newcell = Cell.fromcellpar(np.concatenate([abc, angles]))

54 newcell[npbc:] = 0.0

55 return newcell, C

58def lmn_to_ijk(lmn):

59 if lmn.prod() == 1:

60 ijk = lmn.copy()

61 for idx in range(3):

62 if ijk[idx] == 0:

63 ijk[idx] = 1

64 else:

65 ijk = np.ones(3, dtype=int)

66 if np.any(lmn != -1):

67 r = None

68 for idx in range(3):

69 if lmn[idx] == 1:

70 ijk[idx] = -1

71 elif lmn[idx] == 0:

72 r = idx

73 if ijk.prod() == -1:

74 ijk[r] = -1

75 return ijk

78def _niggli_reduce(g0, eps):

79 I3 = np.eye(3, dtype=int)

80 I6 = np.eye(6, dtype=int)

82 C = I3.copy()

83 D = I6.copy()

85 g = D @ g0

87 def lt(x, y, eps=eps):

88 return x < y - eps

90 def gt(x, y, eps=eps):

91 return lt(y, x, eps)

93 def eq(x, y, eps=eps):

94 return not (lt(x, y, eps) or gt(x, y, eps))

96 for _ in range(10000):

97 if (gt(g[0], g[1])

98 or (eq(g[0], g[1]) and gt(abs(g[3]), abs(g[4])))):

99 C = C @ (-I3[[1, 0, 2]])

100 D = I6[[1, 0, 2, 4, 3, 5]] @ D

101 g = D @ g0

102 continue

103 elif (gt(g[1], g[2])

104 or (eq(g[1], g[2]) and gt(abs(g[4]), abs(g[5])))):

105 C = C @ (-I3[[0, 2, 1]])

106 D = I6[[0, 2, 1, 3, 5, 4]] @ D

107 g = D @ g0

108 continue

110 lmn = np.array(gt(g[3:], 0, eps=eps / 2), dtype=int)

111 lmn -= np.array(lt(g[3:], 0, eps=eps / 2), dtype=int)

113 ijk = lmn_to_ijk(lmn)

115 C *= ijk[np.newaxis]

117 D[3] *= ijk[1] * ijk[2]

118 D[4] *= ijk[0] * ijk[2]

119 D[5] *= ijk[0] * ijk[1]

120 g = D @ g0

122 if (gt(abs(g[3]), g[1])

123 or (eq(g[3], g[1]) and lt(2 * g[4], g[5]))

124 or (eq(g[3], -g[1]) and lt(g[5], 0))):

125 s = int(np.sign(g[3]))

127 A = I3.copy()

128 A[1, 2] = -s

129 C = C @ A

131 B = I6.copy()

132 B[2, 1] = 1

133 B[2, 3] = -s

134 B[3, 1] = -2 * s

135 B[4, 5] = -s

136 D = B @ D

137 g = D @ g0

138 elif (gt(abs(g[4]), g[0])

139 or (eq(g[4], g[0]) and lt(2 * g[3], g[5]))

140 or (eq(g[4], -g[0]) and lt(g[5], 0))):

141 s = int(np.sign(g[4]))

143 A = I3.copy()

144 A[0, 2] = -s

145 C = C @ A

147 B = I6.copy()

148 B[2, 0] = 1

149 B[2, 4] = -s

150 B[3, 5] = -s

151 B[4, 0] = -2 * s

152 D = B @ D

153 g = D @ g0

154 elif (gt(abs(g[5]), g[0])

155 or (eq(g[5], g[0]) and lt(2 * g[3], g[4]))

156 or (eq(g[5], -g[0]) and lt(g[4], 0))):

157 s = int(np.sign(g[5]))

159 A = I3.copy()

160 A[0, 1] = -s

161 C = C @ A

163 B = I6.copy()

164 B[1, 0] = 1

165 B[1, 5] = -s

166 B[3, 4] = -s

167 B[5, 0] = -2 * s

168 D = B @ D

169 g = D @ g0

170 elif (lt(g[[0, 1, 3, 4, 5]].sum(), 0)

171 or (eq(g[[0, 1, 3, 4, 5]].sum(), 0)

172 and gt(2 * (g[0] + g[4]) + g[5], 0))):

173 A = I3.copy()

174 A[:, 2] = 1

175 C = C @ A

177 B = I6.copy()

178 B[2, :] = 1

179 B[3, 1] = 2

180 B[3, 5] = 1

181 B[4, 0] = 2

182 B[4, 5] = 1

183 D = B @ D

184 g = D @ g0

185 else:

186 break

187 else:

188 raise RuntimeError('Niggli reduction not done in 10000 steps!\n'

189 'g={}\n'

190 'operation={}'

191 .format(g.tolist(), C.tolist()))

193 return g, C