Coverage for /builds/kinetik161/ase/ase/build/niggli.py: 98.52%

1import numpy as np 

2 

3 

4def cellvector_products(cell): 

5 cell = _pad_nonpbc(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 

14 

15 

16def _pad_nonpbc(cell): 

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

18 # perpendicular to the periodic ones. 

19 maxlen = max(cell.lengths()) 

20 mask = cell.any(1) 

21 cell = cell.complete() 

22 cell[~mask] *= 2 * maxlen 

23 return cell 

24 

25 

26def niggli_reduce_cell(cell, epsfactor=None): 

27 from ase.cell import Cell 

28 

29 cell = Cell.new(cell) 

30 npbc = cell.rank 

31 

32 if epsfactor is None: 

33 epsfactor = 1e-5 

34 

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 

38 

39 g0 = cellvector_products(cell) 

40 g, C = _niggli_reduce(g0, eps) 

41 

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 

47 

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

49 # We re-zero their lengths before returning: 

50 abc[npbc:] = 0.0 

51 

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

53 

54 newcell[npbc:] = 0.0 

55 return newcell, C 

56 

57 

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 

76 

77 

78def _niggli_reduce(g0, eps): 

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

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

81 

82 C = I3.copy() 

83 D = I6.copy() 

84 

85 g = D @ g0 

86 

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

88 return x < y - eps 

89 

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

91 return lt(y, x, eps) 

92 

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

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

95 

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 

109 

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) 

112 

113 ijk = lmn_to_ijk(lmn) 

114 

115 C *= ijk[np.newaxis] 

116 

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 

121 

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

126 

127 A = I3.copy() 

128 A[1, 2] = -s 

129 C = C @ A 

130 

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

142 

143 A = I3.copy() 

144 A[0, 2] = -s 

145 C = C @ A 

146 

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

158 

159 A = I3.copy() 

160 A[0, 1] = -s 

161 C = C @ A 

162 

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 

176 

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

192 

193 return g, C