Coverage for mchammer/observers/structure_factor

1from collections.abc import Sequence

2from typing import Dict, List, Tuple

3import numpy as np

4from ase import Atoms

5from ase.data import chemical_symbols

6from mchammer.observers.base_observer import BaseObserver

7from icet.input_output.logging_tools import logger

8logger = logger.getChild('structure_factor_observer')

11class StructureFactorObserver(BaseObserver):

12 r"""This class represents a structure factor observer.

14 This observer allows one to compute structure factors along the

15 trajectory sampled by a Monte Carlo (MC) simulation. Structure

16 factors are convenient for monitoring long-range order. The

17 `structure factor <https://en.wikipedia.org/wiki/Structure_factor>`_

18 is defined as:

20 .. math::

22 S(\vec{q}) = \frac{1}{\sum_{j=1}^N f_j^2}

23 \sum_{j,k}^N e^{-i \vec{q} \cdot (\vec{R}_k - \vec{R}_j)}

25 In addition to this "total" structure factor, this observer

26 calculates pair-specific structure factors, which correspond to parts

27 of the summation defined above, with the summation restricted to pairs

28 of specific types, e.g., Au-Au, Au-Cu and Cu-Cu in the example below.

30 Parameters

31 ----------

32 structure

33 Prototype for the structures for which the structure factor will

34 be computed later. The supercell size (but not its decoration)

35 must be identical. This structure is also used to determine the

36 the possible pairs if :attr:`symbol_pairs=None`.

37 q_points

38 Array of q-points at which to evaluate the structure factor.

39 The q-points need to be compatible with the supercell in order

40 for the structure factor to be real.

41 symbol_pairs

42 List of symbol pairs for which structure factors will be computed,

43 e.g., ``[('Al', 'Cu'), ('Al', 'Al')]``. If ``None`` (default) use all

44 pairs possible based on the input structure.

45 form_factors

46 Form factors for each atom type. This can be used to

47 (coarsely) simulate X-ray or neutron spectra. Note that in

48 general the form factors are q-dependent, see, e.g., `here

49 <https://wiki.fysik.dtu.dk/ase/ase/xrdebye.html>`__ and `here

50 <https://wiki.fysik.dtu.dk/ase/ase/

51 xrdebye.html#ase.utils.xrdebye.XrDebye.get_waasmaier>`__. By

52 default (``None``) all form factors are set to 1.

53 interval

54 Observation interval. Defaults to ``None`` meaning that if the

55 observer is used in a Monte Carlo simulation, the :class:`Ensemble` object

56 will determine the interval.

58 Raises

59 ------

60 ValueError

61 If q-point is not consistent with metric of input structure.

63 Example

64 -------

65 The following snippet illustrates how to use the structure factor

66 observer in a simulated annealing run of dummy Cu-Au model to

67 observe the emergence of a long-range ordered :math:`L1_2` structure::

69 >>> import numpy as np

70 >>> from ase.build import bulk

71 >>> from icet import ClusterSpace, ClusterExpansion

72 >>> from mchammer.calculators import ClusterExpansionCalculator

73 >>> from mchammer.ensembles import CanonicalAnnealing

74 >>> from mchammer.observers import StructureFactorObserver

76 >>> # parameters

77 >>> size = 3

78 >>> alat = 4.0

79 >>> symbols = ['Cu', 'Au']

81 >>> # setup

82 >>> prim = bulk('Cu', a=alat, cubic=True)

83 >>> cs = ClusterSpace(prim, [0.9*alat], symbols)

84 >>> ce = ClusterExpansion(cs, [0, 0, 0.2])

86 >>> # make supercell

87 >>> supercell = prim.repeat(size)

88 >>> ns = int(0.25 * len(supercell))

89 >>> supercell.symbols[0:ns] = 'Au'

90 >>> np.random.shuffle(supercell.symbols)

92 >>> # define q-points to sample

93 >>> q_points = []

94 >>> q_points.append(2 * np.pi / alat * np.array([1, 0, 0]))

95 >>> q_points.append(2 * np.pi / alat * np.array([0, 1, 0]))

96 >>> q_points.append(2 * np.pi / alat * np.array([0, 0, 1]))

98 >>> # set up structure factor observer

99 >>> sfo = StructureFactorObserver(supercell, q_points)

100

101 >>> # run simulation

102 >>> calc = ClusterExpansionCalculator(supercell, ce)

103 >>> mc = CanonicalAnnealing(supercell, calc,

104 ... T_start=900, T_stop=500, cooling_function='linear',

105 ... n_steps=400*len(supercell))

106 >>> mc.attach_observer(sfo)

107 >>> mc.run()

108

109 After having run this snippet, one can access the structure factors via the data

110 container::

111

112 >>> dc = mc.data_container

113 >>> print(dc.data)

114

115 The emergence of the ordered low-temperature structure can be monitored by

116 following the temperature dependence of any of the pair-specific structure

117 factors.

118

119 """

