Source code for mchammer.calculators.cluster_expansion_calculator

from typing import List, Union

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
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 for a particular (supercell) structure and commonly employed when setting up a Monte Carlo simulation, see :ref:`ensembles`. Cluster expansions, e.g., of the energy, typically yield property values *per site*. When running a Monte Carlo simulation one, 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. Parameters ---------- structure Structure for which to set up the calculator. cluster_expansion Cluster expansion from which to build calculator. name Human-readable identifier for this calculator. scaling Scaling factor applied to the property value predicted by the cluster expansion. use_local_energy_calculator Evaluate energy changes using only the local environment; this method is generally *much* faster. Unless you know what you are doing do *not* set this option to ``False``. """ def __init__(self, structure: Atoms, cluster_expansion: ClusterExpansion, name: str = 'Cluster Expansion Calculator', scaling: Union[float, int] = None, use_local_energy_calculator: bool = True) -> None: super().__init__(name=name) structure_cpy = structure.copy() cluster_expansion.prune() if cluster_expansion._cluster_space.is_supercell_self_interacting(structure): 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.') self.use_local_energy_calculator = use_local_energy_calculator self.cpp_calc = _ClusterExpansionCalculator( cluster_space=cluster_expansion.get_cluster_space_copy(), structure=Structure.from_atoms(structure_cpy), fractional_position_tolerance=cluster_expansion.fractional_position_tolerance) self._cluster_expansion = cluster_expansion if scaling is None: self._property_scaling = len(structure) else: self._property_scaling = scaling self._sublattices = self.cluster_expansion._cluster_space.get_sublattices(structure) @property def cluster_expansion(self) -> ClusterExpansion: """ Cluster expansion from which calculator was set up. """ return self._cluster_expansion
[docs] def calculate_total(self, *, occupations: List[int]) -> float: """ Calculates and returns the total property value of the current configuration. Parameters ---------- occupations The entire occupation vector (i.e., list of atomic species). """ cv = self.cpp_calc.get_cluster_vector(occupations) return np.dot(cv, self.cluster_expansion.parameters) * self._property_scaling
[docs] def calculate_change(self, *, sites: List[int], current_occupations: List[int], new_site_occupations: List[int]) -> float: """ Calculates and returns the sum of the contributions to the property due to the sites specified in :attr:`sites`. Parameters ---------- sites Indices of sites at which occupations will be changed. current_occupations Entire occupation vector (atomic numbers) before change. new_site_occupations Atomic numbers after change at the sites defined by :attr:`sites`. """ occupations = np.array(current_occupations) if not self.use_local_energy_calculator: 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 change = 0.0 try: for index, new_occupation in zip(sites, new_site_occupations): change += self._calculate_partial_change(occupations=occupations, flip_index=index, new_occupation=new_occupation) occupations[index] = new_occupation # Safe because we work with a copy except Exception as e: msg = 'Caught exception {}. Try setting parameter '.format(e) msg += 'use_local_energy_calculator to False in init' raise RuntimeError(msg) return change * self._property_scaling
def _calculate_partial_change(self, occupations: List[int], flip_index: int, new_occupation: int): """ Internal method to calculate the local contribution for one site index. Parameters ---------- occupations Entire occupation vector. flip_index Lattice index for site where a flipped has occurred. new_occupation Atomic number of new occupation at site :attr:`flip_index`. """ cv_change = self.cpp_calc.get_cluster_vector_change(occupations=occupations, flip_index=flip_index, new_occupation=new_occupation) return np.dot(cv_change, self.cluster_expansion.parameters) @property def sublattices(self) -> Sublattices: """ Sublattices of the calculator structure. """ return self._sublattices