from ase import Atoms
from mchammer.calculators.target_vector_calculator import TargetVectorCalculator
from .canonical_ensemble import CanonicalEnsemble
from .canonical_annealing import _cooling_exponential
import numpy as np
import random
from typing import Any
from icet.input_output.logging_tools import logger
from icet.input_output.repr_tools import html_representation, text_representation
logger = logger.getChild('target_cluster_vector_annealing')
[docs]
class TargetClusterVectorAnnealing:
"""
Instances of this class carry out simulated annealing towards a target
cluster vector.
Each trial step swaps two sites in one of the supercells and accepts the
swap with the Metropolis criterion applied to the objective function of the
:class:`TargetVectorCalculator
<mchammer.calculators.TargetVectorCalculator>` of that supercell.
The artificial temperature decreases exponentially from :attr:`T_start` to
:attr:`T_stop` over the course of the run, and the structure with the best
score seen during the run is kept.
Since it is impossible to know beforehand which supercell shape
accommodates the best match, several supercells can be annealed at the
same time.
The functions described in the section on :ref:`special quasirandom
structures <advanced_topics_sqs_structures>`, such as
:func:`generate_sqs <icet.tools.structure_generation.generate_sqs>` and
:func:`generate_target_structure
<icet.tools.structure_generation.generate_target_structure>`, set up the
calculators and this class for the common cases.
Parameters
----------
structure
Supercells to anneal.
Their occupations define the initial configurations.
calculators
One calculator per supercell in :attr:`structure`, in the same order.
T_start
Artificial temperature at which the annealing starts.
T_stop
Artificial temperature at which the annealing stops.
random_seed
Seed for the random number generator used in the Monte Carlo
simulation.
Raises
------
ValueError
If :attr:`structure` is a single structure instead of a list, or if the
number of structures differs from the number of calculators.
Example
-------
The following snippet generates a special quasirandom structure for a
binary alloy at equal concentrations.
The cluster vector of the ideal random alloy has all elements beyond the
zerolet equal to zero, which is the target::
>>> from ase.build import bulk
>>> from icet import ClusterSpace
>>> from icet.tools.structure_generation import occupy_structure_randomly
>>> from mchammer.calculators import TargetVectorCalculator
>>> from mchammer.ensembles import TargetClusterVectorAnnealing
>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[6.0], chemical_symbols=['Ag', 'Au'])
>>> target_vector = [1.0] + [0.0] * (len(cs) - 1)
>>> # two supercell shapes with random starting occupations
>>> supercells = [prim.repeat((2, 2, 2)), prim.repeat((4, 2, 1))]
>>> for supercell in supercells:
... occupy_structure_randomly(supercell, cs, {'Ag': 0.5, 'Au': 0.5})
>>> calculators = [TargetVectorCalculator(supercell, cs, target_vector)
... for supercell in supercells]
>>> annealing = TargetClusterVectorAnnealing(supercells, calculators,
... T_start=5.0, T_stop=0.001)
>>> sqs = annealing.generate_structure(number_of_trial_steps=2000)
>>> print(annealing.best_score)
"""
def __init__(self, structure: list[Atoms],
calculators: list[TargetVectorCalculator],
T_start: float = 5.0, T_stop: float = 0.001,
random_seed: int | None = None) -> None:
if isinstance(structure, Atoms):
raise ValueError(
'A list of ASE Atoms (supercells) must be provided')
if len(structure) != len(calculators):
raise ValueError('There must be as many supercells as there '
'are calculators ({} != {})'.format(len(structure),
len(calculators)))
logger.info('Initializing target cluster vector annealing '
'with {} supercells'.format(len(structure)))
# random number generator
if random_seed is None:
self._random_seed = random.randint(0, int(1e16))
else:
self._random_seed = random_seed
random.seed(a=self._random_seed)
# Initialize an ensemble for each supercell
sub_ensembles = []
for ens_id, (supercell, calculator) in enumerate(zip(structure, calculators)):
sub_ensembles.append(CanonicalEnsemble(structure=supercell,
calculator=calculator,
random_seed=random.randint(
0, int(1e16)),
user_tag='ensemble_{}'.format(
ens_id),
temperature=T_start,
dc_filename=None))
self._sub_ensembles = sub_ensembles
self._current_score = self._sub_ensembles[0].calculator.calculate_total(
occupations=self._sub_ensembles[0].configuration.occupations)
self._best_score = self._current_score
self._best_structure = structure[0].copy()
self._temperature = T_start
self._T_start = T_start
self._T_stop = T_stop
self._total_trials = 0
self._accepted_trials = 0
self._n_steps = 42
def _get_rows(self) -> list[tuple[str, Any]]:
"""Returns the label and value pairs from which the string and HTML
representations are built.
"""
return [('number of supercells', len(self._sub_ensembles)),
('T_start', self.T_start),
('T_stop', self.T_stop),
('temperature', self.temperature),
('n_steps', self.n_steps),
('total_trials', self.total_trials),
('accepted_trials', self.accepted_trials),
('current_score', self.current_score),
('best_score', self.best_score),
('random_seed', self._random_seed)]
def __str__(self) -> str:
""" String representation. """
return text_representation(title=self.__class__.__name__,
rows=self._get_rows(), width=60, label_width=22)
def _repr_html_(self) -> str:
""" HTML representation. Used, e.g., in jupyter notebooks. """
return html_representation(title=self.__class__.__name__, rows=self._get_rows())
def __repr__(self) -> str:
""" Representation. """
s = type(self).__name__ + '('
s += f'number_of_supercells={len(self._sub_ensembles)}'
s += f', temperature={self.temperature}'
s += ')'
return s
[docs]
def generate_structure(self, number_of_trial_steps: int | None = None) -> Atoms:
"""
Runs the annealing and returns the best structure found.
Parameters
----------
number_of_trial_steps
Total number of trial steps over all supercells.
By default 3000 steps per supercell.
Returns
-------
ase.Atoms
The structure with the best score seen during the run.
"""
if number_of_trial_steps is None:
self._n_steps = 3000 * len(self._sub_ensembles)
else:
self._n_steps = number_of_trial_steps
self._temperature = self._T_start
self._total_trials = 0
self._accepted_trials = 0
while self.total_trials < self.n_steps:
if self._total_trials % 1000 == 0:
logger.info('MC step {}/{} ({} accepted trials, '
'temperature {:.3f}), '
'best score: {:.3f}'.format(self.total_trials,
self.n_steps,
self.accepted_trials,
self.temperature,
self.best_score))
self._do_trial_step()
return self.best_structure
def _do_trial_step(self):
""" Carries out one Monte Carlo trial step. """
self._temperature = _cooling_exponential(
self.total_trials, self.T_start, self.T_stop, self.n_steps)
self._total_trials += 1
# Choose a supercell
ensemble = random.choice(self._sub_ensembles)
# Choose two sites and swap
sublattice_index = ensemble.get_random_sublattice_index(
ensemble._swap_sublattice_probabilities)
sites, species = ensemble.configuration.get_swapped_state(
sublattice_index)
# Update occupations so that the cluster vector (and its score)
# can be calculated
ensemble.configuration.update_occupations(sites, species)
new_score = ensemble.calculator.calculate_total(
occupations=ensemble.configuration.occupations)
if self._acceptance_condition(new_score - self.current_score):
self._current_score = new_score
self._accepted_trials += 1
# Since we are looking for the best structures we want to
# keep track of the best one we have found as yet (the
# current one may have a worse score)
if self._current_score < self._best_score:
self._best_structure = ensemble.structure
self._best_score = self._current_score
else:
ensemble.configuration.update_occupations(
sites, list(reversed(species)))
def _acceptance_condition(self, potential_diff: float) -> bool:
"""
Evaluates Metropolis acceptance criterion.
Parameters
----------
potential_diff
Change in the thermodynamic potential associated with the trial step.
"""
if potential_diff < 0:
return True
elif abs(self.temperature) < 1e-6: # temperature is numerically zero
return False
else:
p = np.exp(-potential_diff / self.temperature)
return p > random.random()
@property
def temperature(self) -> float:
""" Current artificial temperature. """
return self._temperature
@property
def T_start(self) -> float:
""" Artificial temperature at which the annealing starts. """
return self._T_start
@property
def T_stop(self) -> float:
""" Artificial temperature at which the annealing stops. """
return self._T_stop
@property
def n_steps(self) -> int:
""" Total number of trial steps of the run. """
return self._n_steps
@property
def total_trials(self) -> int:
""" Number of trial steps carried out so far. """
return self._total_trials
@property
def accepted_trials(self) -> int:
""" Number of trial steps accepted so far. """
return self._accepted_trials
@property
def current_score(self) -> float:
""" Score of the current configuration. """
return self._current_score
@property
def best_score(self) -> float:
""" Best score found so far. """
return self._best_score
@property
def best_structure(self) -> Atoms:
""" Structure with the best score found so far. """
return self._best_structure