120

121 def __init__(self,

122 structure: Atoms,

123 q_points: List[Sequence],

124 symbol_pairs: List[Tuple[str, str]] = None,

125 form_factors: Dict[str, float] = None,

126 interval: int = None) -> None:

127 super().__init__(interval=interval, return_type=dict, tag='StructureFactorObserver')

128

129 self._original_structure = structure.copy()

130 self._q_points = q_points

131 self._Sq_lookup = self._get_Sq_lookup(structure, q_points)

132

133 if symbol_pairs is not None:

134 self._pairs = symbol_pairs

135 else:

136 self._pairs = [(e1, e2)

137 for e1 in set(structure.symbols)

138 for e2 in set(structure.symbols)]

139 self._pairs = sorted(set([tuple(sorted(p)) for p in self._pairs]))

140 if len(self._pairs) == 0:

141 raise ValueError('The StructureFactorObserver requires '

142 'at least one pair to be defined')

143 if len(self._pairs) == 1 and self._pairs[0][0] == self._pairs[0][1]:

144 logger.warning(f'Only one pair requested {self._pairs[0]}; '

145 'use either a structure occupied by more than '

146 'one species or use the symbol_pairs argument '

147 'to specify more pairs.')

148

149 self._unique_symbols = set(sorted([s for p in self._pairs for s in p]))

150 if form_factors is not None:

151 self._form_factors = form_factors

152 else:

153 self._form_factors = {s: 1 for s in self._unique_symbols}

154 self._form_factors = dict(sorted(self._form_factors.items()))

155

156 for symbol in self._unique_symbols:

157 if symbol not in chemical_symbols:

158 raise ValueError(f'Unknown species {symbol}')

159 if symbol not in self._form_factors:

160 raise ValueError(f'Form factor missing for {symbol}')

161

162 for symbol, ff in self._form_factors.items():

163 if ff == 0:

164 raise ValueError(f'Form factor for {symbol} is zero')

165

166 def _get_Sq_lookup(self,

167 structure: Atoms,

168 q_points: List[Sequence]) -> np.ndarray:

169 """ Get S(q) lookup data for a given supercell and q-points. """

170 n_atoms = len(structure)

171 Sq_lookup = np.zeros((n_atoms, n_atoms, len(self._q_points)), dtype=np.complex128)

172 for i in range(n_atoms):

173 inds = [f for f in range(i, n_atoms)]

174 dist_vectors = structure.get_distances(i, inds, mic=True, vector=True)

175 for j in range(i, n_atoms):

176 Rij = dist_vectors[j - i]

177 Sq_lookup[i, j, :] = np.exp(-1j * np.dot(q_points, Rij))

178 Sq_lookup[j, i, :] = np.exp(-1j * np.dot(q_points, -Rij))

179 return Sq_lookup

180

