Coverage for /builds/kinetik161/ase/ase/geometry/bravais_type_engine.py: 89.19%
74 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import itertools
3import numpy as np
5from ase.cell import Cell
6from ase.lattice import UnconventionalLattice, bravais_lattices, bravais_names
8"""This module implements a crude method to recognize most Bravais lattices.
10Suppose we use the ase.lattice module to generate many
11lattices of some particular type, say, BCT(a, c), and then we
12Niggli-reduce all of them. The Niggli-reduced forms are not
13immediately recognizable, but we know the mapping from each reduced
14form back to the original form. As it turns out, there are apparently
155 such mappings (the proof is left as an exercise for the reader).
17Hence, presumably, every BCT lattice (whether generated by BCT(a, c)
18or in some other form) Niggli-reduces to a form which, through the
19inverse of one of those five operations, is mapped back to the
20recognizable one. Knowing all five operations (or equivalence
21classes), we can characterize any BCT lattice. Same goes for the
22other lattices of sufficiently low dimension.
24For MCL, MCLC, and TRI, we may not recognize all forms correctly,
25but we aspire that this will work for all common inputs."""
28niggli_op_table = { # Generated by generate_niggli_op_table()
29 'LINE': [(1, 0, 0, 0, 1, 0, 0, 0, 1)],
30 'SQR': [(1, 0, 0, 0, 1, 0, 0, 0, 1)],
31 'RECT': [(1, 0, 0, 0, 1, 0, 0, 0, 1)],
32 'CRECT': [(-1, 1, 0, 1, 0, 0, 0, 0, -1),
33 (1, 0, 0, 0, -1, 0, 0, 0, -1)],
34 'HEX2D': [(1, 0, 0, 0, 1, 0, 0, 0, 1)],
35 'OBL': [(-1, 1, 0, 1, 0, 0, 0, 0, -1),
36 (1, -1, 0, 0, 1, 0, 0, 0, 1),
37 (1, 0, 0, 0, -1, 0, 0, 0, -1),
38 (-1, -1, 0, 1, 0, 0, 0, 0, 1),
39 (1, 1, 0, 0, -1, 0, 0, 0, -1)],
40 'BCC': [(1, 0, 0, 0, 1, 0, 0, 0, 1)],
41 'BCT': [(1, 0, 0, 0, 1, 0, 0, 0, 1),
42 (0, 1, 0, 0, 0, 1, 1, 0, 0),
43 (0, 1, 0, 1, 0, 0, 1, 1, -1),
44 (-1, 0, 1, 0, 1, 0, -1, 1, 0),
45 (1, 1, 0, 1, 0, 0, 0, 0, -1)],
46 'CUB': [(1, 0, 0, 0, 1, 0, 0, 0, 1)],
47 'FCC': [(1, 0, 0, 0, 1, 0, 0, 0, 1)],
48 'HEX': [(1, 0, 0, 0, 1, 0, 0, 0, 1), (0, 1, 0, 0, 0, 1, 1, 0, 0)],
49 'ORC': [(1, 0, 0, 0, 1, 0, 0, 0, 1)],
50 'ORCC': [(1, 0, 0, 0, 1, 0, 0, 0, 1),
51 (1, 0, -1, 1, 0, 0, 0, -1, 0),
52 (-1, 1, 0, -1, 0, 0, 0, 0, 1),
53 (0, 1, 0, 0, 0, 1, 1, 0, 0),
54 (0, -1, 1, 0, -1, 0, 1, 0, 0)],
55 'ORCF': [(0, -1, 0, 0, 1, -1, 1, 0, 0), (-1, 0, 0, 1, 0, 1, 1, 1, 0)],
56 'ORCI': [(0, 0, -1, 0, -1, 0, -1, 0, 0),
57 (0, 0, 1, -1, 0, 0, -1, -1, 0),
58 (0, 1, 0, 1, 0, 0, 1, 1, -1),
59 (0, -1, 0, 1, 0, -1, 1, -1, 0)],
60 'RHL': [(0, -1, 0, 1, 1, 1, -1, 0, 0),
61 (1, 0, 0, 0, 1, 0, 0, 0, 1),
62 (1, -1, 0, 1, 0, -1, 1, 0, 0)],
63 'TET': [(1, 0, 0, 0, 1, 0, 0, 0, 1), (0, 1, 0, 0, 0, 1, 1, 0, 0)],
64 'MCL': [(0, 0, 1, -1, -1, 0, 1, 0, 0),
65 (-1, 0, 0, 0, 1, 0, 0, 0, -1),
66 (0, 0, -1, 1, 1, 0, 0, -1, 0),
67 (0, -1, 0, 1, 0, 1, -1, 0, 0),
68 (0, 1, 0, -1, 0, -1, 0, 0, 1),
69 (-1, 0, 0, 0, 1, 1, 0, 0, -1),
70 (0, 1, 0, 1, 0, -1, -1, 0, 0),
71 (0, 0, 1, 1, -1, 0, 0, 1, 0),
72 (0, 1, 0, -1, 0, 0, 0, 0, 1),
73 (0, 0, -1, -1, 1, 0, 1, 0, 0),
74 (1, 0, 0, 0, 1, -1, 0, 0, 1),
75 (0, -1, 0, -1, 0, 1, 0, 0, -1),
76 (-1, 0, 0, 0, -1, 1, 0, 1, 0),
77 (1, 0, 0, 0, -1, -1, 0, 1, 0),
78 (0, 0, -1, 1, 0, 0, 0, -1, 0)],
79 'MCLC': [(1, 1, 1, 1, 0, 1, 0, 0, -1),
80 (1, 1, 1, 1, 1, 0, -1, 0, 0),
81 (1, -1, 1, -1, 0, 1, 0, 0, -1),
82 (-1, 1, 0, 1, 0, 0, 0, 0, -1),
83 (1, 0, 0, 0, 1, 0, 0, 0, 1),
84 (-1, 0, -1, 1, -1, -1, 0, 0, 1),
85 (1, -1, -1, 1, -1, 0, -1, 0, 0),
86 (-1, -1, 0, -1, 0, -1, 1, 0, 0),
87 (1, 0, 1, 1, 0, 0, 0, 1, 0),
88 (-1, 1, 0, -1, 0, 1, 1, 0, 0),
89 (0, -1, 1, -1, 0, 1, 0, 0, -1),
90 (-1, -1, 0, -1, 0, 0, 0, 0, -1),
91 (-1, -1, 1, -1, 0, 1, 0, 0, -1),
92 (1, 0, 0, 0, -1, 1, 0, 0, -1),
93 (-1, 0, -1, 0, -1, -1, 0, 0, 1),
94 (1, 0, -1, -1, 1, -1, 0, 0, 1),
95 (1, -1, 1, 1, -1, 0, 0, 1, 0),
96 (0, -1, 0, 1, 0, -1, 0, 0, 1),
97 (-1, 0, 0, 1, 1, 1, 0, 0, -1),
98 (1, 0, -1, 0, 1, -1, 0, 0, 1),
99 (-1, 1, 0, 1, 1, -1, 0, -1, 0),
100 (1, 1, -1, 1, -1, 0, -1, 0, 0),
101 (-1, -1, -1, -1, -1, 0, 0, 1, 0),
102 (-1, 1, 1, 1, 0, 1, 0, 0, -1),
103 (-1, 0, 0, 0, -1, 0, 0, 0, 1),
104 (-1, -1, 1, 1, -1, 0, 0, 1, 0),
105 (1, 1, 0, -1, 0, -1, 0, 0, 1)],
106 'TRI': [(0, -1, 0, -1, 0, 0, 0, 0, -1),
107 (0, 1, 0, 0, 0, 1, 1, 0, 0),
108 (0, 0, -1, 0, -1, 0, -1, 1, 0),
109 (0, 0, 1, 0, 1, 0, -1, 0, 0),
110 (0, -1, 0, 0, 0, -1, 1, 1, 1),
111 (0, 1, 0, 0, 0, 1, 1, -1, 0),
112 (0, 0, -1, 0, -1, 0, -1, 0, 0),
113 (-1, 1, 0, 0, 0, -1, 0, -1, 0),
114 (0, 0, 1, 1, -1, 0, 0, 1, 0),
115 (0, 0, -1, 1, 1, 1, 0, -1, 0),
116 (-1, 0, 0, 0, 1, 0, 0, -1, -1),
117 (0, 0, 1, 1, 0, 0, 0, 1, 0),
118 (0, 0, 1, 0, 1, 0, -1, -1, -1),
119 (-1, 0, 0, 0, 0, -1, 0, -1, 0),
120 (0, -1, 0, 0, 0, -1, 1, 0, 0),
121 (1, 0, 0, 0, 1, 0, 0, 0, 1),
122 (0, 0, -1, -1, 0, 0, 1, 1, 1),
123 (0, 0, -1, -1, 0, 0, 0, 1, 0),
124 (-1, -1, -1, 0, 0, 1, 0, 1, 0)],
125}
127# XXX The TRI list was generated by looping over all TRI structures in
128# the COD (Crystallography Open Database) and seeing what operations
129# were necessary to map all those to standard form. Hence if the
130# data does not cover all possible inputs, we could miss something.
131#
132# Looping over all possible TRI lattices in general would generate
133# 100+ operations, we don't want to tabulate that.
136def lattice_loop(latcls, length_grid, angle_grid):
137 """Yield all lattices defined by the length and angle grids."""
138 param_grids = []
139 for varname in latcls.parameters:
140 # Actually we could choose one parameter, a, to always be 1,
141 # reducing the dimension of the problem by 1. The lattice
142 # recognition code should do something like that as well, but
143 # it doesn't. This could affect the impact of the eps value
144 # on lattice determination, so we just loop over the whole
145 # thing in order not to worry.
146 if latcls.name in ['MCL', 'MCLC']:
147 special_var = 'c'
148 else:
149 special_var = 'a'
150 if varname == special_var:
151 values = np.ones(1)
152 elif varname in 'abc':
153 values = length_grid
154 elif varname in ['alpha', 'beta', 'gamma']:
155 values = angle_grid
156 else:
157 raise ValueError(varname)
158 param_grids.append(values)
160 for latpars in itertools.product(*param_grids):
161 kwargs = dict(zip(latcls.parameters, latpars))
162 try:
163 lat = latcls(**kwargs)
164 except (UnconventionalLattice, AssertionError):
165 # XXX assertion error can happen because cellpar_to_cell
166 # makes certain assumptions. Should be investigated.
167 # {'b': 0.1, 'gamma': 60.0, 'c': 0.1, 'a': 1.0,
168 # 'alpha': 30.0, 'beta': 30.0} <-- this won't work
169 pass
170 else:
171 yield lat
174def find_niggli_ops(latcls, length_grid, angle_grid):
175 niggli_ops = {}
177 for lat in lattice_loop(latcls, length_grid, angle_grid):
178 cell = lat.tocell()
180 try:
181 rcell, op = cell.niggli_reduce()
182 except RuntimeError:
183 print('Niggli reduce did not converge')
184 continue
185 assert op.dtype == int
186 op_key = tuple(op.ravel())
188 if op_key in niggli_ops:
189 niggli_ops[op_key] += 1
190 else:
191 niggli_ops[op_key] = 1
193 rcell_test = Cell(op.T @ cell)
194 rcellpar_test = rcell_test.cellpar()
195 rcellpar = rcell.cellpar()
196 err = np.abs(rcellpar_test - rcellpar).max()
197 assert err < 1e-7, err
199 return niggli_ops
202def find_all_niggli_ops(length_grid, angle_grid, lattices=None):
203 all_niggli_ops = {}
204 if lattices is None:
205 lattices = [name for name in bravais_names
206 if name not in ['MCL', 'MCLC', 'TRI']]
208 for latname in lattices:
209 latcls = bravais_lattices[latname]
210 print(f'Working on {latname}...')
211 niggli_ops = find_niggli_ops(latcls, length_grid, angle_grid)
212 print(f'Found {len(niggli_ops)} ops for {latname}')
213 for key, count in niggli_ops.items():
214 print(f' {str(np.array(key)):>40}: {count}')
215 print()
216 all_niggli_ops[latname] = niggli_ops
217 return all_niggli_ops
220def generate_niggli_op_table(lattices=None,
221 length_grid=None,
222 angle_grid=None):
224 if length_grid is None:
225 length_grid = np.logspace(-0.5, 1.5, 50).round(3)
226 if angle_grid is None:
227 angle_grid = np.linspace(10, 179, 50).round()
228 all_niggli_ops_and_counts = find_all_niggli_ops(length_grid, angle_grid,
229 lattices=lattices)
231 niggli_op_table = {}
232 for latname, ops in all_niggli_ops_and_counts.items():
233 ops = [op for op in ops if np.abs(op).max() < 2]
234 niggli_op_table[latname] = ops
236 import pprint
237 print(pprint.pformat(niggli_op_table))
238 return niggli_op_table
241# For generation of the table, please see the test_bravais_type_engine
242# unit test. In case there's any trouble, some legacy code can be
243# found also in 6e2b1c6cae0ae6ee04638a9887821e7b1a1f2f3f .