1import itertools 

2import numpy as np 

3from ase.lattice import bravais_lattices, UnconventionalLattice, bravais_names 

4from ase.cell import Cell 

5 

6"""This module implements a crude method to recognize most Bravais lattices. 

7 

8Suppose we use the ase.lattice module to generate many 

9lattices of some particular type, say, BCT(a, c), and then we 

10Niggli-reduce all of them. The Niggli-reduced forms are not 

11immediately recognizable, but we know the mapping from each reduced 

12form back to the original form. As it turns out, there are apparently 

135 such mappings (the proof is left as an exercise for the reader). 

14 

15Hence, presumably, every BCT lattice (whether generated by BCT(a, c) 

16or in some other form) Niggli-reduces to a form which, through the 

17inverse of one of those five operations, is mapped back to the 

18recognizable one. Knowing all five operations (or equivalence 

19classes), we can characterize any BCT lattice. Same goes for the 

20other lattices of sufficiently low dimension. 

21 

22For MCL, MCLC, and TRI, we may not recognize all forms correctly, 

23but we aspire that this will work for all common inputs.""" 

24 

25 

26niggli_op_table = { # Generated by generate_niggli_op_table() 

27 'LINE': [(1, 0, 0, 0, 1, 0, 0, 0, 1)], 

28 'SQR': [(1, 0, 0, 0, 1, 0, 0, 0, 1)], 

29 'RECT': [(1, 0, 0, 0, 1, 0, 0, 0, 1)], 

30 'CRECT': [(-1, 1, 0, 1, 0, 0, 0, 0, -1), 

31 (1, 0, 0, 0, -1, 0, 0, 0, -1)], 

32 'HEX2D': [(1, 0, 0, 0, 1, 0, 0, 0, 1)], 

33 'OBL': [(-1, 1, 0, 1, 0, 0, 0, 0, -1), 

34 (1, -1, 0, 0, 1, 0, 0, 0, 1), 

35 (1, 0, 0, 0, -1, 0, 0, 0, -1), 

36 (-1, -1, 0, 1, 0, 0, 0, 0, 1), 

37 (1, 1, 0, 0, -1, 0, 0, 0, -1)], 

38 'BCC': [(1, 0, 0, 0, 1, 0, 0, 0, 1)], 

39 'BCT': [(1, 0, 0, 0, 1, 0, 0, 0, 1), 

40 (0, 1, 0, 0, 0, 1, 1, 0, 0), 

41 (0, 1, 0, 1, 0, 0, 1, 1, -1), 

42 (-1, 0, 1, 0, 1, 0, -1, 1, 0), 

43 (1, 1, 0, 1, 0, 0, 0, 0, -1)], 

44 'CUB': [(1, 0, 0, 0, 1, 0, 0, 0, 1)], 

45 'FCC': [(1, 0, 0, 0, 1, 0, 0, 0, 1)], 

46 'HEX': [(1, 0, 0, 0, 1, 0, 0, 0, 1), (0, 1, 0, 0, 0, 1, 1, 0, 0)], 

47 'ORC': [(1, 0, 0, 0, 1, 0, 0, 0, 1)], 

48 'ORCC': [(1, 0, 0, 0, 1, 0, 0, 0, 1), 

49 (1, 0, -1, 1, 0, 0, 0, -1, 0), 

50 (-1, 1, 0, -1, 0, 0, 0, 0, 1), 

51 (0, 1, 0, 0, 0, 1, 1, 0, 0), 

52 (0, -1, 1, 0, -1, 0, 1, 0, 0)], 

53 'ORCF': [(0, -1, 0, 0, 1, -1, 1, 0, 0), (-1, 0, 0, 1, 0, 1, 1, 1, 0)], 

54 'ORCI': [(0, 0, -1, 0, -1, 0, -1, 0, 0), 

55 (0, 0, 1, -1, 0, 0, -1, -1, 0), 

56 (0, 1, 0, 1, 0, 0, 1, 1, -1), 

57 (0, -1, 0, 1, 0, -1, 1, -1, 0)], 

58 'RHL': [(0, -1, 0, 1, 1, 1, -1, 0, 0), 

59 (1, 0, 0, 0, 1, 0, 0, 0, 1), 

60 (1, -1, 0, 1, 0, -1, 1, 0, 0)], 

61 'TET': [(1, 0, 0, 0, 1, 0, 0, 0, 1), (0, 1, 0, 0, 0, 1, 1, 0, 0)], 

62 'MCL': [(0, 0, 1, -1, -1, 0, 1, 0, 0), 

63 (-1, 0, 0, 0, 1, 0, 0, 0, -1), 

64 (0, 0, -1, 1, 1, 0, 0, -1, 0), 

65 (0, -1, 0, 1, 0, 1, -1, 0, 0), 

66 (0, 1, 0, -1, 0, -1, 0, 0, 1), 

67 (-1, 0, 0, 0, 1, 1, 0, 0, -1), 

68 (0, 1, 0, 1, 0, -1, -1, 0, 0), 

69 (0, 0, 1, 1, -1, 0, 0, 1, 0), 

70 (0, 1, 0, -1, 0, 0, 0, 0, 1), 

71 (0, 0, -1, -1, 1, 0, 1, 0, 0), 

72 (1, 0, 0, 0, 1, -1, 0, 0, 1), 

73 (0, -1, 0, -1, 0, 1, 0, 0, -1), 

74 (-1, 0, 0, 0, -1, 1, 0, 1, 0), 

75 (1, 0, 0, 0, -1, -1, 0, 1, 0), 

76 (0, 0, -1, 1, 0, 0, 0, -1, 0)], 

77 'MCLC': [(1, 1, 1, 1, 0, 1, 0, 0, -1), 

78 (1, 1, 1, 1, 1, 0, -1, 0, 0), 

79 (1, -1, 1, -1, 0, 1, 0, 0, -1), 

80 (-1, 1, 0, 1, 0, 0, 0, 0, -1), 

81 (1, 0, 0, 0, 1, 0, 0, 0, 1), 

82 (-1, 0, -1, 1, -1, -1, 0, 0, 1), 

83 (1, -1, -1, 1, -1, 0, -1, 0, 0), 

84 (-1, -1, 0, -1, 0, -1, 1, 0, 0), 

85 (1, 0, 1, 1, 0, 0, 0, 1, 0), 

86 (-1, 1, 0, -1, 0, 1, 1, 0, 0), 

87 (0, -1, 1, -1, 0, 1, 0, 0, -1), 

88 (-1, -1, 0, -1, 0, 0, 0, 0, -1), 

89 (-1, -1, 1, -1, 0, 1, 0, 0, -1), 

90 (1, 0, 0, 0, -1, 1, 0, 0, -1), 

91 (-1, 0, -1, 0, -1, -1, 0, 0, 1), 

92 (1, 0, -1, -1, 1, -1, 0, 0, 1), 

93 (1, -1, 1, 1, -1, 0, 0, 1, 0), 

94 (0, -1, 0, 1, 0, -1, 0, 0, 1), 

95 (-1, 0, 0, 1, 1, 1, 0, 0, -1), 

96 (1, 0, -1, 0, 1, -1, 0, 0, 1), 

97 (-1, 1, 0, 1, 1, -1, 0, -1, 0), 

98 (1, 1, -1, 1, -1, 0, -1, 0, 0), 

99 (-1, -1, -1, -1, -1, 0, 0, 1, 0), 

100 (-1, 1, 1, 1, 0, 1, 0, 0, -1), 

101 (-1, 0, 0, 0, -1, 0, 0, 0, 1), 

102 (-1, -1, 1, 1, -1, 0, 0, 1, 0), 

103 (1, 1, 0, -1, 0, -1, 0, 0, 1)], 

104 'TRI': [(0, -1, 0, -1, 0, 0, 0, 0, -1), 

105 (0, 1, 0, 0, 0, 1, 1, 0, 0), 

106 (0, 0, -1, 0, -1, 0, -1, 1, 0), 

107 (0, 0, 1, 0, 1, 0, -1, 0, 0), 

108 (0, -1, 0, 0, 0, -1, 1, 1, 1), 

109 (0, 1, 0, 0, 0, 1, 1, -1, 0), 

110 (0, 0, -1, 0, -1, 0, -1, 0, 0), 

111 (-1, 1, 0, 0, 0, -1, 0, -1, 0), 

112 (0, 0, 1, 1, -1, 0, 0, 1, 0), 

113 (0, 0, -1, 1, 1, 1, 0, -1, 0), 

114 (-1, 0, 0, 0, 1, 0, 0, -1, -1), 

115 (0, 0, 1, 1, 0, 0, 0, 1, 0), 

116 (0, 0, 1, 0, 1, 0, -1, -1, -1), 

117 (-1, 0, 0, 0, 0, -1, 0, -1, 0), 

118 (0, -1, 0, 0, 0, -1, 1, 0, 0), 

119 (1, 0, 0, 0, 1, 0, 0, 0, 1), 

120 (0, 0, -1, -1, 0, 0, 1, 1, 1), 

121 (0, 0, -1, -1, 0, 0, 0, 1, 0), 

122 (-1, -1, -1, 0, 0, 1, 0, 1, 0)], 

123} 

