from collections.abc import Sequence
from typing import Dict, List, Tuple
import numpy as np
from ase import Atoms
from ase.data import chemical_symbols
from mchammer.observers.base_observer import BaseObserver
from icet.input_output.logging_tools import logger
logger = logger.getChild('structure_factor_observer')
[docs]class StructureFactorObserver(BaseObserver):
r"""This class represents a structure factor observer.
This observer allows one to compute structure factors along the
trajectory sampled by a Monte Carlo (MC) simulation. Structure
factors are convenient for monitoring long-range order. The
`structure factor <https://en.wikipedia.org/wiki/Structure_factor>`_
is defined as:
.. math::
S(\vec{q}) = \frac{1}{\sum_{j=1}^N f_j^2}
\sum_{j,k}^N e^{-i \vec{q} \cdot (\vec{R}_k - \vec{R}_j)}
In addition to this "total" structure factor, this observer
calculates pair-specific structure factors, which correspond to parts
of the summation defined above, with the summation restricted to pairs
of specific types, e.g., Au-Au, Au-Cu and Cu-Cu in the example below.
Parameters
----------
structure : Atoms
prototype for the structures for which the structure factor will
be computed later; the supercell size (but not its decoration)
must be identical; this structure is also used to determine the
the possible pairs if `symbol_pairs=None`
q_points : List[np.ndarray]
array of q-points at which to evaluate the structure factor;
the q-points need to be compatible with the supercell in order
for the structure factor to be real
symbol_pairs : Optional[List[Tuple[str, str]]]
list of symbol pairs for which structure factors will be computed,
e.g. [('Al', 'Cu'), ('Al', 'Al')]; if `None` (default) use all
pairs possible based on the input structure
form_factors : Optional[Dict[str, float]]
form factors for each atom type; this can be used to (very
coarsely) simulate X-ray or neutron spectra; note that in general
the form factors are q-dependent, see, e.g.,
`here <https://wiki.fysik.dtu.dk/ase/ase/xrdebye.html>`__ and
`here <https://wiki.fysik.dtu.dk/ase/ase/
xrdebye.html#ase.utils.xrdebye.XrDebye.get_waasmaier>`__;
by default (`None`) all form factors are set to 1
interval : int
the observation interval, defaults to None meaning that if the
observer is used in a Monte Carlo simulation, the `Ensemble` object
will set the interval.
Raises
------
ValueError
if q-point is not consistent with metric of input structure
Attributes
----------
tag : str
name of observer
interval : int
observation interval
Example
-------
The following snippet illustrates how to use the structure factor
observer in a simulated annealing run of dummy Cu-Au model to
observe the emergence of a long-range ordered L1_2 structure::
>>> import numpy as np
>>> from ase.build import bulk
>>> from icet import ClusterSpace, ClusterExpansion
>>> from mchammer.calculators import ClusterExpansionCalculator
>>> from mchammer.ensembles import CanonicalAnnealing
>>> from mchammer.observers import StructureFactorObserver
>>> # parameters
>>> size = 3
>>> alat = 4.0
>>> symbols = ['Cu', 'Au']
>>> # setup
>>> prim = bulk('Cu', a=alat, cubic=True)
>>> cs = ClusterSpace(prim, [0.9*alat], symbols)
>>> ce = ClusterExpansion(cs, [0, 0, 0.2])
>>> # make supercell
>>> supercell = prim.repeat(size)
>>> ns = int(0.25 * len(supercell))
>>> supercell.symbols[0:ns] = 'Au'
>>> np.random.shuffle(supercell.symbols)
>>> # define q-points to sample
>>> q_points = []
>>> q_points.append(2 * np.pi / alat * np.array([1, 0, 0]))
>>> q_points.append(2 * np.pi / alat * np.array([0, 1, 0]))
>>> q_points.append(2 * np.pi / alat * np.array([0, 0, 1]))
>>> # set up structure factor observer
>>> sfo = StructureFactorObserver(supercell, q_points)
>>> # run simulation
>>> calc = ClusterExpansionCalculator(supercell, ce)
>>> mc = CanonicalAnnealing(supercell, calc,
... T_start=900, T_stop=500, cooling_function='linear',
... n_steps=400*len(supercell))
>>> mc.attach_observer(sfo)
>>> mc.run()
After having run this snippet, one can access the structure factors via the data
container::
>>> dc = mc.data_container
>>> print(dc.data)
The emergence of the ordered low-temperature structure can be monitored by
following the temperature dependence of any of the pair-specific structure
factors.
"""
def __init__(self,
structure: Atoms,
q_points: List[Sequence],
symbol_pairs: List[Tuple[str, str]] = None,
form_factors: Dict[str, float] = None,
interval: int = None) -> None:
super().__init__(interval=interval, return_type=dict, tag='StructureFactorObserver')
self._original_structure = structure.copy()
self._q_points = q_points
self._Sq_lookup = self._get_Sq_lookup(structure, q_points)
if symbol_pairs is not None:
self._pairs = symbol_pairs
else:
self._pairs = [(e1, e2)
for e1 in set(structure.symbols)
for e2 in set(structure.symbols)]
self._pairs = sorted(set([tuple(sorted(p)) for p in self._pairs]))
if len(self._pairs) == 0:
raise ValueError('The StructureFactorObserver requires '
'at least one pair to be defined')
if len(self._pairs) == 1 and self._pairs[0][0] == self._pairs[0][1]:
logger.warning(f'Only one pair requested {self._pairs[0]}; '
'use either a structure occupied by more than '
'one species or use the symbol_pairs argument '
'to specify more pairs.')
self._unique_symbols = set(sorted([s for p in self._pairs for s in p]))
if form_factors is not None:
self._form_factors = form_factors
else:
self._form_factors = {s: 1 for s in self._unique_symbols}
self._form_factors = dict(sorted(self._form_factors.items()))
for symbol in self._unique_symbols:
if symbol not in chemical_symbols:
raise ValueError(f'Unknown species {symbol}')
if symbol not in self._form_factors:
raise ValueError(f'Form factor missing for {symbol}')
for symbol, ff in self._form_factors.items():
if ff == 0:
raise ValueError(f'Form factor for {symbol} is zero')
def _get_Sq_lookup(self,
structure: Atoms,
q_points: List[Sequence]) -> np.ndarray:
""" Get S(q) lookup data for a given supercell and q-points. """
n_atoms = len(structure)
dist_vectors = structure.get_all_distances(mic=True, vector=True)
Sq_lookup = np.zeros((n_atoms, n_atoms, len(self._q_points)), dtype=np.complex128)
for i in range(n_atoms):
for j in range(n_atoms):
Rij = dist_vectors[i][j]
Sq_lookup[i, j, :] = np.exp(-1j * np.dot(q_points, Rij))
return Sq_lookup
def _get_indices(self,
structure: Atoms) -> Dict[str, List[int]]:
"""Returns the indices of all atoms by type considering all unique
atom types in the structure used to instantiate the observer object.
"""
indices = dict()
symbols = structure.get_chemical_symbols()
for symbol in self._unique_symbols:
indices[symbol] = np.where(np.array(symbols) == symbol)[0]
return indices
def _compute_structure_factor(self,
structure: Atoms) -> Dict[Tuple[str, str], float]:
indices = self._get_indices(structure)
norm = np.sum([self._form_factors[sym] ** 2 * len(lst) for sym, lst in indices.items()])
if norm <= 0:
prefactor = 0 # occurs if none of the specified atom types are present in the structure
else:
prefactor = 1 / norm
Sq_dict = dict()
for sym1, sym2 in self._pairs:
inds1 = indices[sym1]
inds2 = indices[sym2]
norm = prefactor * self._form_factors[sym1] * self._form_factors[sym2]
Sq = np.zeros(len(self._q_points), dtype=np.complex128)
if len(inds1) != 0 and len(inds2) != 0:
for i in inds1:
Sq += self._Sq_lookup[i, inds2].sum(axis=0) * norm
if sym1 == sym2:
# the imaginary part must be zero because both ij and ji appear in the sum
Sq_dict[(sym1, sym2)] = Sq.real
else:
# since we only consider A-B (not B-A; see __init__) we explicitly add ij and ji,
# which comes down to adding the complex conjugate
Sq_dict[(sym1, sym2)] = np.real(Sq + Sq.conj())
return Sq_dict
[docs] def get_observable(self,
structure: Atoms) -> Dict[str, float]:
"""
Returns the structure factors for a given atomic configuration.
Parameters
----------
structure
input atomic structure
Raises
------
ValueError
if input structure is incompatible with structure used for initialization
"""
if len(structure) != len(self._original_structure):
raise ValueError('Input structure incompatible with structure used for initialization\n'
f' n_input: {len(structure)}\n'
f' n_original: {len(self._original_structure)}')
Sq_dict = self._compute_structure_factor(structure)
return_dict = dict()
for pair, Sq in Sq_dict.items():
for i in range(len(Sq)):
tag = 'sfo_{}_{}_q{}'.format(*pair, i)
return_dict[tag] = Sq[i]
return_dict[f'total_q{i}'] = return_dict.get(f'total_q{i}', 0) + Sq[i]
return_dict = dict(sorted(return_dict.items()))
return return_dict
@property
def form_factors(self) -> Dict[str, float]:
""" Form factors used in structure factor calculation """
return self._form_factors.copy()
@property
def q_points(self) -> List[np.ndarray]:
""" q-points for which structure factor is calculated """
return self._q_points.copy()
def __str__(self) -> str:
""" string representation of observer object. """
width = 60
name = self.__class__.__name__
s = [' {} '.format(name).center(width, '=')]
fmt = '{:15} : {}'
s += [fmt.format('tag', self.tag)]
s += [fmt.format('return_type', self.return_type)]
s += [fmt.format('interval', self.interval)]
s += [fmt.format('pairs', self._pairs)]
s += [fmt.format('form_factors', self.form_factors)]
for k, qpt in enumerate(self.q_points):
s += [fmt.format(f'q-point{k}', qpt)]
return '\n'.join(s)