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

1# fmt: off 

2 

3import numpy as np 

4 

5 

6def cellvector_products(cell): 

7 cell = _pad_nonpbc(cell) 

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

9 g0[0] = cell[0] @ cell[0] 

10 g0[1] = cell[1] @ cell[1] 

11 g0[2] = cell[2] @ cell[2] 

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

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

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

15 return g0 

16 

17 

18def _pad_nonpbc(cell): 

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

20 # perpendicular to the periodic ones. 

21 maxlen = max(cell.lengths()) 

22 mask = cell.any(1) 

23 cell = cell.complete() 

24 cell[~mask] *= 2 * maxlen 

25 return cell 

26 

27 

28def niggli_reduce_cell(cell, epsfactor=None): 

29 from ase.cell import Cell 

30 

31 cell = Cell.new(cell) 

32 npbc = cell.rank 

33 

34 if epsfactor is None: 

35 epsfactor = 1e-5 

36 

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

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

39 eps = epsfactor * vol_normalization 

40 

41 g0 = cellvector_products(cell) 

42 g, C = _niggli_reduce(g0, eps) 

43 

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

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

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

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

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

49 

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

51 # We re-zero their lengths before returning: 

52 abc[npbc:] = 0.0 

53 

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

55 

56 newcell[npbc:] = 0.0 

57 return newcell, C 

58 

59 

60def lmn_to_ijk(lmn): 

61 if lmn.prod() == 1: 

62 ijk = lmn.copy() 

63 for idx in range(3): 

64 if ijk[idx] == 0: 

65 ijk[idx] = 1 

66 else: 

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

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

69 r = None 

70 for idx in range(3): 

71 if lmn[idx] == 1: 

72 ijk[idx] = -1 

73 elif lmn[idx] == 0: 

74 r = idx 

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

76 ijk[r] = -1 

77 return ijk 

78 

79 

80def _niggli_reduce(g0, eps): 

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

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

83 

84 C = I3.copy() 

85 D = I6.copy() 

86 

87 g = D @ g0 

88 

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

90 return x < y - eps 

91 

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

93 return lt(y, x, eps) 

94 

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

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

97 

98 for _ in range(10000): 

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

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

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

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

103 g = D @ g0 

104 continue 

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

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

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

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

109 g = D @ g0 

110 continue 

111 

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

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

114 

115 ijk = lmn_to_ijk(lmn) 

116 

117 C *= ijk[np.newaxis] 

118 

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

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

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

122 g = D @ g0 

123 

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

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

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

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

128 

129 A = I3.copy() 

130 A[1, 2] = -s 

131 C = C @ A 

132 

133 B = I6.copy() 

134 B[2, 1] = 1 

135 B[2, 3] = -s 

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

137 B[4, 5] = -s 

138 D = B @ D 

139 g = D @ g0 

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

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

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

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

144 

145 A = I3.copy() 

146 A[0, 2] = -s 

147 C = C @ A 

148 

149 B = I6.copy() 

150 B[2, 0] = 1 

151 B[2, 4] = -s 

152 B[3, 5] = -s 

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

154 D = B @ D 

155 g = D @ g0 

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

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

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

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

160 

161 A = I3.copy() 

162 A[0, 1] = -s 

163 C = C @ A 

164 

165 B = I6.copy() 

166 B[1, 0] = 1 

167 B[1, 5] = -s 

168 B[3, 4] = -s 

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

170 D = B @ D 

171 g = D @ g0 

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

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

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

175 A = I3.copy() 

176 A[:, 2] = 1 

177 C = C @ A 

178 

179 B = I6.copy() 

180 B[2, :] = 1 

181 B[3, 1] = 2 

182 B[3, 5] = 1 

183 B[4, 0] = 2 

184 B[4, 5] = 1 

185 D = B @ D 

186 g = D @ g0 

187 else: 

188 break 

189 else: 

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

191 'g={}\n' 

192 'operation={}' 

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

194 

195 return g, C