124 

125# XXX The TRI list was generated by looping over all TRI structures in 

126# the COD (Crystallography Open Database) and seeing what operations 

127# were necessary to map all those to standard form. Hence if the 

128# data does not cover all possible inputs, we could miss something. 

129# 

130# Looping over all possible TRI lattices in general would generate 

131# 100+ operations, we don't want to tabulate that. 

132 

133 

134def lattice_loop(latcls, length_grid, angle_grid): 

135 """Yield all lattices defined by the length and angle grids.""" 

136 param_grids = [] 

137 for varname in latcls.parameters: 

138 # Actually we could choose one parameter, a, to always be 1, 

139 # reducing the dimension of the problem by 1. The lattice 

140 # recognition code should do something like that as well, but 

141 # it doesn't. This could affect the impact of the eps value 

142 # on lattice determination, so we just loop over the whole 

143 # thing in order not to worry. 

144 if latcls.name in ['MCL', 'MCLC']: 

145 special_var = 'c' 

146 else: 

147 special_var = 'a' 

148 if varname == special_var: 

149 values = np.ones(1) 

150 elif varname in 'abc': 

151 values = length_grid 

152 elif varname in ['alpha', 'beta', 'gamma']: 

153 values = angle_grid 

154 else: 

155 raise ValueError(varname) 

