Coverage for mchammer/ensembles/canonical_ensemble.py: 93%
24 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1"""Definition of the canonical ensemble class."""
3from ase import Atoms
4from ase.units import kB
5from typing import List
7from .. import DataContainer
8from ..calculators.base_calculator import BaseCalculator
9from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble
12class CanonicalEnsemble(ThermodynamicBaseEnsemble):
13 r"""Instances of this class allow one to simulate systems in the
14 canonical ensemble (:math:`N_iVT`), i.e. at constant temperature
15 (:math:`T`), number of atoms of each species (:math:`N_i`), and
16 volume (:math:`V`).
18 The probability for a particular state in the canonical ensemble is
19 proportional to the well-known Boltzmann factor,
21 .. math::
23 \rho_{\text{C}} \propto \exp [ - E / k_B T ].
25 Since the concentrations or equivalently the number of atoms of each
26 species is held fixed in the canonical ensemble, a trial step must
27 conserve the concentrations. This is accomplished by randomly picking two
28 unlike atoms and swapping their identities. The swap is accepted with
29 probability
31 .. math::
33 P = \min \{ 1, \, \exp [ - \Delta E / k_B T ] \},
35 where :math:`\Delta E` is the change in potential energy caused by the
36 swap.
38 The canonical ensemble provides an ideal framework for studying the
39 properties of a system at a specific concentration. Properties such as
40 potential energy or phenomena such as chemical ordering at a specific
41 temperature can conveniently be studied by simulating at that temperature.
42 The canonical ensemble is also a convenient tool for "optimizing" a
43 system, i.e., finding its lowest energy chemical ordering. In practice,
44 this is usually achieved by simulated annealing, i.e., the system is
45 equilibrated at a high temperature, after which the temperature is
46 continuously lowered until the acceptance probability is almost zero. In a
47 well-behaved system, the chemical ordering at that point corresponds to a
48 low-energy structure, possibly the global minimum at that particular
49 concentration.
51 Parameters
52 ----------
53 structure
54 Stomic configuration to be used in the Monte Carlo simulation;
55 also defines the initial occupation vector.
56 calculator
57 Calculator to be used for calculating the potential changes
58 that enter the evaluation of the Metropolis criterion.
59 temperature
60 Temperature :math:`T` in appropriate units, commonly Kelvin.
61 boltzmann_constant
62 Boltzmann constant :math:`k_B` in appropriate
63 units, i.e., units that are consistent
64 with the underlying cluster expansion
65 and the temperature units. Default: eV/K.
66 user_tag
67 Human-readable tag for ensemble. Default: ``None``.
68 random_seed
69 Seed for the random number generator used in the Monte Carlo simulation.
70 dc_filename
71 Name of file the data container associated with the ensemble
72 will be written to. If the file exists it will be read, the
73 data container will be appended, and the file will be
74 updated/overwritten.
75 data_container_write_period
76 Period in units of seconds at which the data container is
77 written to file. Writing periodically to file provides both
78 a way to examine the progress of the simulation and to back up
79 the data. Default: 600 s.
80 ensemble_data_write_interval
81 Interval at which data is written to the data container. This
82 includes for example the current value of the calculator
83 (i.e., usually the energy) as well as ensembles specific fields
84 such as temperature or the number of atoms of different species.
85 Default: Number of sites in the :attr:`structure`.
86 trajectory_write_interval
87 Interval at which the current occupation vector of the atomic
88 configuration is written to the data container.
89 Default: Number of sites in the :attr:`structure`.
90 sublattice_probabilities
91 Probability for picking a sublattice when doing a random swap.
92 This should be as long as the number of sublattices and should
93 sum up to 1.
96 Example
97 -------
98 The following snippet illustrate how to carry out a simple Monte Carlo
99 simulation in the canonical ensemble. Here, the parameters of the cluster
100 expansion are set to emulate a simple Ising model in order to obtain an
101 example that can be run without modification. In practice, one should of
102 course use a proper cluster expansion::
104 >>> from ase.build import bulk
105 >>> from icet import ClusterExpansion, ClusterSpace
106 >>> from mchammer.calculators import ClusterExpansionCalculator
108 >>> # prepare cluster expansion
109 >>> # the setup emulates a second nearest-neighbor (NN) Ising model
110 >>> # (zerolet and singlet ECIs are zero; only first and second neighbor
111 >>> # pairs are included)
112 >>> prim = bulk('Au')
113 >>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
114 >>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])
116 >>> # prepare initial configuration
117 >>> structure = prim.repeat(3)
118 >>> for k in range(5):
119 >>> structure[k].symbol = 'Ag'
121 >>> # set up and run MC simulation
122 >>> calc = ClusterExpansionCalculator(structure, ce)
123 >>> mc = CanonicalEnsemble(structure=structure, calculator=calc,
124 ... temperature=600,
125 ... dc_filename='myrun_canonical.dc')
126 >>> mc.run(100) # carry out 100 trial swaps
127 """
129 def __init__(self,
130 structure: Atoms,
131 calculator: BaseCalculator,
132 temperature: float,
133 user_tag: str = None,
134 boltzmann_constant: float = kB,
135 random_seed: int = None,
136 dc_filename: str = None,
137 data_container: str = None,
138 data_container_write_period: float = 600,
139 ensemble_data_write_interval: int = None,
140 trajectory_write_interval: int = None,
141 sublattice_probabilities: List[float] = None) -> None:
143 self._ensemble_parameters = dict(temperature=temperature)
145 # add species count to ensemble parameters
146 symbols = set([symbol for sub in calculator.sublattices
147 for symbol in sub.chemical_symbols])
148 for symbol in symbols:
149 key = 'n_atoms_{}'.format(symbol)
150 count = structure.get_chemical_symbols().count(symbol)
151 self._ensemble_parameters[key] = count
153 super().__init__(
154 structure=structure,
155 calculator=calculator,
156 user_tag=user_tag,
157 random_seed=random_seed,
158 data_container=data_container,
159 dc_filename=dc_filename,
160 data_container_class=DataContainer,
161 data_container_write_period=data_container_write_period,
162 ensemble_data_write_interval=ensemble_data_write_interval,
163 trajectory_write_interval=trajectory_write_interval,
164 boltzmann_constant=boltzmann_constant)
166 if sublattice_probabilities is None: 166 ↛ 169line 166 didn't jump to line 169, because the condition on line 166 was never false
167 self._swap_sublattice_probabilities = self._get_swap_sublattice_probabilities()
168 else:
169 self._swap_sublattice_probabilities = sublattice_probabilities
171 @property
172 def temperature(self) -> float:
173 """ Current temperature. """
174 return self._ensemble_parameters['temperature']
176 def _do_trial_step(self):
177 """ Carries out one Monte Carlo trial step. """
178 sublattice_index = self.get_random_sublattice_index(self._swap_sublattice_probabilities)
179 return self.do_canonical_swap(sublattice_index=sublattice_index)