Source code for mchammer.observers.cluster_count_observer
from collections.abc import Iterable
import pandas as pd
from ase import Atoms
from icet import ClusterSpace
from icet.core.local_orbit_list_generator import LocalOrbitListGenerator
from icet.core.structure import structure_to_arrays
from icet.tools.geometry import chemical_symbols_to_numbers
from mchammer.observers.base_observer import BaseObserver
[docs]
class ClusterCountObserver(BaseObserver):
"""This class represents a cluster count observer.
A cluster count observer keeps track of how often each decorated cluster
occurs along the trajectory sampled by a Monte Carlo (MC) simulation.
Given several canonical MC simulations at different temperatures, for
example, it gives access to the temperature dependence of the number of
nearest neighbors of a particular species.
The counts are stored in the data container in one column per decorated
cluster.
The columns are named ``0_Al``, ``0_Cu``, ``1_Al_Al``, ``1_Al_Cu``,
``1_Cu_Al``, ``1_Cu_Cu`` and so on, where the number is the orbit index and
the symbols are the species on the sites of the cluster.
The count is the number of clusters of that orbit in the structure that
carry that decoration.
Parameters
----------
cluster_space
Cluster space that defines the clusters to be counted.
structure
Supercell the observer works on.
interval
Observation interval in trial steps.
By default the ensemble the observer is attached to sets the interval
to the number of sites in the structure.
orbit_indices
Indices of the orbits whose clusters are counted.
By default all orbits are included.
Example
-------
The following snippet counts nearest and next-nearest neighbor pairs along
a canonical Monte Carlo trajectory of an Ising-like cluster expansion::
>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace
>>> from mchammer.calculators import ClusterExpansionCalculator
>>> from mchammer.ensembles import CanonicalEnsemble
>>> from mchammer.observers import ClusterCountObserver
>>> # 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 initial configuration
>>> structure = prim.repeat(3)
>>> for k in range(5):
... structure[k].symbol = 'Ag'
>>> # set up the simulation and attach the observer
>>> calculator = ClusterExpansionCalculator(structure, ce)
>>> mc = CanonicalEnsemble(structure=structure, calculator=calculator,
... temperature=600)
>>> observer = ClusterCountObserver(cs, structure, interval=len(structure))
>>> mc.attach_observer(observer)
>>> mc.run(1000)
>>> # the counts are available from the data container, for example the
>>> # number of nearest neighbor Ag-Au pairs along the trajectory
>>> counts = mc.data_container.get('1_Ag_Au')
The counts of a single structure are available as a data frame::
>>> print(observer.get_cluster_counts(structure))
"""
def __init__(self, cluster_space: ClusterSpace,
structure: Atoms,
interval: int | None = None,
orbit_indices: list[int] | None = None) -> None:
super().__init__(interval=interval, return_type=dict, tag='ClusterCountObserver')
self._cluster_space = cluster_space
local_orbit_list_generator = LocalOrbitListGenerator(
orbit_list=cluster_space.orbit_list,
**structure_to_arrays(structure),
fractional_position_tolerance=cluster_space.fractional_position_tolerance)
self._full_orbit_list = local_orbit_list_generator.generate_full_orbit_list()
if orbit_indices is None:
self._orbit_indices = list(range(len(self._full_orbit_list)))
elif not isinstance(orbit_indices, Iterable):
raise ValueError('Argument orbit_indices should be a list of integers, '
f'not {type(orbit_indices)}')
else:
self._orbit_indices = orbit_indices
self._possible_occupations = self._get_possible_occupations()
def _get_possible_occupations(self) -> dict[int, list[tuple[str]]]:
""" Returns a dictionary containing the possible occupations for each orbit. """
possible_occupations = {}
for i in self._orbit_indices:
possible_occupations_orbit = self._cluster_space.get_possible_orbit_occupations(i)
order = self._full_orbit_list.get_orbit(i).order
assert order == len(possible_occupations_orbit[0]), \
f'Order (n={order}) does not match possible occupations' \
f' (n={len(possible_occupations[0])}, {possible_occupations}).'
possible_occupations[i] = possible_occupations_orbit
return possible_occupations
[docs]
def get_cluster_counts(self, structure: Atoms) -> pd.DataFrame:
"""
Returns the number of clusters of each orbit and decoration in a
structure.
Parameters
----------
structure
Atomic configuration to count clusters in.
Returns
-------
pandas.DataFrame
One row per decorated cluster with the columns ``dc_tag`` (the
column name used in the data container), ``occupation`` (the
species on the sites of the cluster), ``cluster_count``,
``orbit_index``, and ``order``.
"""
rows = []
occupations = structure.get_atomic_numbers()
for i in self._orbit_indices:
orbit = self._full_orbit_list.get_orbit(i)
cluster_counts = orbit.get_cluster_counts(occupations)
for chemical_symbols in self._possible_occupations[i]:
count = cluster_counts.get(tuple(chemical_symbols_to_numbers(chemical_symbols)), 0)
row = {}
row['dc_tag'] = '{}_{}'.format(i, '_'.join(chemical_symbols))
row['occupation'] = chemical_symbols
row['cluster_count'] = count
row['orbit_index'] = i
row['order'] = orbit.order
rows.append(row)
return pd.DataFrame(rows)
[docs]
def get_observable(self, structure: Atoms) -> dict:
"""
Returns the cluster counts of a configuration.
Parameters
----------
structure
Atomic configuration to observe.
Returns
-------
dict
The number of clusters for each decorated cluster, keyed by the
column name used in the data container.
"""
counts = self.get_cluster_counts(structure)
count_dict = {row['dc_tag']: row['cluster_count']
for i, row in counts.iterrows()}
return count_dict
