from typing import Any
import numpy as np
from _icet import _ClusterExpansionCalculator
from icet.input_output.logging_tools import logger
from ase import Atoms
from icet import ClusterExpansion
from icet.core.structure import structure_to_arrays
from icet.core.sublattices import Sublattices
from mchammer.calculators.base_calculator import BaseCalculator
[docs]
class ClusterExpansionCalculator(BaseCalculator):
"""
A :class:`ClusterExpansionCalculator` object enables the efficient
calculation of properties described by a cluster expansion.
It is specific to a particular supercell and is commonly used when setting
up a Monte Carlo simulation, see :ref:`ensembles`.
Cluster expansions, for example of the energy, typically yield property
values *per site*.
A Monte Carlo simulation, however, considers changes in the *total* energy
of the system.
The default behavior is therefore to multiply the output of the cluster
expansion by the number of sites.
This behavior can be changed via the :attr:`scaling` keyword parameter.
The calculator holds the occupations of its supercell and evaluates
property changes against them.
They are taken from the structure it is constructed with, so every site of
that structure has to carry a species the cluster expansion defines.
The held occupations change only through an explicit accept.
The ensemble calls :func:`set_occupations` once at setup and
:func:`accept_change` for every accepted move, and
:func:`calculate_change` leaves the held occupations untouched, including
when it raises.
The compiled calculator behind this class, reachable as :attr:`cpp_calc`,
states its thread contract in its own docstring.
The section on the :ref:`calculator interface
<advanced_topics_calculator_interface>` describes how to drive the
calculator without an ensemble.
Parameters
----------
structure
Structure for which to set up the calculator.
cluster_expansion
Cluster expansion from which to build the calculator.
name
Human-readable identifier for this calculator.
scaling
Scaling factor applied to the property value predicted by the
cluster expansion.
By default the number of sites in :attr:`structure`.
use_local_energy_calculator
Evaluate changes using only the local environment of the changed sites.
This is generally *much* faster than evaluating the whole supercell
twice.
Unless you know what you are doing do *not* set this option to
``False``.
Example
-------
The following snippet sets up a calculator for an Ising-like cluster
expansion and evaluates the energy of a configuration as well as the
energy change of a swap.
The parameters are made up to keep the example self-contained::
>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace
>>> from mchammer.calculators import ClusterExpansionCalculator
>>> # prepare cluster expansion
>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
>>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])
>>> # prepare a configuration
>>> structure = prim.repeat(3)
>>> for k in range(5):
... structure[k].symbol = 'Ag'
>>> # the total energy of the configuration the calculator holds
>>> calculator = ClusterExpansionCalculator(structure, ce)
>>> occupations = structure.get_atomic_numbers()
>>> energy = calculator.calculate_total(occupations=occupations)
>>> # the energy change of swapping the species on sites 0 and 10,
>>> # which leaves the held configuration where it is
>>> sites = [0, 10]
>>> new_species = [occupations[10], occupations[0]]
>>> change = calculator.calculate_change(sites=sites,
... current_occupations=occupations,
... new_site_occupations=new_species)
>>> # adopt the swap
>>> calculator.accept_change(sites=sites, species=new_species)
"""
def __init__(self,
structure: Atoms, cluster_expansion: ClusterExpansion,
name: str = 'Cluster Expansion Calculator',
scaling: float | None = None,
use_local_energy_calculator: bool = True) -> None:
super().__init__(name=name)
structure_cpy = structure.copy()
self.use_local_energy_calculator = use_local_energy_calculator
# The calculator describes a model of its own, so that pruning it
# reaches nothing the caller holds.
self._cluster_expansion = cluster_expansion.copy()
self._cluster_expansion.prune()
# The compiled calculator takes a snapshot of the cluster space it is
# given, which is why it is built after the prune.
self.cpp_calc = _ClusterExpansionCalculator(
cluster_space=self._cluster_expansion._cluster_space,
**structure_to_arrays(structure_cpy),
fractional_position_tolerance=self._cluster_expansion.fractional_position_tolerance)
if self.cpp_calc.is_self_interacting:
logger.warning('The ClusterExpansionCalculator self-interacts, '
'which may lead to erroneous results. To avoid '
'self-interaction, use a larger supercell or a '
'cluster space with shorter cutoffs.')
# The parameter vector is the one the cluster expansion above holds, so
# that this calculator evaluates exactly the model it reports.
self._parameters = self._cluster_expansion._parameters
if scaling is None:
self._property_scaling = len(structure)
else:
self._property_scaling = scaling
self._sublattices = self._cluster_expansion._cluster_space.get_sublattices(structure)
def _get_rows(self) -> list[tuple[str, Any]]:
""" Label and value pairs describing this calculator. """
rows = super()._get_rows()
rows += [('number of sites', sum(len(sl.indices) for sl in self._sublattices)),
('number of parameters', len(self._parameters)),
('cutoffs', ' '.join('{:.4f}'.format(c)
for c in self._cluster_expansion.cutoffs)),
('scaling', self._property_scaling),
('use_local_energy_calculator', self.use_local_energy_calculator),
('self-interacting', self.cpp_calc.is_self_interacting)]
return rows
@property
def cluster_expansion(self) -> ClusterExpansion:
"""Cluster expansion this calculator evaluates (copy).
The orbits whose parameters are all zero are left out of it, see
:func:`ClusterExpansion.prune`.
It is handed out as a copy, so the calculator cannot be steered through
it.
Editing its parameters, pruning it, or editing a value nested in its
metadata changes neither what the calculator computes nor what the next
read of this property reports.
"""
return self._cluster_expansion.copy()
[docs]
def calculate_total(self, *, occupations: list[int]) -> float:
"""
Returns the total property value of a configuration.
The configuration is given in full and is not adopted, so the
configuration the calculator holds is left where it is.
Parameters
----------
occupations
The entire occupation vector by atomic number.
"""
cv = self.cpp_calc.get_cluster_vector(occupations)
return np.dot(cv, self._parameters) * self._property_scaling
[docs]
def set_occupations(self, occupations: list[int]) -> None:
"""
Sets the configuration this calculator describes.
An ensemble calls this method once at setup and keeps the calculator in
step afterwards through :func:`accept_change`.
Parameters
----------
occupations
The entire occupation vector by atomic number.
"""
self.cpp_calc.set_occupations(occupations)
[docs]
def accept_change(self, *, sites: list[int] | None = None,
species: list[int] | None = None) -> None:
"""
Advances the configuration this calculator describes by an accepted
change, which is the only way it moves during a simulation.
Calling this method without stating which change was accepted raises
instead of doing nothing.
Silently not advancing would leave the calculator evaluating every
subsequent change against a stale configuration, with no error
anywhere.
An empty list of sites is a legitimate empty change and does nothing.
Parameters
----------
sites
Indices of the sites whose occupations changed.
species
New occupations by atomic number on those sites.
Raises
------
TypeError
If :attr:`sites` or :attr:`species` is omitted.
ValueError
If :attr:`sites` and :attr:`species` differ in length.
"""
if sites is None or species is None:
raise TypeError('accept_change requires the sites and species of the '
'accepted change, since this calculator holds the '
'configuration and advances it only through accepted '
'changes.')
if len(sites) != len(species):
raise ValueError('sites and species must have the same length.')
if len(sites) == 0:
return
self.cpp_calc.apply_moves([[(int(site), int(occupation))
for site, occupation in zip(sites, species)]])
[docs]
def calculate_change(self, *, sites: list[int],
current_occupations: list[int],
new_site_occupations: list[int]) -> float:
"""
Returns the change of the property caused by changing the occupations
of the sites in :attr:`sites`.
The local calculation evaluates the change against the occupations this
calculator holds, which the ensemble keeps equal to
:attr:`current_occupations` through :func:`set_occupations` and
:func:`accept_change`.
When calling this method outside an ensemble, call
:func:`set_occupations` first.
The held occupations are left untouched, including when this method
raises.
Parameters
----------
sites
Indices of the sites whose occupations change.
current_occupations
The entire occupation vector by atomic number before the change.
The local calculation does not read it.
It is part of the calculator interface and is used when
:attr:`use_local_energy_calculator` is ``False``.
new_site_occupations
Atomic numbers after the change on the sites in :attr:`sites`.
Raises
------
ValueError
If :attr:`sites` and :attr:`new_site_occupations` differ in length.
"""
if len(sites) != len(new_site_occupations):
raise ValueError('sites and new_site_occupations must have the same'
' length, not {} and {}.'
.format(len(sites), len(new_site_occupations)))
if not self.use_local_energy_calculator:
occupations = np.array(current_occupations)
e_before = self.calculate_total(occupations=occupations)
occupations[sites] = np.array(new_site_occupations)
e_after = self.calculate_total(occupations=occupations)
return e_after - e_before
# The trial move is one move: an ordered sequence of
# (site, new occupation) pairs, evaluated sequentially by the
# batched entry point against the held occupations.
# Errors from the evaluation are raised as they are rather than
# rewritten as a suggestion to evaluate the whole supercell instead.
# They report an invalid site or a species the model does not define,
# which the full evaluation would reject in the same way.
move = [(int(site), int(occupation))
for site, occupation in zip(sites, new_site_occupations)]
changes = self.cpp_calc.get_cluster_vector_changes([move])
change = np.dot(changes[0], self._parameters)
return change * self._property_scaling
@property
def sublattices(self) -> Sublattices:
""" Sublattices of the structure the calculator describes. """
return self._sublattices
