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