156 param_grids.append(values) 

157 

158 for latpars in itertools.product(*param_grids): 

159 kwargs = dict(zip(latcls.parameters, latpars)) 

160 try: 

161 lat = latcls(**kwargs) 

162 except (UnconventionalLattice, AssertionError): 

163 # XXX assertion error can happen because cellpar_to_cell 

164 # makes certain assumptions. Should be investigated. 

165 # {'b': 0.1, 'gamma': 60.0, 'c': 0.1, 'a': 1.0, 

166 # 'alpha': 30.0, 'beta': 30.0} <-- this won't work 

167 pass 

168 else: 

169 yield lat 

170 

171 

172def find_niggli_ops(latcls, length_grid, angle_grid): 

173 niggli_ops = {} 

174 

175 for lat in lattice_loop(latcls, length_grid, angle_grid): 

176 cell = lat.tocell() 

177 

178 try: 

179 rcell, op = cell.niggli_reduce() 

180 except RuntimeError: 

181 print('Niggli reduce did not converge') 

182 continue 

183 assert op.dtype == int 

184 op_key = tuple(op.ravel()) 

185 

186 if op_key in niggli_ops: 

187 niggli_ops[op_key] += 1 

188 else: 

189 niggli_ops[op_key] = 1 

190 

191 rcell_test = Cell(op.T @ cell) 

192 rcellpar_test = rcell_test.cellpar() 

193 rcellpar = rcell.cellpar() 

194 err = np.abs(rcellpar_test - rcellpar).max() 

195 assert err < 1e-7, err 

196 

197 return niggli_ops 

198 

199 

200def find_all_niggli_ops(length_grid, angle_grid, lattices=None): 

201 all_niggli_ops = {} 

202 if lattices is None: 

203 lattices = [name for name in bravais_names 

204 if name not in ['MCL', 'MCLC', 'TRI']] 

205 

206 for latname in lattices: 

207 latcls = bravais_lattices[latname] 

208 print('Working on {}...'.format(latname)) 

209 niggli_ops = find_niggli_ops(latcls, length_grid, angle_grid) 

210 print('Found {} ops for {}'.format(len(niggli_ops), latname)) 

211 for key, count in niggli_ops.items(): 

212 print(' {:>40}: {}'.format(str(np.array(key)), count)) 

213 print() 

214 all_niggli_ops[latname] = niggli_ops 

215 return all_niggli_ops 

216 

217 

218def generate_niggli_op_table(lattices=None, 

219 length_grid=None, 

220 angle_grid=None): 

221 

222 if length_grid is None: 

223 length_grid = np.logspace(-0.5, 1.5, 50).round(3) 

224 if angle_grid is None: 

225 angle_grid = np.linspace(10, 179, 50).round() 

226 all_niggli_ops_and_counts = find_all_niggli_ops(length_grid, angle_grid, 

227 lattices=lattices) 

228 

229 niggli_op_table = {} 

230 for latname, ops in all_niggli_ops_and_counts.items(): 

231 ops = [op for op in ops if np.abs(op).max() < 2] 

232 niggli_op_table[latname] = ops 

233 

234 import pprint 

235 print(pprint.pformat(niggli_op_table)) 

236 return niggli_op_table 

237 

238 

239# For generation of the table, please see the test_bravais_type_engine 

240# unit test. In case there's any trouble, some legacy code can be 

241# found also in 6e2b1c6cae0ae6ee04638a9887821e7b1a1f2f3f .