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