Coverage for mchammer/ensembles/hybrid_ensemble.py: 96%
138 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1import random
2from typing import Dict, List, Any
4from math import isclose
6from ase import Atoms
7from ase.units import kB
8from ase.data import atomic_numbers, chemical_symbols
10from .. import DataContainer
11from ..calculators.base_calculator import BaseCalculator
12from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble
13from .vcsgc_ensemble import get_phis
14from .semi_grand_canonical_ensemble import get_chemical_potentials
17class HybridEnsemble(ThermodynamicBaseEnsemble):
18 r"""
19 Instances of this class allows one to combine multiple ensembles.
20 In particular, a dictionary should be provided for each ensemble,
21 which must include the type (:attr:`ensemble`) as well as the index of
22 the sublattice (:attr:`sublattice_index`). In addition, it is possible
23 to provide a list of allowed symbols (:attr:`allowed_symbols`), which
24 must represent a subset of the elements that can occupy the sites
25 on the specified sublattice. Note that additional arguments are
26 required for the SGC and VCSGC ensembles, namely chemical
27 potentials (:attr:`chemical_potentials`) for the former and constraint
28 parameters (:attr:`phis` and :attr:`kappa`) for the latter. For more detailed
29 information regarding the different ensembles, please see
30 :class:`CanonicalEnsemble <mchammer.ensembles.CanonicalEnsemble>`,
31 :class:`SemiGrandCanonicalEnsemble
32 <mchammer.ensembles.SemiGrandCanonicalEnsemble>`, and
33 :class:`VCSGCEnsemble <mchammer.ensembles.VCSGCEnsemble>`.
35 This class is particularly useful for effectively sampling complex
36 multi-component systems with several active sublattices, in which
37 case different ensembles can be defined for each of the latter.
38 The fact that it is possible to set the allowed chemical symbols
39 means that one can vary the concentrations of a few selected
40 species, with the help of a VCSGC or semi-grand canonical
41 ensemble, while still allowing swaps between any two sites, using
42 a canonical ensemble (see also the below example). It is advisable
43 to carefully consider how to define the ensemble probabilities. By
44 default the ensembles are weighted by the sizes of the
45 corresponding sublattices, which should give suitable
46 probabilities in most cases. As is shown in the example below, it
47 might be prudent to provide different values if allowed symbols
48 are provided as well as for cases where there are several
49 ensembles that are active on different sublattices.
51 Parameters
52 ----------
54 structure
55 Atomic configuration to be used in the Monte Carlo simulation;
56 also defines the initial occupation vector.
57 calculator
58 Calculator to be used for calculating the potential changes
59 that enter the evaluation of the Metropolis criterion.
60 temperature
61 Temperature :math:`T` in appropriate units, commonly Kelvin.
62 ensemble_specs
63 A list of dictionaries, which should contain the following items:
65 * :attr:`ensemble`: Can be either ``'vcsgc'``, ``'semi-grand'``
66 or ``'canonical'``, lowercase and uppercase letters or any
67 combination thereof are accepted (required).
68 * :attr:`sublattice_index`: Index for the sublattice of
69 interest (required).
70 * :attr:`allowed_symbols`: List of allowed chemical symbols.
71 Default: read from :class:`ClusterSpace`.
72 * :attr:`chemical_potentials`: Dictionary of chemical
73 potentials for each species :math:`\mu_i`.
74 The key denotes the species, the value
75 specifies the chemical potential in units that are
76 consistent with the underlying cluster expansion.
77 Only applicable and required for SGC ensembles.
78 * :attr:`phis`: Dictionary with average constraint parameters
79 :math:`\phi_i`. The key denotes the species. For a
80 :math:`N`-component sublattice, there should be :math:`N - 1`
81 different :math:`\phi_i` (referred to as :math:`\bar{\phi}` in [SadErh12]_).
82 Only applicable and required for VCSGC ensembles.
83 * :attr:`kappa`: Parameter that constrains the variance of the
84 concentration (referred to as :math:`\bar{\kappa}` in [SadErh12]_).
85 Only applicable and required for VCSGC ensembles.
87 probabilities
88 Probabilities for choosing a
89 particular ensemble with the same length as ensemble specs.
90 If left unspecified the probabilties are scaled to match the
91 sizes of the associated sublattices.
92 boltzmann_constant
93 Boltzmann constant :math:`k_B` in appropriate units, i.e.,
94 units that are consistent with the underlying cluster
95 expansion and the temperature units. Default: eV/K.
96 user_tag
97 Human-readable tag for ensemble. Default: ``None``.
98 random_seed
99 Seed for the random number generator used in the Monte Carlo simulation.
100 dc_filename
101 Name of file the data container associated with the
102 ensemble will be written to. If the file
103 exists it will be read, the data container will be appended,
104 and the file will be updated/overwritten.
105 data_container_write_period
106 Period in units of seconds at which the data container is
107 written to file. Writing periodically to file provides both
108 a way to examine the progress of the simulation and to
109 back up the data. Default: 600 s.
110 ensemble_data_write_interval
111 Interval at which data is written to the data container.
112 This includes for example the current value of the
113 calculator (i.e., usually the energy) as well as ensembles
114 specific fields such as temperature or the number of atoms
115 of different species.
116 Default: Number of sites in the :attr:`structure`.
117 trajectory_write_interval
118 Interval at which the current occupation vector of the
119 atomic configuration is written to the data container.
120 Default: Number of sites in the :attr:`structure`.
122 Example
123 -------
124 The following snippet illustrates how to carry out a simple Monte Carlo
125 simulation using a combination of one canonical and one VCSGC ensemble.
126 Specifically, the concentration of one species (Au) is kept constant
127 while the others (Ag and Pd) are varied, while swaps are still allowed.
128 Here, the parameters of the cluster expansion are set to emulate a simple
129 Ising model in order to obtain an example that can be run without
130 modification. In practice, one should of course use a proper cluster
131 expansion::
133 >>> from ase.build import bulk
134 >>> from icet import ClusterExpansion, ClusterSpace
135 >>> from mchammer.calculators import ClusterExpansionCalculator
137 >>> # prepare cluster expansion
138 >>> # the setup emulates a second nearest-neighbor (NN) Ising model
139 >>> # (zerolet and singlet ECIs are zero; only first and second neighbor
140 >>> # pairs are included)
141 >>> prim = bulk('Au')
142 >>> cs = ClusterSpace(prim, cutoffs=[4.3],
143 ... chemical_symbols=['Ag', 'Au', 'Pd'])
144 >>> ce = ClusterExpansion(
145 ... cs, [0, 0, 0, 0.1, 0.1, 0.1, -0.02, -0.02, -0.02])
147 >>> # define structure object
148 >>> structure = prim.repeat(3)
149 >>> for i, atom in enumerate(structure):
150 >>> if i % 2 == 0:
151 >>> atom.symbol = 'Ag'
152 >>> elif i % 3 == 0:
153 >>> atom.symbol = 'Pd'
155 >>> # the default probabilities for this case would be [0.5, 0.5], but
156 >>> # since the VCSGC ensemble only performs flips on a subset of all
157 >>> # sites on the sublattice, namely those originally occupied by Ag
158 >>> # and Pd atoms, specific values will be provided
159 >>> weights = [len(structure),
160 ... len([s for s in structure.get_chemical_symbols() if s != 'Au'])]
161 >>> norm = sum(weights)
162 >>> probabilities = [w / norm for w in weights]
164 >>> # set up and run MC simulation
165 >>> calc = ClusterExpansionCalculator(structure, ce)
166 >>> ensemble_specs = [
167 ... {'ensemble': 'canonical', 'sublattice_index': 0},
168 ... {'ensemble': 'vcsgc', 'sublattice_index': 0,
169 ... 'phis': {'Ag': -0.2}, 'kappa': 200,
170 ... 'allowed_symbols':['Ag', 'Pd']}]
171 >>> mc = HybridEnsemble(structure=structure, calculator=calc,
172 ... ensemble_specs=ensemble_specs,
173 ... temperature=600, probabilities=probabilities,
174 ... dc_filename='myrun_hybrid.dc')
175 >>> mc.run(100) # carry out 100 trial steps
176 """
178 def __init__(self,
179 structure: Atoms,
180 calculator: BaseCalculator,
181 temperature: float,
182 ensemble_specs: List[Dict],
183 probabilities: List[float] = None,
184 boltzmann_constant: float = kB,
185 user_tag: str = None,
186 random_seed: int = None,
187 dc_filename: str = None,
188 data_container: str = None,
189 data_container_write_period: float = 600,
190 ensemble_data_write_interval: int = None,
191 trajectory_write_interval: int = None) -> None:
193 # define available ensembles
194 self._ensemble_trial_steps = dict([
195 ('canonical', self.do_canonical_swap),
196 ('semi-grand', self.do_sgc_flip),
197 ('vcsgc', self.do_vcsgc_flip),
198 ])
200 self._ensemble_parameters = dict(temperature=temperature)
202 self._trial_steps_per_ensemble = {'ensemble_{}'.format(i): 0 for i in
203 range(len(ensemble_specs))}
205 # process the list of ensembles and parameters
206 self._process_ensemble_specs(ensemble_specs)
208 super().__init__(
209 structure=structure,
210 calculator=calculator,
211 user_tag=user_tag,
212 random_seed=random_seed,
213 boltzmann_constant=boltzmann_constant,
214 dc_filename=dc_filename,
215 data_container=data_container,
216 data_container_class=DataContainer,
217 data_container_write_period=data_container_write_period,
218 ensemble_data_write_interval=ensemble_data_write_interval,
219 trajectory_write_interval=trajectory_write_interval)
221 # postprocess the list of ensembles and parameters
222 self._postprocess_ensemble_args()
224 # set the probabilities
225 self._process_probabilities(probabilities)
227 @property
228 def temperature(self) -> float:
229 """ Current temperature. """
230 return self._ensemble_parameters['temperature']
232 @property
233 def probabilities(self) -> List[float]:
234 """ Ensemble propabilities. """
235 return self._probabilities
237 @property
238 def trial_steps_per_ensemble(self) -> Dict[str, int]:
239 """ Number of Monte Carlo trial steps for each ensemble. """
240 return self._trial_steps_per_ensemble
242 def _process_ensemble_specs(self, ensemble_specs: List[Dict]) -> None:
243 r"""Process the list of ensembles and parameters.
245 Parameters
246 ----------
247 ensemble_specs: List[Dict]
248 A list of dictionaries, which should contain the following items:
249 * 'ensemble', which could be either 'vcsgc'; 'semi-grand' or 'canonical', lowercase
250 and upercase letters or any combination thereof are accepted
251 * 'sublattice_index', index for the sublattice of interest
252 * 'allowed_symbols', list of allowed chemical symbols
253 * 'chemical_potentials', a dictionary of chemical potentials for each species
254 :math:`\mu_i`; the key denotes the species, the value specifies the chemical
255 potential in units that are consistent with the underlying cluster expansion
256 * 'phis ', dictionary with average constraint parameters :math:`\phi_i`; the key
257 denotes the species; for a N-component sublattice, there should be N-1
258 different `\phi_i`
259 * 'kappa', parameter that constrains the variance of the concentration
260 """
262 ensemble_args = []
264 for ind, ensemble_spec in enumerate(ensemble_specs):
266 ensemble_arg: Dict[str, Any] = {}
267 tag = f'ensemble_{ind}'
268 ensemble_arg['tag'] = tag
270 # check the ensemble name
271 if 'ensemble' not in ensemble_spec:
272 raise ValueError(
273 f"The dictionary {ensemble_spec} lacks the required key 'ensemble'")
274 ensemble = ensemble_spec['ensemble'].lower()
275 if ensemble not in self._ensemble_trial_steps.keys():
276 msg = ['Ensemble not available']
277 msg += ['Please choose one of the following:']
278 for key in self._ensemble_trial_steps.keys():
279 msg += [' * ' + key]
280 raise ValueError('\n'.join(msg))
281 ensemble_arg['ensemble'] = ensemble
282 self._ensemble_parameters[tag] = ensemble
284 # check that all required keys, and no unknown keys, are present
285 keys = ['ensemble', 'sublattice_index', 'allowed_symbols']
286 if ensemble == 'semi-grand':
287 keys = ['chemical_potentials'] + keys
288 elif ensemble == 'vcsgc':
289 keys = ['phis', 'kappa'] + keys
290 for key in keys[:-1]:
291 if key not in ensemble_spec:
292 raise ValueError(f"The dictionary {ensemble_spec} lacks the key '{key}',"
293 f' which is required for {ensemble} ensembles')
294 for key in ensemble_spec.keys():
295 if key not in keys:
296 raise ValueError(f"Unknown key '{key}', for a {ensemble} ensemble,"
297 f' in the dictionary {ensemble_spec}')
299 # record the sublattice index
300 ensemble_arg['sublattice_index'] = ensemble_spec['sublattice_index']
302 # process chemical potentials
303 if 'chemical_potentials' in ensemble_spec:
304 chemical_potentials = get_chemical_potentials(ensemble_spec['chemical_potentials'])
305 ensemble_arg['chemical_potentials'] = chemical_potentials
306 for atnum, chempot in chemical_potentials.items():
307 mu_sym = '{}_mu_{}'.format(tag, chemical_symbols[atnum])
308 self._ensemble_parameters[mu_sym] = chempot
310 # process phis
311 if 'phis' in ensemble_spec:
312 phis = get_phis(ensemble_spec['phis'])
313 ensemble_arg['phis'] = phis
314 for sym, phi in phis.items():
315 if isinstance(sym, str): 315 ↛ 316line 315 didn't jump to line 316, because the condition on line 315 was never true
316 chemical_symbol = sym
317 else:
318 chemical_symbol = chemical_symbols[sym]
319 phi_sym = '{}_phi_{}'.format(tag, chemical_symbol)
320 self._ensemble_parameters[phi_sym] = phi
322 # process kappa
323 if 'kappa' in ensemble_spec:
324 ensemble_arg['kappa'] = ensemble_spec['kappa']
325 self._ensemble_parameters['{}_kappa'.format(tag)] = ensemble_spec['kappa']
327 # record the allowed chemical symbols
328 if 'allowed_symbols' in ensemble_spec:
329 ensemble_arg['allowed_symbols'] = ensemble_spec['allowed_symbols']
331 ensemble_args.append(ensemble_arg)
333 self._ensemble_args = ensemble_args
335 def _postprocess_ensemble_args(self):
336 """Process the list of dictionaries with ensemble specific parameters.
337 """
339 for i in range(len(self._ensemble_args)):
341 # check the sublattice index
342 self._check_sublattice_index(self._ensemble_args[i]['sublattice_index'])
344 # extract the allowed species
345 if 'allowed_symbols' in self._ensemble_args[i]:
346 self._ensemble_args[i]['allowed_species'] =\
347 self._extract_allowed_species(self._ensemble_args[i]['allowed_symbols'],
348 self._ensemble_args[i]['sublattice_index'])
349 del self._ensemble_args[i]['allowed_symbols']
350 else:
351 self._ensemble_args[i]['allowed_species'] = None
353 if self._ensemble_args[i]['ensemble'] == 'vcsgc':
354 # Check that each sublattice has N - 1 phis
355 count_specified_elements = 0
356 if self._ensemble_args[i]['allowed_species'] is None: 356 ↛ 360line 356 didn't jump to line 360, because the condition on line 356 was never false
357 allowed_species =\
358 self.sublattices[self._ensemble_args[i]['sublattice_index']].atomic_numbers
359 else:
360 allowed_species = self._ensemble_args[i]['allowed_species']
361 for number in allowed_species:
362 if number in self._ensemble_args[i]['phis'].keys():
363 count_specified_elements += 1
364 if count_specified_elements != len(allowed_species) - 1:
365 raise ValueError('phis must be set for N - 1 elements on a sublattice with'
366 ' N elements')
368 def _check_sublattice_index(self, sublattice_index: int):
369 """Check the :attr:`sublattice_index` item in the :attr:`ensemble_spec` dictionary.
371 Parameters
372 ----------
373 sublattice_index
374 Specific sublattice to consider provided as as an index or a symbol.
375 """
377 if not isinstance(sublattice_index, int):
378 raise TypeError("'sublattice_index' must be an integer, not"
379 f' {format(type(sublattice_index))}')
381 # check that the sublattice exists
382 if sublattice_index not in range(len(self.sublattices)):
383 raise ValueError('There is no sublattice with index {}'.format(sublattice_index))
385 # check that the sublattice is active
386 if len(self.sublattices[sublattice_index].chemical_symbols) == 1:
387 raise ValueError('The sublattice {} is inactive'.format(sublattice_index))
389 def _extract_allowed_species(
390 self,
391 allowed_symbols: List[str],
392 sublattice_index: int,
393 ) -> List[int]:
394 """Check and extract the allowed species from the
395 :attr:`allowed_symbols` in the :attr:`ensemble_spec`
396 dictionary.
398 Parameters
399 ----------
400 allowed_symbols
401 List of allowed chemical symbols.
402 sublattice_index
403 Index for the relevant sublattice.
405 """
407 if not isinstance(allowed_symbols, list) or not all(
408 [isinstance(i, str) for i in allowed_symbols]):
409 raise TypeError(
410 f"'allowed_symbols' must be a List[str], not {type(allowed_symbols)}")
411 for symbol in allowed_symbols:
412 if symbol not in self.sublattices[sublattice_index].chemical_symbols:
413 raise ValueError('The species {} is not allowed on sublattice'
414 ' {}'.format(symbol, sublattice_index))
416 return [atomic_numbers[s] for s in allowed_symbols]
418 def _process_probabilities(self, probabilities: List[float]):
419 """Process the list of probabilities.
421 Parameters
422 ----------
423 probabilities
424 Probabilities for choosing a particular ensemble with the same
425 length as :attr:`self._ensemble_args`.
426 """
428 if probabilities is None:
429 # use the sizes of the associated sublattices when calculating the ensemble
430 # probabilities
431 weights = [len(self.sublattices[ensemble_arg['sublattice_index']].indices) for
432 ensemble_arg in self._ensemble_args]
433 norm = sum(weights)
434 probabilities = [w / norm for w in weights]
435 else:
436 if len(probabilities) != len(self._ensemble_args):
437 raise ValueError('The number of probabilities must be match the number of'
438 ' ensembles')
440 if not isclose(sum(probabilities), 1.0):
441 raise ValueError('The sum of all probabilities must be equal to 1')
443 self._probabilities = probabilities
445 def _do_trial_step(self):
446 """ Carries out one Monte Carlo trial step. """
447 # randomly pick an ensemble
448 ensemble_arg = random.choices(self._ensemble_args, weights=self._probabilities)[0]
450 # count number of trial steps for each ensemble
451 self._trial_steps_per_ensemble[ensemble_arg['tag']] += 1
453 if ensemble_arg['ensemble'] == 'canonical' and not self.configuration.is_swap_possible( 453 ↛ 455line 453 didn't jump to line 455, because the condition on line 453 was never true
454 ensemble_arg['sublattice_index'], ensemble_arg['allowed_species']):
455 return 0
456 else:
457 arguments = {key: val for key, val in ensemble_arg.items() if key not in
458 ['ensemble', 'tag']}
459 return self._ensemble_trial_steps[ensemble_arg['ensemble']](**arguments)
461 def _get_ensemble_data(self) -> Dict:
462 """
463 Returns a dict with the default data of the ensemble. This includes
464 atom counts and free energy derivative.
465 """
466 data = super()._get_ensemble_data()
468 ensemble_types = [e['ensemble'] for e in self._ensemble_args]
470 # free energy derivative
471 if 'vcsgc' in ensemble_types: 471 ↛ 479line 471 didn't jump to line 479, because the condition on line 471 was never false
472 for ensemble_arg in self._ensemble_args:
473 if 'vcsgc' == ensemble_arg['ensemble']:
474 data.update(self._get_vcsgc_free_energy_derivatives(
475 ensemble_arg['phis'], ensemble_arg['kappa'],
476 ensemble_arg['sublattice_index']))
478 # species counts
479 if any([e in ensemble_types for e in ['vcsgc', 'semi-grand']]): 479 ↛ 482line 479 didn't jump to line 482, because the condition on line 479 was never false
480 data.update(self._get_species_counts())
482 return data