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.lattice import bravais_lattices, UnconventionalLattice, bravais_names

4from ase.cell import Cell

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

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

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.

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

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

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}

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.

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)

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

172def find_niggli_ops(latcls, length_grid, angle_grid):

173 niggli_ops = {}

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

176 cell = lat.tocell()

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

186 if op_key in niggli_ops:

187 niggli_ops[op_key] += 1

188 else:

189 niggli_ops[op_key] = 1

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

197 return niggli_ops

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

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

218def generate_niggli_op_table(lattices=None,

219 length_grid=None,

220 angle_grid=None):

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)

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

234 import pprint

235 print(pprint.pformat(niggli_op_table))

236 return niggli_op_table

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 .