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