import hashlib
import json
import random
import numpy as np
from ase import Atoms
from icet.core.sublattices import Sublattices
from icet.tools.geometry import atomic_number_to_chemical_symbol
class SwapNotPossibleError(Exception):
""" Raised when no swap trial move is possible on a sublattice. """
# Atomic number that ASE assigns to the vacancy species 'X'. A site is
# considered occupied when its occupation differs from this value.
VACANCY_ATOMIC_NUMBER = 0
def get_constraint_hash(neighbor_sites_to_avoid: dict[int, list[int]]) -> str:
"""Returns a short digest that identifies a neighbor constraint.
The digest lets a simulation that is restarted from a data container be
checked against the constraint the data container was written with, without
storing a mapping whose size grows with the system.
It is computed from a canonical form, so that a mapping written in a
different order gives the same digest.
Sites without neighbors are dropped, so that two mappings which impose the
same rules agree even when one of them lists unconstrained sites
explicitly.
Parameters
----------
neighbor_sites_to_avoid
Sites that must not be occupied simultaneously, keyed by site index.
"""
canonical = sorted((int(site), sorted(int(n) for n in neighbors))
for site, neighbors in neighbor_sites_to_avoid.items() if neighbors)
return hashlib.sha1(json.dumps(canonical).encode()).hexdigest()[:16]
[docs]
class ConfigurationManager(object):
"""
The configuration manager owns the configuration a Monte Carlo simulation
samples.
It holds the occupation vector together with the sites each species
occupies on each sublattice, which lets it propose swap and flip trial
moves in constant time through :func:`get_swapped_state` and
:func:`get_flip_state`.
Ensembles apply accepted moves through :func:`update_occupations` and
restore a saved configuration through :func:`set_occupations`.
The manager also enforces the neighbor constraint set through
:func:`set_neighbor_sites_to_avoid`.
Note
----
As a user you will usually not interact directly with objects of this type.
Parameters
----------
structure
Configuration to be handled.
sublattices
Sublattices that define the allowed occupations of the sites.
"""
def __init__(self, structure: Atoms, sublattices: Sublattices) -> None:
self._structure = structure.copy()
self._occupations = self._structure.numbers
self._sublattices = sublattices
self._sites_by_species = self._get_sites_by_species()
self._neighbor_sites_to_avoid: dict[int, list[int]] | None = None
def _get_sites_by_species(self) -> list[dict[int, list[int]]]:
"""
Returns the sites each species occupies, as one dictionary per
sublattice that maps an atomic number to the list of sites on that
sublattice which carry it.
"""
sites_by_species = []
for sl in self._sublattices:
species_dict = {key: [] for key in sl.atomic_numbers}
for site in sl.indices:
species_dict[self._occupations[site]].append(site)
sites_by_species.append(species_dict)
return sites_by_species
@property
def occupations(self) -> np.ndarray:
""" Occupation vector of the configuration (copy). """
return self._occupations.copy()
@property
def sublattices(self) -> Sublattices:
""" Sublattices of the configuration. """
return self._sublattices
@property
def structure(self) -> Atoms:
""" Atomic structure associated with configuration (copy). """
structure = self._structure.copy()
structure.set_atomic_numbers(self.occupations)
return structure
[docs]
def get_occupations_on_sublattice(self, sublattice_index: int) -> list[int]:
"""
Returns the occupations on one sublattice by atomic number.
Parameters
----------
sublattice_index
Index of the sublattice.
"""
sl = self.sublattices[sublattice_index]
return list(self.occupations[sl.indices])
[docs]
def is_swap_possible(self, sublattice_index: int,
allowed_species: list[int] | None = None) -> bool:
"""
Returns whether a swap trial move is possible on a sublattice, which
requires at least two different species on it.
Parameters
----------
sublattice_index
Index of the sublattice.
allowed_species
Atomic numbers of the species that may take part in the swap.
By default all species on the sublattice.
"""
sl = self.sublattices[sublattice_index]
if allowed_species is None:
swap_symbols = set(self.occupations[sl.indices])
else:
swap_symbols = set([o for o in self.occupations[sl.indices] if o in
allowed_species])
return len(swap_symbols) > 1
[docs]
def get_swapped_state(self, sublattice_index: int,
allowed_species: list[int] | None = None,
allowed_sites: list[int] | None = None
) -> tuple[list[int], list[int]]:
"""
Proposes a swap of two randomly chosen sites with different species.
The two sites are drawn from one sublattice, so the proposed
configuration respects the allowed occupations of the sites.
The configuration itself is not changed.
Parameters
----------
sublattice_index
Index of the sublattice from which to pick the sites.
allowed_species
Atomic numbers of the species that may take part in the swap.
By default all species on the sublattice.
allowed_sites
Indices of the sites that may take part in the swap.
By default all sites on the sublattice.
Returns
-------
tuple[list[int], list[int]]
The indices of the two sites and their occupations after the swap.
Raises
------
SwapNotPossibleError
If the sublattice holds no site to swap or only one species.
"""
# pick the first site
if allowed_species is None:
available_sites = self.sublattices[sublattice_index].indices
else:
available_sites = [
s for Z in allowed_species for s in
self._get_sites_by_species()[sublattice_index][Z]]
# only include allowed sites
if allowed_sites is not None:
available_sites = list(set(available_sites).intersection(allowed_sites))
try:
site1 = random.choice(available_sites)
except IndexError:
raise SwapNotPossibleError(f'Sublattice {sublattice_index} is empty.')
# pick the second site
if allowed_species is None:
possible_swap_species = \
set(self._sublattices.get_allowed_numbers_on_site(site1)) - \
set([self._occupations[site1]])
else:
possible_swap_species = \
set(allowed_species) - set([self._occupations[site1]])
possible_swap_sites = []
for Z in possible_swap_species:
possible_swap_sites.extend(self._sites_by_species[sublattice_index][Z])
# only include allowed sites
if allowed_sites is not None:
possible_swap_sites = list(set(possible_swap_sites).intersection(allowed_sites))
possible_swap_sites = np.array(possible_swap_sites)
try:
site2 = random.choice(possible_swap_sites)
except IndexError:
raise SwapNotPossibleError(
'Cannot swap on sublattice {} since it is full of {} species .'
.format(sublattice_index,
atomic_number_to_chemical_symbol([self._occupations[site1]])[0]))
return ([site1, site2], [self._occupations[site2], self._occupations[site1]])
[docs]
def get_flip_state(
self,
sublattice_index: int,
allowed_species: list[int] | None = None,
allowed_sites: list[int] | None = None
) -> tuple[int, int]:
"""
Proposes a flip of a randomly chosen site to another species allowed
on it.
The configuration itself is not changed.
Parameters
----------
sublattice_index
Index of the sublattice from which to pick the site.
allowed_species
Atomic numbers of the species that may take part in the flip.
By default all species on the sublattice.
allowed_sites
Indices of the sites that may take part in the flip.
By default all sites on the sublattice.
Returns
-------
tuple[int, int]
The index of the site and its occupation after the flip.
"""
if allowed_species is None:
available_sites = self._sublattices[sublattice_index].indices
else:
available_sites = [s for Z in allowed_species for s in
self._get_sites_by_species()[sublattice_index][Z]]
# only include allowed sites
if allowed_sites is not None:
available_sites = list(set(available_sites).intersection(allowed_sites))
site = random.choice(available_sites)
if allowed_species is not None:
species = random.choice(list(
set(allowed_species) - set([self._occupations[site]])))
else:
species = random.choice(list(
set(self._sublattices[sublattice_index].atomic_numbers) -
set([self._occupations[site]])))
return site, species
[docs]
def update_occupations(self, sites: list[int], species: list[int]) -> None:
"""
Changes the occupations of the given sites.
The whole update is checked before any of it is applied, so a rejected
update leaves the configuration as it was.
Parameters
----------
sites
Indices of the sites to change.
species
New occupations of those sites by atomic number.
Raises
------
ValueError
If the two lists differ in length, if a site appears twice, if a
site index is out of range, or if a species is not allowed on the
site it is given for.
"""
# The whole update is checked before any of it is applied, so a
# rejected update leaves the occupations and the sites by species in
# the state they were in. Applying as we go would leave the two
# disagreeing about every site the loop had already reached.
if len(sites) != len(species):
raise ValueError('sites and species must have the same length.')
if len(set(sites)) != len(sites):
raise ValueError('The same site must not appear more than once'
' in an update: {}'.format(sites))
for site, new_Z in zip(sites, species):
if site < 0 or site >= len(self._occupations):
raise ValueError('Site {} is not a valid site index'.format(site))
sublattice_index = self.sublattices.get_sublattice_index_from_site_index(site)
if new_Z not in self.sublattices[sublattice_index].atomic_numbers:
raise ValueError('Invalid new species {} on site {}'.format(new_Z, site))
for site, new_Z in zip(sites, species):
old_Z = self._occupations[site]
sublattice_index = self.sublattices.get_sublattice_index_from_site_index(site)
# Move the site from the list of sites for the old species to the
# list for the new one.
self._sites_by_species[sublattice_index][old_Z].remove(site)
self._sites_by_species[sublattice_index][new_Z].append(site)
# Update occupation vector itself
self._occupations[sites] = species
[docs]
def set_occupations(self, occupations: list[int]) -> None:
"""
Replaces the occupations of the whole configuration.
This is the absolute counterpart of :func:`update_occupations`, which
expresses a change relative to the current configuration.
It is what a restart needs, since restoring a saved configuration is
not the acceptance of a move.
The occupations are checked before anything changes, so a rejected
input leaves the configuration as it was.
Parameters
----------
occupations
New occupations by atomic number, one per site of the configuration.
Raises
------
ValueError
If the length does not match the configuration, or if a species is
not allowed on the site it is given for.
"""
if len(occupations) != len(self._occupations):
raise ValueError('occupations must have length {}, not {}'
.format(len(self._occupations), len(occupations)))
for site, new_Z in enumerate(occupations):
sublattice_index = self.sublattices.get_sublattice_index_from_site_index(site)
if new_Z not in self.sublattices[sublattice_index].atomic_numbers:
raise ValueError('Invalid new species {} on site {}'.format(new_Z, site))
# Assign in place, since the occupations are the atomic numbers of the
# structure this manager was built from.
self._occupations[:] = occupations
self._sites_by_species = self._get_sites_by_species()
@property
def neighbor_sites_to_avoid(self) -> dict[int, list[int]] | None:
""" Sites that must not be occupied simultaneously, keyed by site index (copy). """
if self._neighbor_sites_to_avoid is None:
return None
return {site: list(neighbors)
for site, neighbors in self._neighbor_sites_to_avoid.items()}
[docs]
def set_neighbor_sites_to_avoid(self,
neighbor_sites_to_avoid: dict[int, list[int]] | None
) -> None:
"""Sets the neighbor constraint that trial moves have to respect.
The constraint belongs to the configuration instead of to an individual
trial move, so that it cannot apply to some moves and not to others.
It governs the trial steps of :class:`ThermodynamicBaseEnsemble
<mchammer.ensembles.ThermodynamicBaseEnsemble>` that carry out a swap,
an SGC flip or a VCSGC flip.
Thermodynamic integration, Wang-Landau sampling and target cluster
vector annealing generate their trial moves differently and are not
subject to it.
Rejecting a trial move leaves the proposal unchanged and therefore
symmetric, so the acceptance criterion stays exact and the constraint
introduces no bias of its own.
That alone does not guarantee that the whole constrained space is
sampled: a restrictive mapping can leave the allowed configurations
disconnected under the available trial moves, in which case the
simulation only reaches the part it starts in.
Parameters
----------
neighbor_sites_to_avoid
Sites that must not be occupied simultaneously, keyed by site
index.
Sites that are absent from the mapping are unconstrained.
``None`` removes the constraint.
The mapping is copied, so changing it afterwards does not change
the constraint.
Raises
------
ValueError
If the constraint is not usable, see :func:`validate_constraint`.
"""
if neighbor_sites_to_avoid is None:
self._neighbor_sites_to_avoid = None
return
# a copy, so that mutating the mapping afterwards cannot get around the
# validation below; the indices are normalized on the way in, which also
# accepts the numpy integers that a neighbor list produces
constraint = {int(site): [int(n) for n in neighbors]
for site, neighbors in neighbor_sites_to_avoid.items()}
self.validate_constraint(constraint)
self._neighbor_sites_to_avoid = constraint
[docs]
def is_constraint_violated(self, sites: list[int], species: list[int]) -> bool:
"""Checks whether a trial move would violate the neighbor constraint.
The constraint forbids two sites that appear in each other's avoid list
from being occupied at the same time, where a site counts as occupied
when it is not held by a vacancy.
Provided that the current configuration satisfies the constraint, only
the sites touched by the trial move can introduce a violation, which is
what allows this check to be local.
Parameters
----------
sites
Indices of the sites that the trial move would change.
species
Occupations by atomic number that the trial move would assign to
:attr:`sites`.
Returns
-------
``True`` if applying the trial move would place occupants on two sites
that appear in each other's avoid list, ``False`` otherwise.
Always ``False`` when no constraint is set.
"""
if self._neighbor_sites_to_avoid is None:
return False
trial_occupations = dict(zip(sites, species))
for site, new_species in trial_occupations.items():
if new_species == VACANCY_ATOMIC_NUMBER:
continue
for neighbor in self._neighbor_sites_to_avoid.get(site, ()):
# a neighbor that the trial move also touches must be evaluated
# in its trial state rather than its current one
occupant = trial_occupations.get(neighbor, self._occupations[neighbor])
if occupant != VACANCY_ATOMIC_NUMBER:
return True
return False
[docs]
def validate_constraint(self, neighbor_sites_to_avoid: dict[int, list[int]]) -> None:
"""Checks that a neighbor constraint is usable for the current configuration.
The configuration has to satisfy the constraint already.
Each trial move is only checked against the sites it touches, which
assumes that the rest of the configuration satisfies the constraint.
A violating configuration is not stuck, since emptying an offending
site is never rejected, but the sampling is biased until the violations
happen to be cleared.
Parameters
----------
neighbor_sites_to_avoid
Sites that must not be occupied simultaneously, keyed by site index.
Sites that are absent from the mapping are unconstrained.
Raises
------
ValueError
If the constraint refers to a site that does not exist, if a site is
listed against itself, if it is not symmetric, or if the current
configuration already violates it.
"""
n_sites = len(self._occupations)
for site, neighbors in neighbor_sites_to_avoid.items():
for index in (site, *neighbors):
if not 0 <= index < n_sites:
raise ValueError(
f'neighbor_sites_to_avoid refers to site {index}, which does not'
f' exist in a configuration with {n_sites} sites.')
if site in neighbors:
raise ValueError(
f'neighbor_sites_to_avoid lists site {site} against itself, which'
' would leave it permanently unoccupiable.')
for neighbor in neighbors:
if site not in neighbor_sites_to_avoid.get(neighbor, ()):
raise ValueError(
'neighbor_sites_to_avoid must be symmetric: site'
f' {site} lists {neighbor}, but {neighbor} does not list {site}.')
for site, neighbors in neighbor_sites_to_avoid.items():
if self._occupations[site] == VACANCY_ATOMIC_NUMBER:
continue
for neighbor in neighbors:
if self._occupations[neighbor] != VACANCY_ATOMIC_NUMBER:
raise ValueError(
'The configuration already violates neighbor_sites_to_avoid:'
f' sites {site} and {neighbor} are both occupied. Each trial'
' move is only checked against the sites it touches, which'
' assumes that the rest of the configuration satisfies the'
' constraint, so starting here would sample a biased transient.')
