Source code for mchammer.observers.site_occupancy_observer

from ase import Atoms
from icet import ClusterSpace
from icet.core.structure import Structure
from mchammer.observers.base_observer import BaseObserver
from typing import List, Dict
import numpy as np

[docs]class SiteOccupancyObserver(BaseObserver): """ This class represents a site occupation factor (SOF) observer. A SOF observer allows to compute the site occupation factors along the trajectory sampled by a Monte Carlo (MC) simulation. Parameters ---------- cluster_space : icet.ClusterSpace cluster space from which the allowed species are extracted structure : ase.Atoms supercell consistent with primitive structure in cluster space; used to determine which species are allowed on each site sites : dict(str, list(int)) dictionary containing lists of sites that are to be considered; the keys will be taken as the names of the sites; the indices refer to the primitive structure associated with the cluster space interval : int the observation interval, defaults to None meaning that if the observer is used in a Monte Carlo simulation, then the Ensemble object will set the interval. Attributes ---------- tag : str name of observer interval : int observation interval Example ------- The following snippet illustrate how to use the site occupancy factor (SOF) observer in a Monte Carlo simulation of a surface slab. Here, the SOF observer is used to monitor the concentrations of different species at the surface, the first subsurface layer, and the remaining 'bulk'. A minimal cluster expansion is used with slightly modified surface interactions in order to obtain an example that can be run without much ado. In practice, one should of course use a proper cluster expansion:: >>> from import fcc111 >>> from icet import ClusterExpansion, ClusterSpace >>> from mchammer.calculators import ClusterExpansionCalculator >>> from mchammer.ensembles import CanonicalEnsemble >>> from mchammer.observers import SiteOccupancyObserver >>> # prepare reference structure >>> prim = fcc111('Au', size=(1, 1, 10), vacuum=10.0) >>> prim.translate((0.1, 0.1, 0.0)) >>> prim.wrap() >>> prim.pbc = True # icet requires pbc in all directions >>> # prepare cluster expansion >>> cs = ClusterSpace(prim, cutoffs=[3.7], chemical_symbols=['Ag', 'Au']) >>> params = [0] + 5 * [0] + 10 * [0.1] >>> params[1] = 0.01 >>> params[6] = 0.12 >>> ce = ClusterExpansion(cs, params) >>> print(ce) >>> # prepare initial configuration based on a 2x2 supercell >>> structure = prim.repeat((2, 2, 1)) >>> for k in range(20): >>> structure[k].symbol = 'Ag' >>> # set up MC simulation >>> calc = ClusterExpansionCalculator(structure, ce) >>> mc = CanonicalEnsemble(structure=structure, calculator=calc, temperature=600, ... dc_filename='myrun_sof.dc') >>> # set up observer and attach it to the MC simulation >>> sites = {'surface': [0, 9], 'subsurface': [1, 8], ... 'bulk': list(range(2, 8))} >>> sof = SiteOccupancyObserver(cs, structure, sites, interval=len(structure)) >>> mc.attach_observer(sof) >>> # run 1000 trial steps >>> After having run this snippet one can access the SOFs via the data container:: >>> print( """ def __init__(self, cluster_space: ClusterSpace, structure: Atoms, sites: Dict[str, List[int]], interval: int = None) -> None: super().__init__(interval=interval, return_type=dict, tag='SiteOccupancyObserver') self._sites = {site: sorted(indices) for site, indices in sites.items()} self._set_allowed_species(cluster_space, structure) def _set_allowed_species(self, cluster_space: ClusterSpace, structure: Atoms): """ Set the allowed species for the selected sites in the Atoms object Parameters ---------- cluster_space Cluster space implicitly defining allowed species structure Specific supercell (consistent with cluster_space) whose allowed species are to be determined """ primitive_structure = Structure.from_atoms(cluster_space.primitive_structure) chemical_symbols = cluster_space.chemical_symbols if len(chemical_symbols) == 1: # If the allowed species are the same for all sites no loop is # required allowed_species = {site: chemical_symbols[0] for site in self._sites.keys()} else: # Loop over the lattice sites to find the allowed species allowed_species = {} for site, indices in self._sites.items(): allowed_species[site] = None positions = structure.positions[np.array(indices)] lattice_sites = primitive_structure.find_lattice_sites_by_positions( positions=positions, fractional_position_tolerance=cluster_space.fractional_position_tolerance) for k, lattice_site in enumerate(lattice_sites): species = chemical_symbols[lattice_site.index] # check that the allowed species are equal for all sites if allowed_species[site] is not None and \ species != allowed_species[site]: raise Exception('The allowed species {} for the site' ' with index {} differs from the' ' result {} for the previous index' ' ({})!'.format(species, indices[k], allowed_species[site], indices[k-1])) allowed_species[site] = species self._allowed_species = allowed_species
[docs] def get_observable(self, structure: Atoms) -> Dict[str, List[float]]: """ Returns the site occupation factors for a given atomic configuration. Parameters ---------- structure input atomic structure. """ chemical_symbols = np.array(structure.get_chemical_symbols()) sofs = {} for site, indices in self._sites.items(): counts = {species: 0 for species in self._allowed_species[site]} symbols, sym_counts = np.unique(chemical_symbols[indices], return_counts=True) for sym, count in zip(symbols, sym_counts): counts[sym] += count for species in counts.keys(): key = 'sof_{}_{}'.format(site, species) sofs[key] = float(counts[species]) / len(indices) return sofs