"""
This module provides the kinetic Monte Carlo ensemble.
"""
import random
import numpy as np
from ase import Atoms
from ase.data import chemical_symbols as ase_chemical_symbols
from ase.units import kB
from .base_ensemble import BaseEnsemble
from ..calculators.cluster_expansion_calculator import ClusterExpansionCalculator
from ..configuration_manager import VACANCY_ATOMIC_NUMBER
from ..data_containers.kinetic_data_container import KineticDataContainer
from ..kinetics.event_table import EventTable
from ..kinetics.transition import Transition
#: Tolerance of the comparison between the running potential and its
#: recomputation at the ensemble data interval, relative to the magnitude
#: of the recomputed potential with an absolute floor of one.
SENTINEL_TOLERANCE = 1e-6
[docs]
class KineticMonteCarloEnsemble(BaseEnsemble):
"""Rejection-free kinetic Monte Carlo ensemble using the residence time
algorithm of Bortz, Kalos and Lebowitz.
Every step selects one hop with probability proportional to its rate,
applies it, and advances the physical time by an exponentially
distributed increment.
The rate of a candidate hop is
:math:`\\Gamma = \\nu \\exp(-E_a / k_B T)` with
:math:`E_a = E_\\mathrm{KRA} + \\Delta E / 2`, where the kinetically
resolved activation energy comes from the barrier model of the
transition and the energy difference between the endpoints comes from
the base cluster expansion.
Because the kinetically resolved activation energy is symmetric in the
endpoints, forward and reverse rates satisfy detailed balance against
the base cluster expansion exactly, for any parameters of the barrier
model, and the sampled equilibrium distribution is that of the base
cluster expansion.
A barrier model can predict a kinetically resolved activation energy
smaller than :math:`|\\Delta E| / 2`, which would give one direction
of the hop a negative barrier.
The kinetically resolved activation energy is therefore clamped from
below at :math:`|\\Delta E| / 2`, which makes the downhill direction
barrierless and gives the uphill direction the barrier
:math:`|\\Delta E|`.
The clamp raises both directions of a hop by the same amount and
preserves detailed balance, which a clamp of each barrier on its own
would break.
Every rate evaluation that hits the clamp is counted in
``clamp_count`` of the ensemble data.
An event is re-evaluated whenever a change nearby invalidates its
rate, and the count therefore measures how far the barrier model
reaches below the clamp instead of the number of clamped hops.
A nonzero count means the barrier model is used outside the range it
was fitted for.
A restart evaluates every event once and adds the clamped evaluations
of that refresh to the count.
A vacant site carries the species ``'X'`` of atomic number 0, which is
the convention of :mod:`mchammer` and the default of
:func:`~icet.core.event.find_hop_events` and
:func:`~icet.tools.local_environment_enumeration.enumerate_event_environments`.
A hop is available when the first site of the event carries the
migrating species and the second carries ``'X'``.
Besides the rate evaluation, the setup builds the invalidation lists
through one neighbor list query, and every step sums the rates of all
events to select one.
Both are linear in the supercell, and the sum dominates the step
beyond a few thousand sites.
Parameters
----------
structure
Atomic configuration to be used in the simulation, which also
defines the initial occupation vector, including the placement of
the migrating species and the vacancies.
calculator
Calculator for the base cluster expansion, used for the energy
difference between the endpoints of a hop.
It has to evaluate through its local energy calculator, since the
rate refresh crosses the boundary once per step with every
candidate hop.
transitions
Transition descriptors, one per hop type, each bundling a barrier
model and an attempt frequency.
temperature
Temperature in K.
boltzmann_constant
Boltzmann constant in eV/K.
user_tag
Human-readable tag for the ensemble.
random_seed
Seed for the random number generator.
dc_filename
Name of the file the data container is written to.
data_container_write_period
Period in seconds at which the data container is written to file.
ensemble_data_write_interval
Interval in steps at which ensemble data is written.
trajectory_write_interval
Interval in steps at which the occupation vector is written.
"""
def __init__(self,
structure: Atoms,
calculator: ClusterExpansionCalculator,
transitions: list[Transition],
temperature: float,
boltzmann_constant: float = kB,
user_tag: str | None = None,
random_seed: int | None = None,
dc_filename: str | None = None,
data_container_write_period: float = 600,
ensemble_data_write_interval: int | None = None,
trajectory_write_interval: int | None = None) -> None:
if not transitions:
raise ValueError('At least one transition is required')
# A temperature of zero would divide a clamped zero barrier by
# zero, and the resulting NaN rates would pass every guard of the
# step, so the temperature is validated here.
if not (np.isfinite(temperature) and temperature > 0):
raise ValueError('temperature must be a positive finite number'
f' ({temperature})')
if not (np.isfinite(boltzmann_constant) and boltzmann_constant > 0):
raise ValueError('boltzmann_constant must be a positive finite'
f' number ({boltzmann_constant})')
if not isinstance(calculator, ClusterExpansionCalculator):
raise TypeError('The calculator must be a'
' ClusterExpansionCalculator, since the rate'
' refresh evaluates every candidate hop through'
' its batched change evaluation'
f' ({type(calculator).__name__})')
if not calculator.use_local_energy_calculator:
raise ValueError('The batched change evaluation requires the'
' local calculator. Construct the calculator'
' with use_local_energy_calculator=True')
self._transitions = list(transitions)
self._temperature = float(temperature)
self._boltzmann_constant = boltzmann_constant
self._ensemble_parameters = dict(
temperature=temperature,
boltzmann_constant=boltzmann_constant)
# The ensemble draws from a generator of its own, which makes the
# sequence of a run a function of the seed alone.
# Other libraries draw from the generator the random module shares,
# among them the database layer of ase while the data container is
# written.
self._generator = random.Random()
# Fresh tracking state, overwritten on restart.
self._time = 0.0
self._clamp_count = 0
self._generator_state = None
self._particle_ids = np.arange(len(structure))
self._displacements = np.zeros((len(structure), 3))
self._hop_counts: dict = {}
self._maximum_time = None
super().__init__(
structure=structure,
calculator=calculator,
user_tag=user_tag,
random_seed=random_seed,
dc_filename=dc_filename,
data_container_class=KineticDataContainer,
data_container_write_period=data_container_write_period,
ensemble_data_write_interval=ensemble_data_write_interval,
trajectory_write_interval=trajectory_write_interval)
if self._generator_state is None:
self._generator.seed(self._random_seed)
else:
self._generator.setstate(self._generator_state)
# Particle i starts on site i, so the species of the particles are
# the occupations of the reference structure of the data container.
self._particle_numbers = self.data_container.structure.numbers
self._validate_supercells(structure)
# The event table treats the augmentation of a barrier model as a
# constant, which the validation establishes.
self._validate_augmentation()
self._events = EventTable(structure, self.calculator,
self._transitions, self._temperature,
self._boltzmann_constant)
self._potential = self.calculator.calculate_total(
occupations=self.configuration.occupations)
# ------------------------------------------------------------------
# setup
# ------------------------------------------------------------------
def _validate_supercells(self, structure: Atoms) -> None:
"""Checks that every calculator was built on the ensemble supercell.
The event catalog indexes the ensemble structure with the site
indices of the barrier calculators, and the energy difference of
every candidate hop comes from the base calculator against the
same indices, so a calculator built on a different supercell would
pair wrong sites without any index error.
"""
cell = np.array(structure.cell)
def matches(supercell) -> bool:
# Positions are compared modulo the lattice, so a supercell
# that differs only by periodic wrapping is accepted.
if len(supercell) != len(structure) \
or not np.allclose(supercell.cell, cell, atol=1e-8):
return False
deltas = supercell.get_scaled_positions(wrap=False) \
- structure.get_scaled_positions(wrap=False)
deltas -= np.rint(deltas)
return np.linalg.norm(deltas @ cell, axis=1).max() <= 1e-6
if not matches(self.calculator.structure):
raise ValueError(
'The calculator of the base cluster expansion was built on'
' a supercell that differs from the ensemble structure, so'
' the energy differences it returns belong to other sites')
for index, transition in enumerate(self._transitions):
if not matches(transition.calculator.structure):
raise ValueError(
f'The calculator of transition {index} was built on a'
' supercell that differs from the ensemble structure,'
' so its event instances do not index the ensemble'
' structure')
def _validate_augmentation(self) -> None:
"""Validates barrier-model sublattices and augmentation invariance.
Every anchor sublattice group of every transition must exist in
the base cluster expansion, which catches a barrier model built
for a different system.
Every augmentation feature must be invariant under the hops of
every transition of this simulation: a concentration feature on a
sublattice group is invariant under a hop exactly when both
anchors lie inside the group or both lie outside, since the
occupation multiset of the group is then conserved.
A single barrier model cannot see the other transitions, so the
check runs here, and the remedy for a violation is refitting the
model with a narrower augmentation.
"""
base_groups = {tuple(sorted(sublattice.chemical_symbols))
for sublattice in self.calculator.sublattices}
anchor_groups = []
for transition in self._transitions:
lcs = transition.calculator.local_cluster_expansion \
.local_cluster_space
symbols = lcs.chemical_symbols
groups = [tuple(sorted(symbols[index]))
for index in lcs.event.site_indices]
anchor_groups.append(groups)
for group in groups:
if group not in base_groups:
raise ValueError(
f'The anchor sublattice {group} of a transition does'
' not exist in the base cluster expansion, so the'
' barrier model and the base cluster expansion'
' describe different systems')
for transition_index, transition in enumerate(self._transitions):
lcs = transition.calculator.local_cluster_expansion \
.local_cluster_space
for feature in lcs.augmentation_features:
for other_index, groups in enumerate(anchor_groups):
if feature['order'] == 1:
# A concentration is conserved when both anchors
# lie inside the group or both lie outside.
group = feature['groups'][0]
inside = [anchor == group for anchor in groups]
violated = any(inside) and not all(inside)
else:
# An orbit correlation is conserved when no site
# of the orbit lies on an anchor group.
violated = bool(set(feature['groups'])
& set(groups))
if violated:
raise ValueError(
'The augmentation feature on sublattice'
f" {feature['sublattices']} of transition"
f' {transition_index} is not invariant under'
f' the hops of transition {other_index}.'
' Refit the barrier model with a narrower'
' augmentation')
# ------------------------------------------------------------------
# the step
# ------------------------------------------------------------------
def _next_random_number(self) -> float:
"""Returns the next random number from the generator of this
ensemble."""
return self._generator.random()
def _do_trial_step(self) -> int:
"""Carries out one kinetic Monte Carlo step.
Every step is accepted, which is the essential difference from the
thermodynamic ensembles.
"""
self._clamp_count += self._events.refresh()
chosen, total_rate = self._events.select(self._next_random_number())
# The complement keeps the argument of the logarithm positive,
# since the random number can be exactly zero.
time_increment = -np.log(1.0 - self._next_random_number()) \
/ total_rate
# Everything the step needs is read before the hop is applied, so
# that the only operation that can fail is the hop itself and a
# failure leaves the occupations, the tracking state and the
# running potential all where they were.
a, b = self._events.sites_of(chosen)
transition = self._events.transition_of(chosen)
vector = self._events.vector_of(chosen)
energy_change = self._events.energy_change_of(chosen)
symbol = ase_chemical_symbols[transition.migrating_number]
vacancy_symbol = ase_chemical_symbols[VACANCY_ATOMIC_NUMBER]
# Apply the hop, which is the only path that advances the
# occupation state of the calculator.
self.update_occupations([a, b], [VACANCY_ATOMIC_NUMBER,
transition.migrating_number])
self._potential += energy_change
self._events.mark_dirty(a)
self._events.mark_dirty(b)
# Track time, particle identity, and unwrapped displacement.
self._time += time_increment
atom = self._particle_ids[a]
vacancy = self._particle_ids[b]
self._displacements[atom] += vector
self._displacements[vacancy] -= vector
self._particle_ids[a], self._particle_ids[b] = vacancy, atom
self._hop_counts[symbol] = self._hop_counts.get(symbol, 0) + 1
self._hop_counts[vacancy_symbol] = \
self._hop_counts.get(vacancy_symbol, 0) + 1
return 1
def _run(self, number_of_trial_steps: int):
"""Runs the simulation for a number of steps without interruption.
The time criterion is enforced here, per step, because every
kinetic Monte Carlo step advances physical time.
A criterion checked only at the observer interval would overshoot
by up to a full interval of hops.
"""
for _ in range(number_of_trial_steps):
if self._terminate_sampling():
break
accepted = self._do_trial_step()
self._step += 1
self._accepted_trials += accepted
[docs]
def run(self, number_of_steps: int | None = None,
maximum_time: float | None = None):
"""Runs the kinetic Monte Carlo simulation.
At least one of the two criteria is required.
The run stops after the step that crosses :attr:`maximum_time`, so
the recorded time overshoots the target by one stochastic
increment.
Parameters
----------
number_of_steps
Number of hops to carry out.
maximum_time
Physical time in seconds at which to stop.
"""
if number_of_steps is None and maximum_time is None:
raise ValueError('Provide number_of_steps, maximum_time, or both')
if number_of_steps is not None \
and not isinstance(number_of_steps, (int, np.integer)):
raise TypeError('number_of_steps must be an integer'
f' ({number_of_steps})')
self._maximum_time = maximum_time
if number_of_steps is None:
number_of_steps = 2 ** 31
super().run(int(number_of_steps))
def _terminate_sampling(self) -> bool:
"""Returns whether the time criterion is met.
:func:`_run` checks it before every step, which bounds the
overshoot to one stochastic increment.
"""
return self._maximum_time is not None \
and self._time >= self._maximum_time
# ------------------------------------------------------------------
# observation, persistence, restart
# ------------------------------------------------------------------
@property
def temperature(self) -> float:
"""Temperature in K."""
return self._temperature
@property
def boltzmann_constant(self) -> float:
"""Boltzmann constant in eV/K."""
return self._boltzmann_constant
@property
def time(self) -> float:
"""Physical time in seconds."""
return self._time
def _get_ensemble_data(self) -> dict:
"""Returns the ensemble data of the current step.
The running potential is compared against recomputation here,
which is the sentinel that catches a desynchronization between the
configuration and the accumulated energy.
"""
recomputed = self.calculator.calculate_total(
occupations=self.configuration.occupations)
scale = max(1.0, abs(recomputed))
if abs(self._potential - recomputed) > SENTINEL_TOLERANCE * scale:
raise RuntimeError(
'The running potential deviates from recomputation'
f' ({self._potential} != {recomputed}), which indicates a'
' desynchronization between the configuration and the'
' accumulated energy')
data = {'potential': recomputed,
'time': self._time,
'clamp_count': self._clamp_count}
for number in np.unique(self._particle_numbers):
symbol = ase_chemical_symbols[number]
mask = self._particle_numbers == number
squared = np.sum(self._displacements[mask] ** 2, axis=1)
data[f'msd_{symbol}'] = float(squared.mean())
data[f'hops_{symbol}'] = self._hop_counts.get(symbol, 0)
return data
[docs]
def write_data_container(self, outfile: str | bytes):
"""Updates the last state of the simulation and writes the data
container to file.
The stored state includes the state of the generator of this
ensemble, from which a restart continues the sequence.
Parameters
----------
outfile
File to which to write.
"""
self._data_container._update_last_state(
last_step=self.step,
occupations=self.configuration.occupations.tolist(),
accepted_trials=self._accepted_trials,
random_state=random.getstate(),
generator_state=self._generator.getstate(),
time=self._time,
particle_ids=self._particle_ids.tolist(),
displacements=self._displacements.tolist(),
hop_counts=dict(self._hop_counts),
clamp_count=self._clamp_count)
self.data_container.write(outfile)
def _restart_ensemble(self):
"""Restarts the ensemble from the last state in the data container.
Time, particle identities and accumulated displacements cannot be
reconstructed from stored occupations and are restored here.
The rates follow from the restored occupations, and the first
step evaluates every event.
"""
super()._restart_ensemble()
last_state = self.data_container._last_state
self._time = float(last_state['time'])
self._particle_ids = np.array(last_state['particle_ids'])
self._displacements = np.array(last_state['displacements'])
self._hop_counts = dict(last_state['hop_counts'])
self._clamp_count = int(last_state['clamp_count'])
state = last_state['generator_state']
# A round trip through json turns the tuple of the state into a
# list, which setstate does not accept.
self._generator_state = (state[0], tuple(state[1]), state[2])
