"""
This module provides the event prototype used to anchor local cluster
expansions, together with a helper that discovers symmetry-distinct hop
events in a primitive structure.
"""
from dataclasses import dataclass
from itertools import product
from typing import Any
import numpy as np
import spglib
from ase import Atoms
from ase.data import chemical_symbols as ase_chemical_symbols
from icet.input_output.repr_tools import (html_representation,
text_representation)
from icet.tools.geometry import call_spglib
[docs]
@dataclass(frozen=True)
class HopEvent:
"""Prototype of a hop event in the primitive structure.
A hop event is an ordered tuple of lattice sites, each given as a
basis index together with a unit cell offset in multiples of the
primitive lattice vectors, plus the species that migrates between
them.
The first site is placed in the (0, 0, 0) cell by convention.
An offset counts lattice vectors from the position of the basis site as
the primitive structure gives it, including a position outside the cell.
Attributes
----------
site_indices
Basis indices of the sites involved in the event.
offsets
Unit cell offsets of the sites, one triple per site.
species
Chemical symbol of the migrating species.
"""
site_indices: tuple
offsets: tuple
species: str
def __post_init__(self):
if len(self.site_indices) != len(self.offsets):
raise ValueError('site_indices and offsets must have the same length')
if len(self.site_indices) < 1:
raise ValueError('an event needs at least one site')
if any(len(offset) != 3 for offset in self.offsets):
raise ValueError('each offset must have three components')
def __len__(self) -> int:
return len(self.site_indices)
@property
def sites(self) -> list:
"""Lattice sites of the event as (basis index, offset) pairs."""
return list(zip(self.site_indices, self.offsets))
def _get_rows(self) -> list[tuple[str, Any]]:
"""Returns the label and value pairs that describe this event."""
return [('sites', self.sites),
('migrating species', self.species)]
def __str__(self) -> str:
""" String representation. """
return text_representation(title='Hop Event', rows=self._get_rows(),
width=60, label_width=17)
def _repr_html_(self) -> str:
""" HTML representation. Used, e.g., in jupyter notebooks. """
return html_representation(title='Hop Event', rows=self._get_rows())
[docs]
def positions(self, primitive_structure: Atoms) -> np.ndarray:
"""Returns the Cartesian positions of the event sites.
Parameters
----------
primitive_structure
Primitive structure in which the event is defined.
"""
return (primitive_structure.positions[list(self.site_indices)]
+ np.dot(self.offsets, primitive_structure.cell))
[docs]
def to_dict(self) -> dict:
"""Returns a dict representation, suitable for serialization."""
return {'site_indices': list(self.site_indices),
'offsets': [list(offset) for offset in self.offsets],
'species': self.species}
[docs]
@classmethod
def from_dict(cls, dct: dict) -> 'HopEvent':
"""Reconstructs an event from its dict representation."""
return cls(site_indices=tuple(dct['site_indices']),
offsets=tuple(tuple(offset) for offset in dct['offsets']),
species=dct['species'])
def get_offset_bounds(cell: np.ndarray, cutoff: float) -> np.ndarray:
"""Returns, per lattice direction, the largest cell offset that can
hold a point within the cutoff of the origin cell.
The bound per direction is the cutoff divided by the spacing of the
lattice planes normal to that direction.
The norm of the shortest cell vector overestimates that spacing for
a skewed cell, which is why it is not a valid bound.
Parameters
----------
cell
Cell matrix with the cell vectors as rows.
cutoff
Radius that the offsets must cover.
"""
inverse_cell = np.linalg.inv(np.array(cell))
plane_spacings = 1.0 / np.linalg.norm(inverse_cell, axis=0)
return np.ceil(cutoff / plane_spacings).astype(int) + 1
def find_sites_within_cutoff(scaled_positions: np.ndarray, cell: np.ndarray,
centers: list, basis_indices: list,
cutoff: float) -> list:
"""Returns the lattice sites of the given basis sites that lie within
the cutoff of the nearest center, excluding the centers.
Parameters
----------
scaled_positions
Scaled positions of the primitive basis, unwrapped.
cell
Cell matrix with the cell vectors as rows.
centers
Lattice sites as (basis index, offset) pairs.
basis_indices
Basis sites whose lattice sites are searched.
cutoff
Maximum distance from the nearest center.
"""
centers_scaled = np.array([scaled_positions[index] + np.asarray(offset)
for index, offset in centers])
center_positions = np.dot(centers_scaled, cell)
excluded = {(index, tuple(offset)) for index, offset in centers}
bounds = get_offset_bounds(cell, cutoff)
sites = []
for index in basis_indices:
# The bounds hold for sites less than one cell from a center.
# The range is taken relative to the offsets that bring this site
# closest to each center, which reduces any site to that case.
nearest = np.rint(centers_scaled - scaled_positions[index]).astype(int)
ranges = [range(low, high + 1) for low, high
in zip(nearest.min(axis=0) - bounds,
nearest.max(axis=0) + bounds)]
for offset in product(*ranges):
if (index, offset) in excluded:
continue
position = np.dot(scaled_positions[index] + offset, cell)
if np.linalg.norm(center_positions - position,
axis=1).min() <= cutoff:
sites.append((index, offset))
return sites
def get_allowed_species(primitive_structure: Atoms,
chemical_symbols: list) -> list:
"""Returns the list of allowed species per site of the primitive structure.
Parameters
----------
primitive_structure
Primitive structure.
chemical_symbols
Either a flat list of chemical symbols, in which case every site
allows the same species, or one list per site.
Raises
------
ValueError
If a site lists the same species twice, since a duplicate would
produce one point function too many and a degenerate basis.
"""
if all(isinstance(symbol, str) for symbol in chemical_symbols):
allowed = [list(chemical_symbols)] * len(primitive_structure)
else:
if len(chemical_symbols) != len(primitive_structure):
raise ValueError('chemical_symbols must have one entry per site'
f' ({len(chemical_symbols)} !='
f' {len(primitive_structure)})')
allowed = [list(symbols) for symbols in chemical_symbols]
for site, symbols in enumerate(allowed):
if len(set(symbols)) != len(symbols):
raise ValueError(f'Site {site} lists a species more than once'
f' ({symbols})')
return allowed
def get_sublattice_decorated_numbers(primitive_structure: Atoms,
chemical_symbols: list) -> np.ndarray:
"""Returns fake atomic numbers that encode the sublattice of each site.
Sites that allow the same set of species receive the same number, so
spglib applied to the decorated structure yields exactly the symmetry
operations that respect the sublattice structure.
Parameters
----------
primitive_structure
Primitive structure.
chemical_symbols
Allowed species, flat or per site.
"""
allowed = get_allowed_species(primitive_structure, chemical_symbols)
sorted_sets = sorted({tuple(sorted(symbols)) for symbols in allowed})
return np.array([sorted_sets.index(tuple(sorted(symbols))) + 1
for symbols in allowed])
def get_symmetry_operations(primitive_structure: Atoms,
chemical_symbols: list,
symprec: float = 1e-5) -> dict:
"""Returns the symmetry operations of the sublattice-decorated primitive structure.
The result is a dict with keys ``rotations`` and ``translations``, both
in fractional coordinates of the primitive cell.
Parameters
----------
primitive_structure
Primitive structure.
chemical_symbols
Allowed species, flat or per site.
symprec
Tolerance in Å passed to spglib.
"""
numbers = get_sublattice_decorated_numbers(primitive_structure, chemical_symbols)
cell = (np.array(primitive_structure.cell),
primitive_structure.get_scaled_positions(),
numbers)
return call_spglib(spglib.get_symmetry, cell, symprec=symprec)
def map_sites(sites: list, rotation: np.ndarray, translation: np.ndarray,
scaled_positions: np.ndarray, cell: np.ndarray,
symprec: float) -> list:
"""Returns the images of lattice sites under a symmetry operation.
A lattice site is a (basis index, offset) pair, and so is each image.
The image of a site is the nearest lattice site.
Parameters
----------
sites
Lattice sites to map.
rotation
Rotation of the operation in fractional coordinates.
translation
Translation of the operation in fractional coordinates.
scaled_positions
Scaled positions of the primitive basis, unwrapped.
cell
Cell matrix with the cell vectors as rows.
symprec
Tolerance of the symmetry analysis in Å.
Raises
------
ValueError
If an image lies farther than :attr:`symprec` from every lattice
site.
"""
images = []
for index, offset in sites:
image = np.dot(rotation, scaled_positions[index] + np.asarray(offset)) \
+ translation
deltas = image - scaled_positions
offsets = np.rint(deltas)
# spglib accepts an operation when it maps every basis site to
# within symprec of a basis site in Cartesian coordinates.
# The rotation maps the offset of a site to an integer offset, so
# the image of every lattice site has the error of its basis site.
distances = np.linalg.norm(np.dot(deltas - offsets, cell), axis=1)
basis = int(np.argmin(distances))
if distances[basis] > symprec:
raise ValueError('The image of a lattice site under a symmetry'
' operation does not lie on a lattice site')
images.append((basis, tuple(int(x) for x in offsets[basis])))
return images
def canonical_key(sites: list, unordered: bool = False) -> tuple:
"""Returns a key of lattice sites that is invariant under lattice
translations.
The sites are translated so that one of them sits in the (0, 0, 0)
cell.
An ordered key places the first site there and keeps the order.
An unordered key tries every site, sorts the translated sites, and
takes the smallest result, which makes it invariant under a
permutation of the sites as well.
Parameters
----------
sites
Lattice sites as (basis index, offset) pairs.
unordered
Whether the key ignores the order of the sites.
Returns
-------
tuple
(key, offset), where offset is the cell offset subtracted from
every site.
"""
candidates = []
for _, origin in (sites if unordered else sites[:1]):
shift = tuple(int(x) for x in origin)
translated = [(int(index),
tuple(int(x) - s for x, s in zip(offset, shift)))
for index, offset in sites]
if unordered:
translated.sort()
candidates.append((tuple(translated), shift))
return min(candidates)
[docs]
def find_hop_events(primitive_structure: Atoms,
chemical_symbols: list,
species: str,
cutoff: float,
vacancy_species: str = 'X',
symprec: float = 1e-5) -> list:
"""Finds the symmetry-distinct hop events of a primitive structure.
A candidate hop is an ordered pair of sites within :attr:`cutoff` of
each other such that both sites allow the migrating species and the
vacancy species.
Both species are required on both sites, since the migrating atom
sits on one site and the vacancy on the other, and the hop exchanges
them.
A binary alloy therefore needs the vacancy species among the allowed
species of every sublattice it migrates on.
Candidates that are equivalent under the space group of the
sublattice-decorated primitive structure, including exchange of the
endpoints, are reported once.
Parameters
----------
primitive_structure
Primitive structure.
chemical_symbols
Allowed species, flat or per site.
species
Chemical symbol of the migrating species.
cutoff
Maximum distance between the endpoints of a hop.
vacancy_species
Chemical symbol used for the vacant site.
The default ``'X'`` is the species of atomic number 0 that
:mod:`mchammer` treats as a vacancy, which is what
:class:`~mchammer.ensembles.KineticMonteCarloEnsemble` requires on
the target site of a hop.
symprec
Tolerance in Å passed to spglib, which also widens the cutoff.
"""
if species not in ase_chemical_symbols:
raise ValueError(f'Unknown species {species}')
if vacancy_species not in ase_chemical_symbols:
raise ValueError(f'Unknown vacancy species {vacancy_species}')
if species == vacancy_species:
raise ValueError(
f'The migrating species and the vacancy species are both'
f' {species}, so a hop would carry a vacancy onto a vacancy'
' and change nothing')
allowed = get_allowed_species(primitive_structure, chemical_symbols)
hosts = [index for index, symbols in enumerate(allowed)
if species in symbols and vacancy_species in symbols]
if not hosts:
raise ValueError(
f'No site of the primitive structure allows both {species} and'
f' the vacancy species {vacancy_species}, so the structure'
' hosts no hop of this species. Add the vacancy species to the'
' allowed species of its sublattice')
cell = np.array(primitive_structure.cell)
scaled_positions = primitive_structure.get_scaled_positions(wrap=False)
symmetry = get_symmetry_operations(primitive_structure, chemical_symbols, symprec)
# Candidate pairs with the first site in the home cell.
candidates = [((i, (0, 0, 0)), site) for i in hosts
for site in find_sites_within_cutoff(
scaled_positions, cell, [(i, (0, 0, 0))], hosts,
cutoff + symprec)]
# Deduplicate under the space group, including endpoint exchange.
distinct = []
seen = set()
for pair in candidates:
if canonical_key(pair, unordered=True)[0] in seen:
continue
# Mark every image of this pair as seen.
for rotation, translation in zip(symmetry['rotations'],
symmetry['translations']):
images = map_sites(pair, rotation, translation, scaled_positions,
cell, symprec)
seen.add(canonical_key(images, unordered=True)[0])
distinct.append(HopEvent(
site_indices=(pair[0][0], pair[1][0]),
offsets=(pair[0][1], pair[1][1]),
species=species))
return distinct
