Coverage for mchammer/ensembles/wang_landau_ensemble.py: 93%
244 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1"""Definition of the Wang-Landau algorithm class."""
3import random
5from collections import OrderedDict
6from typing import Dict, List, Optional, Tuple, Union
8import numpy as np
10from ase import Atoms
12from .. import WangLandauDataContainer
13from ..calculators.base_calculator import BaseCalculator
14from .thermodynamic_base_ensemble import BaseEnsemble
15from icet.input_output.logging_tools import logger
17logger = logger.getChild('wang_landau_ensemble')
20class WangLandauEnsemble(BaseEnsemble):
21 r"""Instances of this class allow one to sample a system using the
22 Wang-Landau (WL) algorithm, see Phys. Rev. Lett. **86**, 2050
23 (2001) [WanLan01a]_. The WL algorithm enables one to acquire the
24 density of states (DOS) as a function of energy, from which one
25 can readily calculate many thermodynamic observables as a function
26 of temperature. To this end, the WL algorithm accumulates both the
27 microcanonical entropy :math:`S(E)` and a histogram :math:`H(E)`
28 on an energy grid with a predefined spacing (:attr:`energy_spacing`).
30 The algorithm is initialized as follows.
32 #. Generate an initial configuration.
33 #. Initialize counters for the microcanonical entropy
34 :math:`S(E)` and the histogram :math:`H(E)` to zero.
35 #. Set the fill factor :math:`f=1`.
37 It then proceeds as follows.
39 #. Propose a new configuration (see :attr:`trial_move`).
40 #. Accept or reject the new configuration with probability
42 .. math::
44 P = \min \{ 1, \, \exp [ S(E_\mathrm{new}) - S(E_\mathrm{cur}) ] \},
46 where :math:`E_\mathrm{cur}` and :math:`E_\mathrm{new}` are the
47 energies of the current and new configurations, respectively.
48 #. Update the microcanonical entropy :math:`S(E)\leftarrow S(E) + f`
49 and histogram :math:`H(E) \leftarrow H(E) + 1` where
50 :math:`E` is the energy of the system at the end of the move.
51 #. Check the flatness of the histogram :math:`H(E)`. If
52 :math:`H(E) > \chi \langle H(E)\rangle\,\forall E` reset the histogram
53 :math:`H(E) = 0` and reduce the fill factor :math:`f \leftarrow f / 2`.
54 The parameter :math:`\chi` is set via :attr:`flatness_limit`.
55 #. If :math:`f` is smaller than :attr:`fill_factor_limit` terminate
56 the loop, otherwise return to 1.
58 The microcanonical entropy :math:`S(E)` and the histogram along
59 with related information are written to the data container every
60 time :math:`f` is updated. Using the density :math:`\rho(E) = \exp S(E)`
61 one can then readily compute various thermodynamic quantities,
62 including, e.g., the average energy:
64 .. math::
66 \left<E\right> = \frac{\sum_E E \rho(E) \exp(-E / k_B T)}{
67 \sum_E \rho(E) \exp(-E / k_B T)}
69 Parameters
70 ----------
71 structure
72 Atomic configuration to be used in the Wang-Landau simulation;
73 also defines the initial occupation vector.
74 calculator
75 Calculator to be used for calculating potential changes.
76 trial_move
77 One can choose between two different trial moves for
78 generating new configurations. In a ``'swap'`` move two sites are
79 selected and their occupations are swapped. In a ``'flip'`` move
80 one site is selected and its occupation is flipped to a
81 different species. While ``'swap'`` moves conserve the
82 concentrations of the species in the system, ``'flip'`` moves
83 allow one in principle to sample the full composition space.
84 energy_spacing
85 Sets the bin size of the energy grid on which the microcanonical
86 entropy :math:`S(E)`, and thus the density :math:`\exp S(E)`, is
87 evaluated. The spacing should be small enough to capture the features
88 of the density of states. Too small values will, however, render the
89 convergence very tedious if not impossible.
90 energy_limit_left
91 Sets the lower limit of the energy range within which the
92 microcanonical entropy :math:`S(E)` will be sampled. By default
93 (``None``) no limit is imposed. Setting limits can be useful if only a
94 part of the density of states is required.
95 energy_limit_right
96 Sets the upper limit of the energy range within which the
97 microcanonical entropy :math:`S(E)` will be sampled. By default
98 (``None``) no limit is imposed. Setting limits can be useful if only a
99 part of the density of states is required.
100 fill_factor_limit
101 If the fill_factor :math:`f` falls below this value, the
102 WL sampling loop is terminated.
103 flatness_check_interval
104 For computational efficiency the flatness condition is only
105 evaluated every :attr:`flatness_check_interval`-th trial step. By
106 default (``None``) :attr:`flatness_check_interval` is set to 1000
107 times the number of sites in :attr:`structure`, i.e., 1000 Monte
108 Carlo sweeps.
109 flatness_limit
110 The histogram :math:`H(E)` is deemed sufficiently flat if
111 :math:`H(E) > \chi \left<H(E)\right>\,\forall
112 E`. :attr:`flatness_limit` sets the parameter :math:`\chi`.
113 window_search_penalty
114 If :attr:`energy_limit_left` and/or :attr:`energy_limit_right` have been
115 set, a modified acceptance probability,
116 :math:`P=\min\{1,\,\exp[C_\mathrm{WSP}(d_\mathrm{new}-
117 d_\mathrm{cur})]\}`, will be used until a configuration is
118 found within the interval of interest. This parameter,
119 specifically, corresponds to :math:`C_\mathrm{WSP}`, which
120 controls how strongly moves that lead to an increase in the
121 distance, i.e. difference in energy divided by the energy
122 spacing, to the energy window (:math:`d_\mathrm{new}>
123 d_\mathrm{cur}`) should be penalized. A higher value leads
124 to a lower acceptance probability for such moves.
125 user_tag
126 Human-readable tag for ensemble. Default: ``None``.
127 dc_filename
128 Name of file the data container associated with the ensemble
129 will be written to. If the file exists it will be read, the
130 data container will be appended, and the file will be
131 updated/overwritten.
132 random_seed
133 Seed for the random number generator used in the Monte Carlo
134 simulation.
135 ensemble_data_write_interval
136 Interval at which data is written to the data container. This
137 includes for example the current value of the calculator
138 (i.e., usually the energy) as well as ensemble specific fields
139 such as temperature or the number of atoms of different species.
140 Default: Number of sites in the :attr:`structure`.
141 data_container_write_period
142 Period in units of seconds at which the data container is
143 written to file. Writing periodically to file provides both
144 a way to examine the progress of the simulation and to back up
145 the data. Default: 600 s.
146 trajectory_write_interval
147 Interval at which the current occupation vector of the atomic
148 configuration is written to the data container.
149 Default: Number of sites in the :attr:`structure`.
150 sublattice_probabilities
151 Probability for picking a sublattice when doing a random swap.
152 The list must contain as many elements as there are sublattices
153 and it needs to sum up to 1.
155 Example
156 -------
157 The following snippet illustrates how to carry out a Wang-Landau
158 simulation. For the purpose of demonstration, the parameters of
159 the cluster expansion are set to obtain a two-dimensional square
160 Ising model, one of the systems studied in the original work by
161 Wang and Landau::
163 >>> from ase import Atoms
164 >>> from icet import ClusterExpansion, ClusterSpace
165 >>> from mchammer.calculators import ClusterExpansionCalculator
166 >>> from mchammer.ensembles import WangLandauEnsemble
168 >>> # prepare cluster expansion
169 >>> prim = Atoms('Au', positions=[[0, 0, 0]], cell=[1, 1, 10], pbc=True)
170 >>> cs = ClusterSpace(prim, cutoffs=[1.1], chemical_symbols=['Ag', 'Au'])
171 >>> ce = ClusterExpansion(cs, [0, 0, 2])
173 >>> # prepare initial configuration
174 >>> structure = prim.repeat((4, 4, 1))
175 >>> for k in range(8):
176 ... structure[k].symbol = 'Ag'
178 >>> # set up and run Wang-Landau simulation
179 >>> calculator = ClusterExpansionCalculator(structure, ce)
180 >>> mc = WangLandauEnsemble(structure=structure,
181 ... calculator=calculator,
182 ... energy_spacing=1,
183 ... dc_filename='ising_2d_run.dc')
184 >>> mc.run(number_of_trial_steps=len(structure)*100)
185 >>> # Note: in practice one requires many more steps
187 """
189 def __init__(self,
190 structure: Atoms,
191 calculator: BaseCalculator,
192 energy_spacing: float,
193 energy_limit_left: float = None,
194 energy_limit_right: float = None,
195 trial_move: str = 'swap',
196 fill_factor_limit: float = 1e-6,
197 flatness_check_interval: int = None,
198 flatness_limit: float = 0.8,
199 window_search_penalty: float = 2.0,
200 user_tag: str = None,
201 dc_filename: str = None,
202 data_container: str = None,
203 random_seed: int = None,
204 data_container_write_period: float = 600,
205 ensemble_data_write_interval: int = None,
206 trajectory_write_interval: int = None,
207 sublattice_probabilities: List[float] = None) -> None:
209 # set trial move
210 if trial_move == 'swap':
211 self.do_move = self._do_swap
212 self._get_sublattice_probabilities = self._get_swap_sublattice_probabilities
213 elif trial_move == 'flip':
214 self.do_move = self._do_flip
215 self._get_sublattice_probabilities = self._get_flip_sublattice_probabilities
216 else:
217 raise ValueError('Invalid value for trial_move: {}.'
218 ' Must be either "swap" or "flip".'.format(trial_move))
220 # set default values that are system dependent
221 if flatness_check_interval is None:
222 flatness_check_interval = len(structure) * 1000
224 # parameters pertaining to construction of entropy and histogram
225 self._energy_spacing = energy_spacing
226 self._fill_factor_limit = fill_factor_limit
227 self._flatness_check_interval = flatness_check_interval
228 self._flatness_limit = flatness_limit
230 # energy window
231 self._window_search_penalty = window_search_penalty
232 self._bin_left = self._get_bin_index(energy_limit_left)
233 self._bin_right = self._get_bin_index(energy_limit_right)
234 if self._bin_left is not None and \
235 self._bin_right is not None and self._bin_left >= self._bin_right:
236 raise ValueError('Invalid energy window: left boundary ({}, {}) must be'
237 ' smaller than right boundary ({}, {})'
238 .format(energy_limit_left, self._bin_left,
239 energy_limit_right, self._bin_right))
241 # ensemble parameters
242 self._ensemble_parameters = {}
243 self._ensemble_parameters['energy_spacing'] = energy_spacing
244 self._ensemble_parameters['trial_move'] = trial_move
245 self._ensemble_parameters['energy_limit_left'] = energy_limit_left
246 self._ensemble_parameters['energy_limit_right'] = energy_limit_right
247 # The following parameters are _intentionally excluded_ from
248 # the ensemble_parameters dict as it would prevent users from
249 # changing their values between restarts. The latter is advantageous
250 # as these runs can require restarts and possibly parameter adjustments
251 # to achieve convergence.
252 # * fill_factor_limit
253 # * flatness_check_interval
254 # * flatness_limit
255 # * entropy_write_frequency
256 # * window_search_penalty
258 # add species count to ensemble parameters
259 symbols = set([symbol for sub in calculator.sublattices
260 for symbol in sub.chemical_symbols])
261 for symbol in symbols:
262 key = 'n_atoms_{}'.format(symbol)
263 count = structure.get_chemical_symbols().count(symbol)
264 self._ensemble_parameters[key] = count
266 # set the convergence, which may be updated in case of a restart
267 self._converged: bool = None
269 # the constructor of the parent classes must be called *after*
270 # the ensemble_parameters dict has been populated
271 super().__init__(
272 structure=structure,
273 calculator=calculator,
274 user_tag=user_tag,
275 random_seed=random_seed,
276 dc_filename=dc_filename,
277 data_container=data_container,
278 data_container_class=WangLandauDataContainer,
279 data_container_write_period=data_container_write_period,
280 ensemble_data_write_interval=ensemble_data_write_interval,
281 trajectory_write_interval=trajectory_write_interval)
283 # handle probabilities for swaps on different sublattices
284 if sublattice_probabilities is None:
285 self._sublattice_probabilities = self._get_sublattice_probabilities()
286 else:
287 self._sublattice_probabilities = sublattice_probabilities
289 # initialize Wang-Landau algorithm; in the case of a restart
290 # these quantities are read from the data container file; the
291 # if-conditions prevent these values from being overwritten
292 self._potential = self.calculator.calculate_total(
293 occupations=self.configuration.occupations)
294 self._reached_energy_window = self._inside_energy_window(
295 self._get_bin_index(self._potential))
296 if not hasattr(self, '_fill_factor'):
297 self._fill_factor = 1.0
298 if not hasattr(self, '_fill_factor_history'):
299 if self._reached_energy_window:
300 self._fill_factor_history = {self.step: self._fill_factor}
301 else:
302 self._fill_factor_history = {}
303 if not hasattr(self, '_entropy_history'):
304 self._entropy_history = {}
305 if not hasattr(self, '_histogram'):
306 self._histogram: Dict[int, int] = {}
307 if not hasattr(self, '_entropy'):
308 self._entropy: Dict[int, float] = {}
310 @property
311 def fill_factor(self) -> float:
312 """ current value of fill factor """
313 return self._fill_factor
315 @property
316 def fill_factor_history(self) -> Dict[int, float]:
317 """evolution of the fill factor in the Wang-Landau algorithm (key=MC
318 trial step, value=fill factor)
319 """
320 return self._fill_factor_history
322 @property
323 def converged(self) -> Optional[bool]:
324 """ True if convergence has been achieved """
325 return self._converged
327 @property
328 def flatness_limit(self) -> float:
329 r"""The histogram :math:`H(E)` is deemed sufficiently flat if
330 :math:`H(E) > \chi \left<H(E)\right>\,\forall
331 E` where :attr:`flatness_limit` sets the parameter :math:`\chi`.
332 """
333 return self._flatness_limit
335 @flatness_limit.setter
336 def flatness_limit(self, new_value) -> None:
337 self._flatness_limit = new_value
338 self._converged = None
340 @property
341 def fill_factor_limit(self) -> float:
342 """ If the fill factor :math:`f` falls below this value, the
343 Wang-Landau sampling is terminated. """
344 return self._fill_factor_limit
346 @fill_factor_limit.setter
347 def fill_factor_limit(self, new_value) -> None:
348 self._fill_factor_limit = new_value
349 self._converged = None
351 @property
352 def flatness_check_interval(self) -> int:
353 """ Number of MC trial steps between checking the flatness condition. """
354 return self._flatness_check_interval
356 @flatness_check_interval.setter
357 def flatness_check_interval(self, new_value: int) -> None:
358 self._flatness_check_interval = new_value
360 def run(self, number_of_trial_steps: int):
361 """
362 Samples the ensemble for the given number of trial steps.
364 Parameters
365 ----------
366 number_of_trial_steps
367 Maximum number of MC trial steps to run in total. The
368 run will terminate earlier if :attr:`fill_factor_limit` is reached.
369 reset_step
370 If ``True`` the MC trial step counter and the data container will
371 be reset to zero and empty, respectively.
373 Raises
374 ------
375 TypeError
376 If :attr:`number_of_trial_steps` is not an ``int``.
377 """
378 if self.converged:
379 logger.warning('Convergence has already been reached.')
380 else:
381 super().run(number_of_trial_steps)
383 def _terminate_sampling(self) -> bool:
384 """Returns ``True`` if the Wang-Landau algorithm has converged. This is
385 used in the :func:`run` method implemented in :class:`BaseEnsemble` to
386 evaluate whether the sampling loop should be terminated.
387 """
388 # N.B.: self._converged can be None
389 if self._converged is not None:
390 return self._converged
391 else:
392 return False
394 def _restart_ensemble(self):
395 """Restarts ensemble using the last state saved in the data container
396 file. Note that this method does _not_ use the :attr:`last_state` property of
397 the data container but rather uses the last data written to the data frame.
398 """
399 super()._restart_ensemble()
400 self._fill_factor = self.data_container._last_state['fill_factor']
401 self._fill_factor_history = self.data_container._last_state['fill_factor_history']
402 self._entropy_history = self.data_container._last_state['entropy_history']
403 self._histogram = self.data_container._last_state['histogram']
404 self._entropy = self.data_container._last_state['entropy']
405 histogram = np.array(list(self._histogram.values()))
406 limit = self._flatness_limit * np.average(histogram)
407 self._converged = (self._fill_factor <= self._fill_factor_limit
408 ) & np.all(histogram >= limit)
410 def write_data_container(self, outfile: Union[str, bytes]):
411 """Updates the last state of the Wang-Landau simulation and
412 writes the data container to file.
414 Parameters
415 ----------
416 outfile
417 File to which to write.
418 """
419 self._data_container._update_last_state(
420 last_step=self.step,
421 occupations=self.configuration.occupations.tolist(),
422 accepted_trials=self._accepted_trials,
423 random_state=random.getstate(),
424 fill_factor=self._fill_factor,
425 fill_factor_history=self._fill_factor_history,
426 entropy_history=self._entropy_history,
427 histogram=OrderedDict(sorted(self._histogram.items())),
428 entropy=OrderedDict(sorted(self._entropy.items())))
429 self.data_container.write(outfile)
431 def _acceptance_condition(self, potential_diff: float) -> bool:
432 """Evaluates Metropolis acceptance criterion.
434 Parameters
435 ----------
436 potential_diff
437 Change in the thermodynamic potential associated
438 with the trial step.
439 """
441 # acceptance/rejection step
442 bin_old = self._get_bin_index(self._potential)
443 bin_new = self._get_bin_index(self._potential + potential_diff)
444 bin_cur = bin_old
445 if self._allow_move(bin_cur, bin_new):
446 S_cur = self._entropy.get(bin_cur, 0)
447 S_new = self._entropy.get(bin_new, 0)
448 delta = np.exp(S_cur - S_new)
449 if delta >= 1 or delta >= self._next_random_number():
450 accept = True
451 self._potential += potential_diff
452 bin_cur = bin_new
453 else:
454 accept = False
455 else:
456 accept = False
458 if not self._reached_energy_window:
459 # check whether the target energy window has been reached
460 self._reached_energy_window = self._inside_energy_window(bin_cur)
461 # if the target window has been reached remove unused bins
462 # from histogram and entropy counters
463 if self._reached_energy_window:
464 self._fill_factor_history[self.step] = self._fill_factor
465 # flush data from data container except for initial step
466 self._data_container._data_list = [self._data_container._data_list[0]]
467 self._entropy = {k: self._entropy[k]
468 for k in self._entropy if self._inside_energy_window(k)}
469 self._histogram = {k: self._histogram[k]
470 for k in self._histogram if self._inside_energy_window(k)}
471 else:
472 # then reconsider accept/reject based on whether we
473 # approached the window or not
474 dist_new = np.inf
475 dist_old = np.inf
476 if self._bin_left is not None: 476 ↛ 479line 476 didn't jump to line 479, because the condition on line 476 was never false
477 dist_new = min(dist_new, abs(bin_new - self._bin_left))
478 dist_old = min(dist_old, abs(bin_old - self._bin_left))
479 if self._bin_right is not None:
480 dist_new = min(dist_new, abs(bin_new - self._bin_right))
481 dist_old = min(dist_old, abs(bin_old - self._bin_right))
482 assert dist_new < np.inf and dist_old < np.inf
483 exp_dist = np.exp((dist_old - dist_new) * self._window_search_penalty)
484 if exp_dist >= 1 or exp_dist >= self._next_random_number():
485 # should be accepted
486 if not accept:
487 # reset potential
488 self._potential += potential_diff
489 bin_cur = bin_new
490 accept = True
491 else:
492 # should be rejected
493 if accept: 493 ↛ 496line 493 didn't jump to line 496, because the condition on line 493 was never false
494 # reset potential
495 self._potential -= potential_diff
496 bin_cur = bin_old
497 accept = False
499 # update histograms and entropy counters
500 self._update_entropy(bin_cur)
502 return accept
504 def _update_entropy(self, bin_cur: int) -> None:
505 """Updates counters for histogram and entropy, checks histogram
506 flatness, and updates fill factor if indicated.
507 """
509 # update histogram and entropy
510 self._entropy[bin_cur] = self._entropy.get(bin_cur, 0) + self._fill_factor
511 self._histogram[bin_cur] = self._histogram.get(bin_cur, 0) + 1
513 # check flatness of histogram
514 if self.step % self._flatness_check_interval == 0 and \
515 self.step > 0 and self._reached_energy_window:
517 # shift entropy counter in order to avoid overflow
518 entropy_ref = np.min(list(self._entropy.values()))
519 for k in self._entropy:
520 self._entropy[k] -= entropy_ref
522 # check whether the Wang-Landau algorithm has converged
523 histogram = np.array(list(self._histogram.values()))
524 limit = self._flatness_limit * np.average(histogram)
525 is_flat = np.all(histogram >= limit)
526 self._converged = (self._fill_factor <= self._fill_factor_limit) & is_flat
527 if is_flat and not self._converged:
528 # update fill factor
529 self._fill_factor /= 2
530 self._fill_factor_history[self.step] = self._fill_factor
531 # update entropy history
532 self._entropy_history[self.step] = OrderedDict(
533 sorted(self._entropy.items()))
534 # reset histogram
535 self._histogram = dict.fromkeys(self._histogram, 0)
537 def _get_bin_index(self, energy: float) -> Optional[int]:
538 """ Returns bin index for histogram and entropy dictionaries. """
539 if energy is None or np.isnan(energy):
540 return None
541 return int(round(energy / self._energy_spacing))
543 def _allow_move(self, bin_cur: Optional[int], bin_new: int) -> bool:
544 """Returns ``True`` if the current move is to be included in the
545 accumulation of histogram and entropy. This logic has been
546 moved into a separate function in order to enhance
547 readability.
548 """
549 if self._bin_left is None and self._bin_right is None:
550 # no limits on energy window
551 return True
552 if self._bin_left is not None:
553 if bin_cur < self._bin_left:
554 # not yet in window (left limit)
555 return True
556 if bin_new < self._bin_left:
557 # imposing left limit
558 return False
559 if self._bin_right is not None:
560 if bin_cur > self._bin_right:
561 # not yet in window (right limit)
562 return True
563 if bin_new > self._bin_right:
564 # imposing right limit
565 return False
566 return True
568 def _inside_energy_window(self, bin_k: int) -> bool:
569 """Returns ``True`` if :attr:`bin_k` is inside the energy window specified for
570 this simulation.
571 """
572 if self._bin_left is not None and bin_k < self._bin_left:
573 return False
574 if self._bin_right is not None and bin_k > self._bin_right:
575 return False
576 return True
578 def _do_trial_step(self):
579 """ Carries out one Monte Carlo trial step. """
580 sublattice_index = self.get_random_sublattice_index(self._sublattice_probabilities)
581 return self.do_move(sublattice_index=sublattice_index)
583 def _do_swap(self, sublattice_index: int, allowed_species: List[int] = None) -> int:
584 """Carries out a Monte Carlo trial that involves swapping the species
585 on two sites. This method has been copied from
586 :class:`ThermodynamicBaseEnsemble`.
588 Parameters
589 ---------
590 sublattice_index
591 Index of sublattice the swap will act on.
592 allowed_species
593 List of atomic numbers for allowed species.
595 Returns
596 -------
597 Returns 1 or 0 depending on if trial move was accepted or rejected.
598 """
599 sites, species = self.configuration.get_swapped_state(sublattice_index, allowed_species)
600 potential_diff = self._get_property_change(sites, species)
601 if self._acceptance_condition(potential_diff):
602 self.update_occupations(sites, species)
603 return 1
604 return 0
606 def _do_flip(self, sublattice_index: int, allowed_species: List[int] = None) -> int:
607 """Carries out one Monte Carlo trial step that involves flipping the
608 species on one site. This method has been adapted from
609 :class:`ThermodynamicBaseEnsemble`.
611 Parameters
612 ---------
613 sublattice_index
614 Index of sublattice the flip will act on.
615 allowed_species
616 List of atomic numbers for allowed species.
618 Returns
619 -------
620 Returns 1 or 0 depending on if trial move was accepted or rejected.
621 """
622 index, species = self.configuration.get_flip_state(sublattice_index, allowed_species)
623 potential_diff = self._get_property_change([index], [species])
624 if self._acceptance_condition(potential_diff):
625 self.update_occupations([index], [species])
626 return 1
627 return 0
629 def _get_swap_sublattice_probabilities(self) -> List[float]:
630 """Returns sublattice probabilities suitable for swaps. This method
631 has been copied without modification from :class:`ThermodynamicBaseEnsemble`.
632 """
633 sublattice_probabilities = []
634 for i, sl in enumerate(self.sublattices):
635 if self.configuration.is_swap_possible(i):
636 sublattice_probabilities.append(len(sl.indices))
637 else:
638 sublattice_probabilities.append(0)
639 norm = sum(sublattice_probabilities)
640 if norm == 0:
641 raise ValueError('No swaps are possible on any of the active sublattices.')
642 sublattice_probabilities = [p / norm for p in sublattice_probabilities]
643 return sublattice_probabilities
645 def _get_flip_sublattice_probabilities(self) -> List[float]:
646 """Returns the default sublattice probability which is based on the
647 sizes of a sublattice. This method has been copied without
648 modification from :class:`ThermodynamicBaseEnsemble`.
649 """
650 sublattice_probabilities = []
651 for _, sl in enumerate(self.sublattices):
652 if len(sl.chemical_symbols) > 1:
653 sublattice_probabilities.append(len(sl.indices))
654 else:
655 sublattice_probabilities.append(0)
656 norm = sum(sublattice_probabilities)
657 sublattice_probabilities = [p / norm for p in sublattice_probabilities]
658 return sublattice_probabilities
661def get_bins_for_parallel_simulations(n_bins: int,
662 energy_spacing: float,
663 minimum_energy: float,
664 maximum_energy: float,
665 overlap: int = 4,
666 bin_size_exponent: float = 1.0) -> List[Tuple[float, float]]:
667 """Generates a list of energy bins (lower and upper bound) suitable for
668 parallel Wang-Landau simulations. For the latter, the energy range is
669 split up into a several bins (:attr:`n_bins`). Each bin is then sampled in a
670 separate Wang-Landau simulation. Once the density of states in the
671 individual bins has been converged the total density of states can be
672 constructed by patching the segments back together. To this end, one
673 requires some over overlap between the segments (:attr:`overlap`).
675 The function returns a list of tuples. Each tuple provides the lower
676 (:attr:`energy_limit_left`) and upper (:attr:`energy_limit_right`) bound of one
677 bin, which are then to be used to set :attr:`energy_limit_left` and
678 :attr:`energy_limit_right` when initializing a :class:`WangLandauEnsemble`
679 instance.
681 Note
682 ----
683 The left-most/right-most bin has no lower/upper bound (set to ``None``).
685 Parameters
686 ----------
687 n_bins
688 Number of bins.
689 energy_spacing
690 Sets the bin size of the energy grid used by the Wang-Landau
691 simulation, see :class:`WangLandauEnsemble` for details.
692 minimum_energy
693 An estimate for the lowest energy to be encountered in this system.
694 maximum_energy
695 An estimate for the highest energy to be encountered in this system.
696 overlap
697 Amount of overlap between bins in units of :attr:`energy_spacing`.
698 bin_size_exponent
699 This parameter allows one to generate a non-uniform
700 distribution of bin sizes. If :attr:`bin_size_exponent` is smaller
701 than one bins at the lower and upper end of the energy range
702 (specified via :attr:`minimum_energy` and :attr:`maximum_energy`) will
703 be shrunk relative to the bins in the middle of the energy
704 range. In principle this can be used one to achieve a more
705 even distribution of computational load between the individual
706 Wang-Landau simulations.
708 Note
709 ----
710 This is an option for advanced users. Only use this keyword
711 if you know what you are doing.
712 """
714 limits = np.linspace(-1, 1, n_bins + 1)
715 limits = np.sign(limits) * np.abs(limits) ** bin_size_exponent
716 limits *= 0.5 * (maximum_energy - minimum_energy)
717 limits += 0.5 * (maximum_energy + minimum_energy)
718 limits[0], limits[-1] = None, None
720 bounds = []
721 for k, (energy_limit_left, energy_limit_right) in enumerate(zip(limits[:-1], limits[1:])):
722 if energy_limit_left is not None and energy_limit_right is not None and \
723 (energy_limit_right - energy_limit_left) / energy_spacing < 2 * overlap:
724 raise ValueError('Energy window too small. min/max: {}/{}'
725 .format(energy_limit_right, energy_limit_left) +
726 ' Try decreasing n_bins ({}) and/or overlap ({}).'
727 .format(n_bins, overlap))
728 if energy_limit_left is not None:
729 energy_limit_left -= overlap * energy_spacing
730 if energy_limit_right is not None:
731 energy_limit_right += overlap * energy_spacing
732 bounds.append((energy_limit_left, energy_limit_right))
734 return bounds