Coverage for icet/tools/structure_enumeration.py: 96%
195 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1from itertools import product
2from typing import Dict, List, Tuple, Generator, Union
4import numpy as np
6from ase import Atoms
7from ase.build import make_supercell
8from spglib import get_symmetry
9from spglib import niggli_reduce as spg_nigg_red
11from icet.tools.geometry import ase_atoms_to_spglib_cell
12from .structure_enumeration_support.normal_form_matrices \
13 import HermiteNormalForm, yield_reduced_hnfs, SmithNormalForm, get_unique_snfs
14from .structure_enumeration_support.labeling_generation \
15 import LabelingGenerator
18def _translate_labelings(
19 labeling: tuple,
20 snf: SmithNormalForm,
21 nsites: int,
22 include_self: bool = False) -> Generator[Tuple[int], None, None]:
23 """
24 Yields labelings that are equivalent to the original labeling under
25 translations as dictated by the Smith normal form (SNF) provided.
27 Parameters
28 ----------
29 labeling
30 Labeling to be translated.
31 snf
32 SmithNormalForm object.
33 nsites
34 Number of sites in the primitive cell.
35 include_self
36 If ``True`` original labeling will be included.
37 """
39 # Compute size of each block within which translations occur
40 sizes = [nsites * block for block in snf.blocks]
42 # Loop over all possible translations within group as defined by snf
43 for trans in product(range(snf.S[0]), range(snf.S[1]), range(snf.S[2])):
44 if not include_self and sum(trans) == 0:
45 continue
47 labeling_trans = ()
48 for i in range(snf.S[0]):
49 group = (i + trans[0]) % snf.S[0]
50 block_i = labeling[sizes[0] * group:sizes[0] * (group + 1)]
51 for j in range(snf.S[1]):
52 group = (j + trans[1]) % snf.S[1]
53 block_j = block_i[sizes[1] * group:sizes[1] * (group + 1)]
54 for k in range(snf.S[2]):
55 group = (k + trans[2]) % snf.S[2]
56 labeling_trans += tuple(block_j[sizes[2] * group:
57 sizes[2] * (group + 1)])
58 yield labeling_trans
61def _get_all_labelings(snf: SmithNormalForm,
62 labeling_generator: LabelingGenerator,
63 nsites: int) -> List[tuple]:
64 """
65 Returns inequivalent labelings corresponding to a Smith normal form
66 (SNF) matrix. Superperiodic labelings as well as labelings that are
67 equivalent under translations for this particular SNF will not be
68 included. However, labelings that are equivalent by rotations that
69 leave the cell (but not the labeling) unchanged will still be
70 included, since these have to be removed for each Hermite normal
71 form (HNF) separately.
73 Parameters
74 ----------
75 snf
76 SmithNormalForm object.
77 labeling_generator
78 :class:`LabelingGenerator` object.
79 nsites
80 Number of sites per primitive cell.
81 """
82 labelings = []
83 for labeling in labeling_generator.yield_labelings(snf.ncells):
84 unique = True
85 for labeling_trans in _translate_labelings(labeling, snf, nsites,
86 include_self=False):
87 # Check whether it translates into itself. If so,
88 # then it has been added with a smaller cell.
89 if labeling == labeling_trans:
90 unique = False
91 break
93 # Check with previous labelings,
94 # if labeling can be translated into a previously
95 # added labeling, then it is not unique
96 if labeling_trans in labelings:
97 unique = False
98 break
99 if unique:
100 labelings.append(labeling)
101 return labelings
104def _permute_labeling(labeling: tuple, snf: SmithNormalForm,
105 transformation: List[np.ndarray],
106 nsites: int) -> Tuple[int]:
107 """
108 Returns permuted labeling according to transformations defined by
109 transformation.
111 Parameters
112 ----------
113 labeling
114 Labeling to be rotated.
115 snf
116 SmithNormalForm object.
117 transformation
118 Transformation consisting of rotation, translation and basis
119 shift.
120 nsites
121 Nnumber of sites in the primitive cell.
122 """
124 # Calculate transformation imposed by LRL multiplication
125 new_group_order = np.dot(snf.group_order, transformation[0].T)
127 # Loop over every atom to find its new position
128 labeling_new = [0] * len(labeling)
129 for member_index, member in enumerate(new_group_order):
131 # Transform according to Gp,
132 # but each site in the primitive cell also transforms in its own way
133 for basis in range(nsites):
134 new_cell = member + transformation[1][basis]
136 # Calculate new index, first by finding the right block,
137 # then by adding the basis index to that block
138 new_index = 0
139 for i in range(3):
140 new_index += (new_cell[i] % snf.S[i]) * snf.blocks[i] * nsites
141 new_index += transformation[2][basis]
143 # Add the contribution to the hash key
144 element = labeling[member_index * nsites + basis]
145 labeling_new[new_index] = element
146 return tuple(labeling_new)
149def _yield_unique_labelings(labelings: List[int], snf: SmithNormalForm,
150 hnf: HermiteNormalForm, nsites: int) -> tuple:
151 """
152 Yields labelings that are unique in every imaginable sense.
154 Parameters
155 ----------
156 labelkeys
157 List of hash keys to labelings that may still contain labelings
158 that are equivalent under rotations that leave the supercell
159 shape unchanged.
160 snf
161 SmithNormalForm object.
162 hnf
163 HermiteNormalForm object.
164 nsites
165 Number of sites in the primitive cell.
166 """
167 saved_labelings = []
168 for labeling in labelings:
170 # Check whether labeling is just a rotated version of a previous
171 # labeling. Apply transformation that is specific to the hnf
172 # and check all translations of the transformed labeling.
173 unique = True
174 for transformation in hnf.transformations:
176 labeling_rot = _permute_labeling(labeling, snf, transformation,
177 nsites)
179 # Commonly, the transformation leaves the labeling
180 # unchanged, so check that first as a special case
181 # (yields a quite significant speedup)
182 if labeling_rot == labeling:
183 continue
185 # Translate in all possible ways
186 for labeling_rot_trans in \
187 _translate_labelings(labeling_rot, snf, nsites,
188 include_self=True):
189 if labeling_rot_trans in saved_labelings:
190 # Then we have rotated and translated the labeling
191 # into one that was already yielded
192 unique = False
193 break
194 if not unique:
195 break
196 if unique:
197 # Then we have finally found a unique structure
198 # defined by an HNF matrix and a labeling
199 saved_labelings.append(labeling)
200 yield labeling
203def _labeling_to_ase_atoms(labeling: tuple, hnf: HermiteNormalForm,
204 cell: np.ndarray,
205 new_cell: np.ndarray,
206 basis: np.ndarray, chemical_symbols: List[str],
207 pbc: List[bool]) -> Atoms:
208 """
209 Returns structure object corresponding to the given labeling using
210 labeling, HNF matrix and parent lattice.
212 Parameters
213 ---------
214 labeling
215 Permutation of index of elements.
216 hnf
217 HNF object defining the supercell.
218 cell
219 Basis vectors of primtive cell listed row-wise.
220 new_cell
221 new cell shape
222 basis
223 Scaled coordinates to all sites in the primitive cell.
224 chemical_symbols
225 List of elements, e.g. ``['Au', 'Ag']``.
226 pbc
227 Periodic boundary conditions of the primitive structure.
228 """
229 a = hnf.H[0, 0]
230 b = hnf.H[1, 0]
231 c = hnf.H[1, 1]
232 d = hnf.H[2, 0]
233 e = hnf.H[2, 1]
234 f = hnf.H[2, 2]
235 symbols = []
236 positions = []
237 for z1 in range(a):
238 offset_1 = (b * z1) / a
239 for z2 in range(int(offset_1), int(offset_1 + c)):
240 offset_2 = z1 * (d - (e * b) / c) / a + (e * z2) / c
241 for z3 in range(int(offset_2), int(f + offset_2)):
242 for basis_vector_i, basis_vector in enumerate(basis):
244 # Determine position in Cartesian coordinates
245 pos = np.dot(cell.T, np.array([z1, z2, z3]) + basis_vector)
246 positions.append(pos)
248 # Now determine which species this is
249 g = np.dot(hnf.snf.L, [z1, z2, z3])
250 assert np.allclose(g, np.round(g))
251 g = np.round(g).astype('int64')
252 g = [g[i] % hnf.snf.S[i] for i in range(3)]
253 ind = len(basis) * (g[0] * hnf.snf.S[1] * hnf.snf.S[2] +
254 g[1] * hnf.snf.S[2] + g[2]) + basis_vector_i
255 symbols.append(chemical_symbols[labeling[ind]])
256 structure = Atoms(symbols, positions, cell=new_cell, pbc=(True, True, True))
257 structure.wrap()
258 structure.pbc = pbc
259 return structure
262def _get_symmetry_operations(structure: Atoms,
263 symprec: float,
264 position_tolerance: float) -> Dict[str, list]:
265 """
266 Returns symmetry operations permissable for a given structure as
267 obtained via `spglib <https://atztogo.github.io/spglib/>`_. The
268 symmetry operations consist of three parts: rotation, translation
269 and basis shifts. The latter define the way that sublattices
270 shift upon rotation (correponds to `d_Nd` in [HarFor09]_).
272 Parameters
273 ----------
274 structure
275 Structure for which symmetry operations are sought.
276 symprec
277 Tolerance imposed when analyzing the symmetry using spglib.
278 position_tolerance
279 Tolerance applied when comparing positions in Cartesian coordinates.
280 """
281 symmetries = get_symmetry(ase_atoms_to_spglib_cell(structure), symprec=symprec)
282 assert symmetries, ('spglib.get_symmetry() failed. Please make sure that'
283 ' the structure object is sensible.')
284 rotations = symmetries['rotations']
285 translations = symmetries['translations']
287 basis = structure.get_scaled_positions()
289 # Calculate how atoms within the primitive cell are shifted (from one site
290 # to another) and translated (from one primtive cell to another) upon
291 # operation with rotation matrix. Note that the translations are needed
292 # here because different sites translate differently.
293 basis_shifts = np.zeros((len(rotations), len(structure)), dtype='int64')
294 sites_translations = []
295 for i, rotation in enumerate(rotations):
296 translation = translations[i]
297 site_translations = []
298 for j, basis_element in enumerate(basis):
300 # Calculate how the site is transformed when operated on by
301 # symmetry of parent lattice (rotation and translation)
302 site_rot_trans = np.dot(rotation, basis_element) + translation
304 # The site may now have been moved to a different site in a
305 # different cell. We want to separate the two. (In HarFor09,
306 # basis_shift (site_rot_trans) corresponds to d_Nd and
307 # site_translation to t_Nd)
308 site_translation = [0, 0, 0]
309 for index in range(3):
310 while site_rot_trans[index] < -position_tolerance:
311 site_rot_trans[index] += 1
312 site_translation[index] -= 1
313 while site_rot_trans[index] > 1 - position_tolerance:
314 site_rot_trans[index] -= 1
315 site_translation[index] += 1
316 site_translations.append(site_translation)
318 # Find the basis element that the shifted basis correponds to
319 found = False
320 for basis_index, basis_element_comp in enumerate(basis):
321 distance = site_rot_trans - basis_element_comp
323 # Make sure that they do not differ with a basis vector
324 for dist_comp_i, dist_comp in enumerate(distance):
325 if abs(abs(dist_comp) - 1) < position_tolerance: 325 ↛ 326line 325 didn't jump to line 326, because the condition on line 325 was never true
326 distance[dist_comp_i] = 0
328 if (abs(distance) < position_tolerance).all():
329 assert not found
330 basis_shifts[i, j] = basis_index
331 found = True
332 assert found
334 sites_translations.append(np.array(site_translations))
336 symmetries['translations'] = sites_translations
337 symmetries['basis_shifts'] = basis_shifts
338 return symmetries
341def enumerate_structures(structure: Atoms,
342 sizes: Union[List[int], range],
343 chemical_symbols: list,
344 concentration_restrictions: dict = None,
345 niggli_reduce: bool = None,
346 symprec: float = 1e-5,
347 position_tolerance: float = None) -> Atoms:
348 """
349 Yields a sequence of enumerated structures. The function generates
350 *all* inequivalent structures that are permissible given a certain
351 lattice. Using the :attr:`chemical_symbols` and
352 :attr:`concentration_restrictions` keyword arguments it is possible to
353 specify which chemical_symbols are to be included on which site and
354 in which concentration range.
356 The function is sensitive to the boundary conditions of the input
357 structure. An enumeration of, for example, a surface can thus be
358 performed by setting :attr:`structure.pbc = [True, True, False]`.
360 The algorithm implemented here was developed by Gus L. W. Hart and
361 Rodney W. Forcade in Phys. Rev. B **77**, 224115 (2008)
362 [HarFor08]_ and Phys. Rev. B **80**, 014120 (2009) [HarFor09]_.
364 Parameters
365 ----------
366 structure
367 Primitive structure from which derivative superstructures should
368 be generated.
369 sizes
370 Number of sites (included in enumeration).
371 chemical_symbols
372 Chemical species with which to decorate the structure, e.g.,
373 ``['Au', 'Ag']``; see below for more examples.
374 concentration_restrictions
375 Allowed concentration range for one or more element in
376 :attr:`chemical_symbols`, e.g., ``{'Au': (0, 0.2)}`` will only
377 enumerate structures in which the Au content is between 0 and
378 20 %. Here, concentration is always defined as the number of
379 atoms of the specified kind divided by the number of *all*
380 atoms.
381 niggli_reduction
382 If ``True`` perform a Niggli reduction with spglib for each
383 structure. The default is ``True`` if :attr:`structure` is periodic in
384 all directions, ``False`` otherwise.
385 symprec
386 Tolerance imposed when analyzing the symmetry using spglib.
387 position_tolerance
388 Tolerance applied when comparing positions in Cartesian coordinates;
389 by default this value is set equal to :attr:`symprec`.
391 Examples
392 --------
394 The following code snippet illustrates how to enumerate structures
395 with up to 6 atoms in the unit cell for a binary alloy without any
396 constraints::
398 >>> from ase.build import bulk
399 >>> prim = bulk('Ag')
400 >>> for structure in enumerate_structures(structure=prim,
401 ... sizes=range(1, 5),
402 ... chemical_symbols=['Ag', 'Au']):
403 ... pass # Do something with the structure
405 To limit the concentration range to 10 to 40% Au the code should
406 be modified as follows::
408 >>> conc_restr = {'Au': (0.1, 0.4)}
409 >>> for structure in enumerate_structures(structure=prim,
410 ... sizes=range(1, 5),
411 ... chemical_symbols=['Ag', 'Au'],
412 ... concentration_restrictions=conc_restr):
413 ... pass # Do something with the structure
415 Often one would like to consider mixing on only one
416 sublattice. This can be achieved as illustrated for a
417 Ga(1-x)Al(x)As alloy as follows::
419 >>> prim = bulk('GaAs', crystalstructure='zincblende', a=5.65)
420 >>> for structure in enumerate_structures(structure=prim,
421 ... sizes=range(1, 9),
422 ... chemical_symbols=[['Ga', 'Al'], ['As']]):
423 ... pass # Do something with the structure
425 """
427 if position_tolerance is None: 427 ↛ 430line 427 didn't jump to line 430, because the condition on line 427 was never false
428 position_tolerance = symprec
430 nsites = len(structure)
431 basis = structure.get_scaled_positions()
433 # Construct descriptor of where species are allowed to be
434 if isinstance(chemical_symbols[0], str):
435 iter_chemical_symbols = [tuple(range(len(chemical_symbols)))] * nsites
436 elements = chemical_symbols
437 elif len(chemical_symbols) == nsites:
438 assert isinstance(chemical_symbols[0][0], str)
439 elements = []
440 for site in chemical_symbols:
441 for element in site:
442 if element not in elements:
443 elements.append(element)
444 iter_chemical_symbols = []
445 for site in chemical_symbols:
446 iter_chemical_symbols.append(tuple(elements.index(i)
447 for i in site))
448 else:
449 raise Exception('chemical_symbols needs to be a list of strings '
450 'or a list of list of strings.')
452 # Adapt concentration restrictions to iter_chemical_symbols
453 if concentration_restrictions:
454 concentrations = {}
455 for key, concentration_range in concentration_restrictions.items():
456 assert len(concentration_range) == 2, \
457 ('Each concentration range' +
458 ' needs to be specified as (c_low, c_high)')
459 if key not in elements: 459 ↛ 460line 459 didn't jump to line 460, because the condition on line 459 was never true
460 raise ValueError('{} found in concentration_restrictions but'
461 ' not in chemical_symbols'.format(key))
462 concentrations[elements.index(key)] = concentration_range
463 else:
464 concentrations = None
466 # Construct labeling generator
467 labeling_generator = LabelingGenerator(iter_chemical_symbols,
468 concentrations)
470 # Niggli reduce by default if all directions have
471 # periodic boundary conditions
472 if niggli_reduce is None: 472 ↛ 475line 472 didn't jump to line 475, because the condition on line 472 was never false
473 niggli_reduce = (sum(structure.pbc) == 3)
475 symmetries = _get_symmetry_operations(structure,
476 symprec=symprec,
477 position_tolerance=position_tolerance)
479 # Loop over each cell size
480 for ncells in sizes:
481 if ncells == 0: 481 ↛ 482line 481 didn't jump to line 482, because the condition on line 481 was never true
482 continue
484 hnfs = list(yield_reduced_hnfs(ncells, symmetries, structure.pbc))
485 snfs = get_unique_snfs(hnfs)
487 for snf in snfs:
488 labelings = _get_all_labelings(snf, labeling_generator, nsites)
489 for hnf in snf.hnfs:
490 if niggli_reduce:
491 new_cell = spg_nigg_red(np.dot(hnf.H.T, structure.cell))
492 if new_cell is None: 492 ↛ 493line 492 didn't jump to line 493, because the condition on line 492 was never true
493 new_cell = np.dot(hnf.H.T, structure.cell)
494 else:
495 new_cell = np.dot(hnf.H.T, structure.cell)
496 for labeling in _yield_unique_labelings(labelings, snf, hnf,
497 nsites):
498 yield _labeling_to_ase_atoms(labeling, hnf,
499 structure.cell, new_cell,
500 basis, elements,
501 structure.pbc)
504def enumerate_supercells(structure: Atoms,
505 sizes: Union[List[int], range],
506 niggli_reduce: bool = None,
507 symprec: float = 1e-5,
508 position_tolerance: float = None) -> Atoms:
509 """
510 Yields a sequence of enumerated supercells. The function generates
511 *all* inequivalent supercells that are permissible given a certain
512 lattice. Any supercell can be reduced to one of the supercells
513 generated.
515 The function is sensitive to the boundary conditions of the input
516 structure. An enumeration of, for example, a surface can thus be
517 performed by setting :attr:`structure.pbc = [True, True, False]`.
519 The algorithm is based on Gus L. W. Hart and
520 Rodney W. Forcade in Phys. Rev. B **77**, 224115 (2008)
521 [HarFor08]_ and Phys. Rev. B **80**, 014120 (2009) [HarFor09]_.
523 Parameters
524 ----------
525 structure
526 Primitive structure from which supercells should be generated.
527 sizes
528 Number of sites (included in enumeration).
529 niggli_reduction
530 If ``True`` perform a Niggli reduction with spglib for each
531 supercell. The default is ``True`` if :attr:`structure` is periodic in
532 all directions, ``False`` otherwise.
533 symprec
534 Tolerance imposed when analyzing the symmetry using spglib.
535 position_tolerance
536 Tolerance applied when comparing positions in Cartesian coordinates.
537 By default this value is set equal to :attr:`symprec`.
539 Examples
540 --------
542 The following code snippet illustrates how to enumerate supercells
543 with up to 6 atoms in the unit cell::
545 >>> from ase.build import bulk
546 >>> prim = bulk('Ag')
547 >>> for supercell in enumerate_supercells(structure=prim, sizes=range(1, 7)):
548 ... pass # Do something with the supercell
550 """
552 if position_tolerance is None: 552 ↛ 557line 552 didn't jump to line 557, because the condition on line 552 was never false
553 position_tolerance = symprec
555 # Niggli reduce by default if all directions have
556 # periodic boundary conditions
557 if niggli_reduce is None:
558 niggli_reduce = (sum(structure.pbc) == 3)
560 symmetries = _get_symmetry_operations(structure,
561 symprec=symprec,
562 position_tolerance=position_tolerance)
564 for ncells in sizes:
565 for hnf in yield_reduced_hnfs(ncells, symmetries, structure.pbc):
566 supercell = make_supercell(structure, hnf.H.T)
567 if niggli_reduce:
568 new_cell = spg_nigg_red(np.dot(hnf.H.T, structure.cell))
569 if new_cell is None: # Happens when spglib fails to Niggli reduce 569 ↛ 570line 569 didn't jump to line 570, because the condition on line 569 was never true
570 yield supercell
571 else:
572 Pprim = np.dot(new_cell, np.linalg.inv(structure.cell))
573 yield make_supercell(structure, Pprim)
574 else:
575 yield supercell