Coverage for mchammer/observers/site_occupancy_observer.py: 100%
37 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-09-14 04:08 +0000
1from ase import Atoms
2from icet import ClusterSpace
3from icet.core.structure import Structure
4from mchammer.observers.base_observer import BaseObserver
5import numpy as np
8class SiteOccupancyObserver(BaseObserver):
9 """
10 This class represents a site occupation factor (SOF) observer.
12 A SOF observer allows to compute the site occupation factors along the
13 trajectory sampled by a Monte Carlo (MC) simulation.
15 Parameters
16 ----------
17 cluster_space
18 Cluster space from which the allowed species are extracted.
19 structure
20 Supercell consistent with primitive structure in :attr:`cluster_space`.
21 Used to determine which species are allowed on each site.
22 sites
23 Dictionary containing lists of sites that are to be considered.
24 The keys will be taken as the names of the sites; the indices refer
25 to the primitive structure associated with the cluster space.
26 interval
27 Observation interval. Defaults to ``None`` meaning that if the
28 observer is used in a Monte Carlo simulation, then the
29 :class:`Ensemble` object will determine the interval.
31 Example
32 -------
33 The following snippet illustrates how to use the site occupancy factor (SOF)
34 observer in a Monte Carlo simulation of a surface slab. Here, the SOF
35 observer is used to monitor the concentrations of different species at the
36 surface, the first subsurface layer, and the remaining 'bulk'. A minimal
37 cluster expansion is used with slightly modified surface interactions in
38 order to obtain an example that can be run without much ado. In practice,
39 one should of course use a proper cluster expansion::
41 >>> from ase.build import fcc111
42 >>> from icet import ClusterExpansion, ClusterSpace
43 >>> from mchammer.calculators import ClusterExpansionCalculator
44 >>> from mchammer.ensembles import CanonicalEnsemble
45 >>> from mchammer.observers import SiteOccupancyObserver
47 >>> # prepare reference structure
48 >>> prim = fcc111('Au', size=(1, 1, 10), vacuum=10.0)
49 >>> prim.translate((0.1, 0.1, 0.0))
50 >>> prim.wrap()
51 >>> prim.pbc = True # icet requires pbc in all directions
53 >>> # prepare cluster expansion
54 >>> cs = ClusterSpace(prim, cutoffs=[3.7], chemical_symbols=['Ag', 'Au'])
55 >>> params = [0] + 5 * [0] + 10 * [0.1]
56 >>> params[1] = 0.01
57 >>> params[6] = 0.12
58 >>> ce = ClusterExpansion(cs, params)
59 >>> print(ce)
61 >>> # prepare initial configuration based on a 2x2 supercell
62 >>> structure = prim.repeat((2, 2, 1))
63 >>> for k in range(20):
64 >>> structure[k].symbol = 'Ag'
66 >>> # set up MC simulation
67 >>> calc = ClusterExpansionCalculator(structure, ce)
68 >>> mc = CanonicalEnsemble(structure=structure, calculator=calc, temperature=600,
69 ... dc_filename='myrun_sof.dc')
71 >>> # set up observer and attach it to the MC simulation
72 >>> sites = {'surface': [0, 9], 'subsurface': [1, 8],
73 ... 'bulk': list(range(2, 8))}
74 >>> sof = SiteOccupancyObserver(cs, structure, sites, interval=len(structure))
75 >>> mc.attach_observer(sof)
77 >>> # run 1000 trial steps
78 >>> mc.run(1000)
80 After having run this snippet one can access the SOFs via the data
81 container::
83 >>> print(mc.data_container.data)
84 """
86 def __init__(self, cluster_space: ClusterSpace,
87 structure: Atoms,
88 sites: dict[str, list[int]],
89 interval: int = None) -> None:
90 super().__init__(interval=interval, return_type=dict, tag='SiteOccupancyObserver')
92 self._sites = {site: sorted(indices)
93 for site, indices in sites.items()}
95 self._set_allowed_species(cluster_space, structure)
97 def _set_allowed_species(self,
98 cluster_space: ClusterSpace,
99 structure: Atoms):
100 """
101 Set the allowed species for the selected sites in the Atoms object.
103 Parameters
104 ----------
105 cluster_space
106 Cluster space implicitly defining allowed species.
107 structure
108 Specific supercell (consistent with cluster_space) whose
109 allowed species are to be determined.
110 """
112 primitive_structure = Structure.from_atoms(cluster_space.primitive_structure)
113 chemical_symbols = cluster_space.chemical_symbols
115 if len(chemical_symbols) == 1:
116 # If the allowed species are the same for all sites no loop is
117 # required
118 allowed_species = {site: chemical_symbols[0] for
119 site in self._sites.keys()}
120 else:
121 # Loop over the lattice sites to find the allowed species
122 allowed_species = {}
123 for site, indices in self._sites.items():
124 allowed_species[site] = None
125 positions = structure.positions[np.array(indices)]
126 lattice_sites = primitive_structure.find_lattice_sites_by_positions(
127 positions=positions,
128 fractional_position_tolerance=cluster_space.fractional_position_tolerance)
129 for k, lattice_site in enumerate(lattice_sites):
130 species = chemical_symbols[lattice_site.index]
131 # check that the allowed species are equal for all sites
132 if allowed_species[site] is not None and \
133 species != allowed_species[site]:
134 raise Exception('The allowed species {} for the site'
135 ' with index {} differs from the'
136 ' result {} for the previous index'
137 ' ({})!'.format(species, indices[k],
138 allowed_species[site],
139 indices[k-1]))
140 allowed_species[site] = species
142 self._allowed_species = allowed_species
144 def get_observable(self, structure: Atoms) -> dict[str, list[float]]:
145 """
146 Returns the site occupation factors for a given atomic configuration.
148 Parameters
149 ----------
150 structure
151 Input atomic structure.
152 """
154 chemical_symbols = np.array(structure.get_chemical_symbols())
155 sofs = {}
156 for site, indices in self._sites.items():
157 counts = {species: 0 for species in self._allowed_species[site]}
158 symbols, sym_counts = np.unique(chemical_symbols[indices],
159 return_counts=True)
160 for sym, count in zip(symbols, sym_counts):
161 counts[sym] += count
163 for species in counts.keys():
164 key = 'sof_{}_{}'.format(site, species)
165 sofs[key] = float(counts[species]) / len(indices)
167 return sofs