Source code for icet.tools.local_environment_enumeration
"""
This module provides a helper that enumerates training environments
around an event, for the stage of the workflow in which barriers are
computed with a first-principles or machine-learned method.
"""
import math
from itertools import product
import numpy as np
from ase import Atoms
from icet.core.local_cluster_space import LocalClusterSpace
from icet.tools.geometry import chemical_symbols_to_numbers
[docs]
def find_event_instance(local_cluster_space: LocalClusterSpace,
structure: Atoms) -> tuple:
"""Finds one instance of the event prototype in a supercell.
The instance returned is the one with the lowest site indices.
Parameters
----------
local_cluster_space
Local cluster space whose event prototype is to be located.
structure
Supercell in which to search.
"""
return min(local_cluster_space.find_event_instances(structure))
[docs]
def enumerate_event_environments(local_cluster_space: LocalClusterSpace,
background: Atoms,
event_sites: tuple = None,
enumeration_radius: float = None,
n_structures: int = None,
vacancy_species: str = 'X',
seed: int = 42) -> list:
"""Enumerates symmetrically distinct training environments around an event.
Decorations of the environment sites within :attr:`enumeration_radius`
of the event are enumerated, the background configuration is held
fixed beyond that radius, and duplicates under the point group of the
event are removed.
Exhaustive enumeration grows combinatorially with the number of
enumerated sites, so beyond a size threshold the helper falls back to
symmetry-reduced random sampling.
The anchor sites of the returned structures are decorated with the
migrating species on the first site and the vacancy species on the
remaining sites, which is the initial state of the hop for the NEB
stage of the workflow.
Parameters
----------
local_cluster_space
Local cluster space that defines event, environment, and symmetry.
background
Supercell whose occupations provide everything beyond the
enumerated sites.
event_sites
Indices of the anchor sites in the background.
Located automatically when not given.
enumeration_radius
Maximum distance of an enumerated site from the nearest anchor
site.
Defaults to the environment cutoff of the local cluster space.
n_structures
Maximum number of environments to return.
When the number of distinct environments exceeds this value a
random subset is returned.
Required when the number of candidate decorations exceeds 10^5.
vacancy_species
Chemical symbol used for the vacant anchor sites.
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.
seed
Seed for the random sampling.
Returns
-------
list
Pairs of (structure, event sites).
"""
lcs = local_cluster_space
if event_sites is None:
event_sites = find_event_instance(lcs, background)
else:
event_sites = tuple(int(site) for site in event_sites)
if enumeration_radius is None:
enumeration_radius = lcs.environment_cutoff
# The mapping call also runs the aliasing guard.
environment_indices = lcs.get_environment_site_mapping(background,
event_sites)
# Select the environment sites within the enumeration radius.
anchor_positions = lcs.event.positions(lcs.primitive_structure)
enumerated_sites = []
enumerated_indices = []
for site, supercell_index in zip(lcs.environment_sites,
environment_indices):
position = lcs.get_site_position(site)
distance = np.linalg.norm(anchor_positions - position, axis=1).min()
if distance <= enumeration_radius:
enumerated_sites.append(site)
enumerated_indices.append(supercell_index)
if not enumerated_sites:
raise ValueError('No environment sites within the enumeration radius')
# Permutations of the enumerated sites induced by the anchor operations.
permutations = lcs.get_site_permutations(enumerated_sites)
# Allowed occupations per enumerated site.
chemical_symbols = lcs.chemical_symbols
allowed = [sorted(chemical_symbols_to_numbers(chemical_symbols[site[0]]))
for site in enumerated_sites]
# math.prod computes the exact product, where a fixed-width product
# would wrap for 63 or more binary sites and dodge the guard below.
total = math.prod(len(numbers) for numbers in allowed)
if total > 10 ** 5 and n_structures is None:
raise ValueError(f'{total} candidate decorations. Pass n_structures'
' to sample instead of enumerating exhaustively')
rng = np.random.default_rng(seed)
if total <= 10 ** 5:
distinct = _distinct_labelings(product(*allowed), permutations)
if n_structures is not None and len(distinct) > n_structures:
chosen = rng.choice(len(distinct), size=n_structures,
replace=False)
distinct = [distinct[index] for index in sorted(chosen)]
else:
distinct = _sample_labelings(allowed, permutations, n_structures, rng)
# Build the structures.
results = []
migrating_number, vacancy_number = chemical_symbols_to_numbers(
[lcs.event.species, vacancy_species])
for labeling in distinct:
structure = background.copy()
numbers = structure.get_atomic_numbers()
for supercell_index, occupation in zip(enumerated_indices, labeling):
numbers[supercell_index] = occupation
numbers[event_sites[0]] = migrating_number
for site in event_sites[1:]:
numbers[site] = vacancy_number
structure.set_atomic_numbers(numbers)
results.append((structure, event_sites))
return results
def _canonical_labeling(labeling: tuple, permutations: list) -> tuple:
"""Returns the canonical representative of a labeling under permutations."""
return min(tuple(labeling[index] for index in permutation)
for permutation in permutations)
def _distinct_labelings(candidates, permutations: list) -> list:
"""Returns the distinct canonical labelings of an iterable of candidates."""
distinct = []
seen = set()
for labeling in candidates:
canonical = _canonical_labeling(tuple(labeling), permutations)
if canonical not in seen:
seen.add(canonical)
distinct.append(canonical)
return distinct
def _sample_labelings(allowed: list, permutations: list,
n_structures: int, rng) -> list:
"""Samples distinct canonical labelings at random."""
distinct = []
seen = set()
attempts = 0
while len(distinct) < n_structures and attempts < 100 * n_structures:
attempts += 1
labeling = tuple(int(rng.choice(numbers)) for numbers in allowed)
canonical = _canonical_labeling(labeling, permutations)
if canonical not in seen:
seen.add(canonical)
distinct.append(canonical)
return distinct