181 def _get_indices(self,

182 structure: Atoms) -> Dict[str, List[int]]:

183 """Returns the indices of all atoms by type considering all unique

184 atom types in the structure used to instantiate the observer object.

185 """

186 indices = dict()

187 symbols = structure.get_chemical_symbols()

188 for symbol in self._unique_symbols:

189 indices[symbol] = np.where(np.array(symbols) == symbol)[0]

190 return indices

191

192 def _compute_structure_factor(self,

193 structure: Atoms) -> Dict[Tuple[str, str], float]:

194 indices = self._get_indices(structure)

195 norm = np.sum([self._form_factors[sym] ** 2 * len(lst) for sym, lst in indices.items()])

196 if norm <= 0:

197 prefactor = 0 # occurs if none of the specified atom types are present in the structure

198 else:

199 prefactor = 1 / norm

200

201 Sq_dict = dict()

202 for sym1, sym2 in self._pairs:

203 inds1 = indices[sym1]

204 inds2 = indices[sym2]

205 norm = prefactor * self._form_factors[sym1] * self._form_factors[sym2]

206 Sq = np.zeros(len(self._q_points), dtype=np.complex128)

207 if len(inds1) != 0 and len(inds2) != 0:

208 for i in inds1:

209 Sq += self._Sq_lookup[i, inds2].sum(axis=0) * norm

210 if sym1 == sym2:

211 # the imaginary part must be zero because both ij and ji appear in the sum

212 Sq_dict[(sym1, sym2)] = Sq.real

213 else:

214 # since we only consider A-B (not B-A; see __init__) we explicitly add ij and ji,

215 # which comes down to adding the complex conjugate

216 Sq_dict[(sym1, sym2)] = np.real(Sq + Sq.conj())

217 return Sq_dict

218

219 def get_observable(self,

220 structure: Atoms) -> Dict[str, float]:

221 """Returns the structure factors for a given atomic configuration.

222

223 Parameters

224 ----------

225 structure

226 Input atomic structure.

227

228 Raises

229 ------

230 ValueError

231 If the input structure is incompatible with structure used

232 for initialization of the :class:`StructureFactorObserver`.

233 """

234 if len(structure) != len(self._original_structure):

235 raise ValueError('Input structure incompatible with structure used for initialization\n'

236 f' n_input: {len(structure)}\n'

237 f' n_original: {len(self._original_structure)}')

238 Sq_dict = self._compute_structure_factor(structure)

239 return_dict = dict()

240 for pair, Sq in Sq_dict.items():

241 for i in range(len(Sq)):

242 tag = 'sfo_{}_{}_q{}'.format(*pair, i)

243 return_dict[tag] = Sq[i]

244 return_dict[f'total_q{i}'] = return_dict.get(f'total_q{i}', 0) + Sq[i]

245 return_dict = dict(sorted(return_dict.items()))

246 return return_dict

247

248 @property

249 def form_factors(self) -> Dict[str, float]:

250 """ Form factors used in structure factor calculation. """

251 return self._form_factors.copy()

252

253 @property

254 def q_points(self) -> List[np.ndarray]:

255 """ q-points for which structure factor is calculated. """

256 return self._q_points.copy()

257

258 def __str__(self) -> str:

259 """ String representation of observer object. """

260 width = 60

261 name = self.__class__.__name__

262 s = [' {} '.format(name).center(width, '=')]

263

264 fmt = '{:15} : {}'

265 s += [fmt.format('tag', self.tag)]

266 s += [fmt.format('return_type', self.return_type)]

267 s += [fmt.format('interval', self.interval)]

268 s += [fmt.format('pairs', self._pairs)]

269 s += [fmt.format('form_factors', self.form_factors)]

270 for k, qpt in enumerate(self.q_points):

271 s += [fmt.format(f'q-point{k}', qpt)]

272 return '\n'.join(s)

Coverage for mchammer/observers/structure_factor_observer.py: 100%

102 statements