Coverage for mchammer/observers/structure_factor_observer.py: 100%
101 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-09-14 04:08 +0000
1from collections.abc import Sequence
2import numpy as np
3from ase import Atoms
4from ase.data import chemical_symbols
5from mchammer.observers.base_observer import BaseObserver
6from icet.input_output.logging_tools import logger
7logger = logger.getChild('structure_factor_observer')
10class StructureFactorObserver(BaseObserver):
11 r"""This class represents a structure factor observer.
13 This observer allows one to compute structure factors along the
14 trajectory sampled by a Monte Carlo (MC) simulation. Structure
15 factors are convenient for monitoring long-range order. The
16 `structure factor <https://en.wikipedia.org/wiki/Structure_factor>`_
17 is defined as:
19 .. math::
21 S(\vec{q}) = \frac{1}{\sum_{j=1}^N f_j^2}
22 \sum_{j,k}^N e^{-i \vec{q} \cdot (\vec{R}_k - \vec{R}_j)}
24 In addition to this "total" structure factor, this observer
25 calculates pair-specific structure factors, which correspond to parts
26 of the summation defined above, with the summation restricted to pairs
27 of specific types, e.g., Au-Au, Au-Cu and Cu-Cu in the example below.
29 Parameters
30 ----------
31 structure
32 Prototype for the structures for which the structure factor will
33 be computed later. The supercell size (but not its decoration)
34 must be identical. This structure is also used to determine the
35 the possible pairs if :attr:`symbol_pairs=None`.
36 q_points
37 Array of q-points at which to evaluate the structure factor.
38 The q-points need to be compatible with the supercell in order
39 for the structure factor to be real.
40 symbol_pairs
41 List of symbol pairs for which structure factors will be computed,
42 e.g., ``[('Al', 'Cu'), ('Al', 'Al')]``. If ``None`` (default) use all
43 pairs possible based on the input structure.
44 form_factors
45 Form factors for each atom type. This can be used to
46 (coarsely) simulate X-ray or neutron spectra. Note that in
47 general the form factors are q-dependent, see, e.g., `here
48 <https://wiki.fysik.dtu.dk/ase/ase/xrdebye.html>`__ and `here
49 <https://wiki.fysik.dtu.dk/ase/ase/
50 xrdebye.html#ase.utils.xrdebye.XrDebye.get_waasmaier>`__. By
51 default (``None``) all form factors are set to 1.
52 interval
53 Observation interval. Defaults to ``None`` meaning that if the
54 observer is used in a Monte Carlo simulation, the :class:`Ensemble` object
55 will determine the interval.
57 Raises
58 ------
59 ValueError
60 If q-point is not consistent with metric of input structure.
62 Example
63 -------
64 The following snippet illustrates how to use the structure factor
65 observer in a simulated annealing run of dummy Cu-Au model to
66 observe the emergence of a long-range ordered :math:`L1_2` structure::
68 >>> import numpy as np
69 >>> from ase.build import bulk
70 >>> from icet import ClusterSpace, ClusterExpansion
71 >>> from mchammer.calculators import ClusterExpansionCalculator
72 >>> from mchammer.ensembles import CanonicalAnnealing
73 >>> from mchammer.observers import StructureFactorObserver
75 >>> # parameters
76 >>> size = 3
77 >>> alat = 4.0
78 >>> symbols = ['Cu', 'Au']
80 >>> # setup
81 >>> prim = bulk('Cu', a=alat, cubic=True)
82 >>> cs = ClusterSpace(prim, [0.9*alat], symbols)
83 >>> ce = ClusterExpansion(cs, [0, 0, 0.2])
85 >>> # make supercell
86 >>> supercell = prim.repeat(size)
87 >>> ns = int(0.25 * len(supercell))
88 >>> supercell.symbols[0:ns] = 'Au'
89 >>> np.random.shuffle(supercell.symbols)
91 >>> # define q-points to sample
92 >>> q_points = []
93 >>> q_points.append(2 * np.pi / alat * np.array([1, 0, 0]))
94 >>> q_points.append(2 * np.pi / alat * np.array([0, 1, 0]))
95 >>> q_points.append(2 * np.pi / alat * np.array([0, 0, 1]))
97 >>> # set up structure factor observer
98 >>> sfo = StructureFactorObserver(supercell, q_points)
100 >>> # run simulation
101 >>> calc = ClusterExpansionCalculator(supercell, ce)
102 >>> mc = CanonicalAnnealing(supercell, calc,
103 ... T_start=900, T_stop=500, cooling_function='linear',
104 ... n_steps=400*len(supercell))
105 >>> mc.attach_observer(sfo)
106 >>> mc.run()
108 After having run this snippet, one can access the structure factors via the data
109 container::
111 >>> dc = mc.data_container
112 >>> print(dc.data)
114 The emergence of the ordered low-temperature structure can be monitored by
115 following the temperature dependence of any of the pair-specific structure
116 factors.
118 """
120 def __init__(self,
121 structure: Atoms,
122 q_points: list[Sequence],
123 symbol_pairs: list[tuple[str, str]] = None,
124 form_factors: dict[str, float] = None,
125 interval: int = None) -> None:
126 super().__init__(interval=interval, return_type=dict, tag='StructureFactorObserver')
128 self._original_structure = structure.copy()
129 self._q_points = q_points
130 self._Sq_lookup = self._get_Sq_lookup(structure, q_points)
132 if symbol_pairs is not None:
133 self._pairs = symbol_pairs
134 else:
135 self._pairs = [(e1, e2)
136 for e1 in set(structure.symbols)
137 for e2 in set(structure.symbols)]
138 self._pairs = sorted(set([tuple(sorted(p)) for p in self._pairs]))
139 if len(self._pairs) == 0:
140 raise ValueError('The StructureFactorObserver requires '
141 'at least one pair to be defined')
142 if len(self._pairs) == 1 and self._pairs[0][0] == self._pairs[0][1]:
143 logger.warning(f'Only one pair requested {self._pairs[0]}; '
144 'use either a structure occupied by more than '
145 'one species or use the symbol_pairs argument '
146 'to specify more pairs.')
148 self._unique_symbols = set(sorted([s for p in self._pairs for s in p]))
149 if form_factors is not None:
150 self._form_factors = form_factors
151 else:
152 self._form_factors = {s: 1 for s in self._unique_symbols}
153 self._form_factors = dict(sorted(self._form_factors.items()))
155 for symbol in self._unique_symbols:
156 if symbol not in chemical_symbols:
157 raise ValueError(f'Unknown species {symbol}')
158 if symbol not in self._form_factors:
159 raise ValueError(f'Form factor missing for {symbol}')
161 for symbol, ff in self._form_factors.items():
162 if ff == 0:
163 raise ValueError(f'Form factor for {symbol} is zero')
165 def _get_Sq_lookup(self,
166 structure: Atoms,
167 q_points: list[Sequence]) -> np.ndarray:
168 """ Get S(q) lookup data for a given supercell and q-points. """
169 n_atoms = len(structure)
170 Sq_lookup = np.zeros((n_atoms, n_atoms, len(self._q_points)), dtype=np.complex128)
171 for i in range(n_atoms):
172 inds = [f for f in range(i, n_atoms)]
173 dist_vectors = structure.get_distances(i, inds, mic=True, vector=True)
174 for j in range(i, n_atoms):
175 Rij = dist_vectors[j - i]
176 Sq_lookup[i, j, :] = np.exp(-1j * np.dot(q_points, Rij))
177 Sq_lookup[j, i, :] = np.exp(-1j * np.dot(q_points, -Rij))
178 return Sq_lookup
180 def _get_indices(self,
181 structure: Atoms) -> dict[str, list[int]]:
182 """Returns the indices of all atoms by type considering all unique
183 atom types in the structure used to instantiate the observer object.
184 """
185 indices = dict()
186 symbols = structure.get_chemical_symbols()
187 for symbol in self._unique_symbols:
188 indices[symbol] = np.where(np.array(symbols) == symbol)[0]
189 return indices
191 def _compute_structure_factor(self,
192 structure: Atoms) -> dict[tuple[str, str], float]:
193 indices = self._get_indices(structure)
194 norm = np.sum([self._form_factors[sym] ** 2 * len(lst) for sym, lst in indices.items()])
195 if norm <= 0:
196 prefactor = 0 # occurs if none of the specified atom types are present in the structure
197 else:
198 prefactor = 1 / norm
200 Sq_dict = dict()
201 for sym1, sym2 in self._pairs:
202 inds1 = indices[sym1]
203 inds2 = indices[sym2]
204 norm = prefactor * self._form_factors[sym1] * self._form_factors[sym2]
205 Sq = np.zeros(len(self._q_points), dtype=np.complex128)
206 if len(inds1) != 0 and len(inds2) != 0:
207 for i in inds1:
208 Sq += self._Sq_lookup[i, inds2].sum(axis=0) * norm
209 if sym1 == sym2:
210 # the imaginary part must be zero because both ij and ji appear in the sum
211 Sq_dict[(sym1, sym2)] = Sq.real
212 else:
213 # since we only consider A-B (not B-A; see __init__) we explicitly add ij and ji,
214 # which comes down to adding the complex conjugate
215 Sq_dict[(sym1, sym2)] = np.real(Sq + Sq.conj())
216 return Sq_dict
218 def get_observable(self,
219 structure: Atoms) -> dict[str, float]:
220 """Returns the structure factors for a given atomic configuration.
222 Parameters
223 ----------
224 structure
225 Input atomic structure.
227 Raises
228 ------
229 ValueError
230 If the input structure is incompatible with structure used
231 for initialization of the :class:`StructureFactorObserver`.
232 """
233 if len(structure) != len(self._original_structure):
234 raise ValueError('Input structure incompatible with structure used for initialization\n'
235 f' n_input: {len(structure)}\n'
236 f' n_original: {len(self._original_structure)}')
237 Sq_dict = self._compute_structure_factor(structure)
238 return_dict = dict()
239 for pair, Sq in Sq_dict.items():
240 for i in range(len(Sq)):
241 tag = 'sfo_{}_{}_q{}'.format(*pair, i)
242 return_dict[tag] = Sq[i]
243 return_dict[f'total_q{i}'] = return_dict.get(f'total_q{i}', 0) + Sq[i]
244 return_dict = dict(sorted(return_dict.items()))
245 return return_dict
247 @property
248 def form_factors(self) -> dict[str, float]:
249 """ Form factors used in structure factor calculation. """
250 return self._form_factors.copy()
252 @property
253 def q_points(self) -> list[np.ndarray]:
254 """ q-points for which structure factor is calculated. """
255 return self._q_points.copy()
257 def __str__(self) -> str:
258 """ String representation of observer object. """
259 width = 60
260 name = self.__class__.__name__
261 s = [' {} '.format(name).center(width, '=')]
263 fmt = '{:15} : {}'
264 s += [fmt.format('tag', self.tag)]
265 s += [fmt.format('return_type', self.return_type)]
266 s += [fmt.format('interval', self.interval)]
267 s += [fmt.format('pairs', self._pairs)]
268 s += [fmt.format('form_factors', self.form_factors)]
269 for k, qpt in enumerate(self.q_points):
270 s += [fmt.format(f'q-point{k}', qpt)]
271 return '\n'.join(s)