Coverage for mchammer/ensembles/semi_grand_canonical_ensemble.py: 96%
38 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1"""
2Definition of the semi-grand canonical ensemble class.
3"""
5from ase import Atoms
6from ase.data import atomic_numbers, chemical_symbols
7from ase.units import kB
8from collections import OrderedDict
9from typing import Dict, Union, List
11from .. import DataContainer
12from ..calculators.base_calculator import BaseCalculator
13from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble
16class SemiGrandCanonicalEnsemble(ThermodynamicBaseEnsemble):
17 r"""Instances of this class allow one to simulate systems in the
18 semi-grand canonical (SGC) ensemble (:math:`N\Delta\mu_i VT`), i.e. at
19 constant temperature (:math:`T`), total number of sites (:math:`N=\sum_i
20 N_i`), relative chemical potentials (:math:`\Delta\mu_i=\mu_i - \mu_1`,
21 where :math:`i` denotes the species), and volume (:math:`V`).
23 The probability for a particular state in the SGC ensemble for a
24 :math:`m`-component system can be written
26 .. math::
28 \rho_{\text{SGC}} \propto \exp\Big[ - \big( E
29 + \sum_{i>1}^m \Delta\mu_i N_i \big) \big / k_B T \Big]
31 with the *relative* chemical potentials :math:`\Delta\mu_i = \mu_i -
32 \mu_1` and species counts :math:`N_i`. Unlike the :ref:`canonical ensemble
33 <canonical_ensemble>`, the number of the respective species (or
34 equivalently the concentrations) are allowed to vary in the SGC ensemble.
35 A trial step thus consists of randomly picking an atom and changing its
36 identity with probability
38 .. math::
40 P = \min \Big\{ 1, \, \exp \big[ - \big( \Delta E
41 + \sum_i \Delta \mu_i \Delta N_i \big) \big / k_B T \big]
42 \Big\},
44 where :math:`\Delta E` is the change in potential energy caused by the
45 swap.
47 There exists a simple relation between the differences in chemical
48 potential and the canonical free energy :math:`F`. In a binary system this
49 relationship reads
51 .. math:: \Delta \mu = - \frac{1}{N} \frac{\partial F}{\partial c} (
52 N, V, T, \langle c \rangle).
54 Here :math:`c` denotes concentration (:math:`c=N_i/N`) and :math:`\langle
55 c \rangle` the average concentration observed in the simulation. By
56 recording :math:`\langle c \rangle` while gradually changing
57 :math:`\Delta \mu`, one can thus in principle calculate the difference in
58 canonical free energy between the pure phases (:math:`c=0` or :math:`1`)
59 and any concentration by integrating :math:`\Delta \mu` over that
60 concentration range. In practice this requires that the recorded average
61 concentration :math:`\langle c \rangle` varies continuously with
62 :math:`\Delta \mu`. This is not the case for materials with multiphase
63 regions (such as miscibility gaps), because in such regions :math:`\Delta
64 \mu` maps to multiple concentrations. In a Monte Carlo simulation, this is
65 typically manifested by discontinuous jumps in concentration. Such jumps
66 mark the phase boundaries of a multiphase region and can thus be used to
67 construct the phase diagram. To recover the free energy, however, such
68 systems require sampling in other ensembles, such as the
69 :ref:`variance-constrained semi-grand canonical ensemble <vcsgc_ensemble>`.
71 Parameters
72 ----------
73 structure
74 Atomic configuration to be used in the Monte Carlo simulation;
75 also defines the initial occupation vector.
76 calculator
77 Calculator to be used for calculating the potential changes
78 that enter the evaluation of the Metropolis criterion.
79 temperature
80 Temperature :math:`T` in appropriate units, commonly Kelvin.
81 chemical_potentials
82 Chemical potential for each species :math:`\mu_i`. The key
83 denotes the species, the value specifies the chemical potential in
84 units that are consistent with the underlying cluster expansion.
85 boltzmann_constant
86 Boltzmann constant :math:`k_B` in appropriate units, i.e. units that are consistent
87 with the underlying cluster expansion and the temperature units. Default: eV/K.
88 user_tag
89 Human-readable tag for ensemble. Default: ``None``.
90 random_seed
91 Seed for the random number generator used in the Monte Carlo simulation.
92 dc_filename
93 Name of file the data container associated with the ensemble
94 will be written to. If the file exists it will be read, the
95 data container will be appended, and the file will be
96 updated/overwritten.
97 data_container_write_period
98 Period in units of seconds at which the data container is
99 written to file. Writing periodically to file provides both
100 a way to examine the progress of the simulation and to back up
101 the data. Default: 600 s.
102 ensemble_data_write_interval
103 Interval at which data is written to the data container. This
104 includes for example the current value of the calculator
105 (i.e., usually the energy) as well as ensembles specific fields
106 such as temperature or the number of atoms of different species.
107 Default: Number of sites in the :attr:`structure`.
108 trajectory_write_interval
109 Interval at which the current occupation vector of the atomic
110 configuration is written to the data container.
111 Default: Number of sites in the :attr:`structure`.
112 sublattice_probabilities
113 Probability for picking a sublattice when doing a random flip.
114 This should be as long as the number of sublattices and should
115 sum up to 1.
117 Example
118 -------
119 The following snippet illustrate how to carry out a simple Monte Carlo
120 simulation in the semi-canonical ensemble. Here, the parameters of the
121 cluster expansion are set to emulate a simple Ising model in order to
122 obtain an example that can be run without modification. In practice, one
123 should of course use a proper cluster expansion::
125 >>> from ase.build import bulk
126 >>> from icet import ClusterExpansion, ClusterSpace
127 >>> from mchammer.calculators import ClusterExpansionCalculator
129 >>> # prepare cluster expansion
130 >>> # the setup emulates a second nearest-neighbor (NN) Ising model
131 >>> # (zerolet and singlet ECIs are zero; only first and second neighbor
132 >>> # pairs are included)
133 >>> prim = bulk('Au')
134 >>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
135 >>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])
137 >>> # set up and run MC simulation (T=600 K, delta_mu=0.8 eV/atom)
138 >>> structure = prim.repeat(3)
139 >>> calc = ClusterExpansionCalculator(structure, ce)
140 >>> mc = SemiGrandCanonicalEnsemble(structure=structure, calculator=calc,
141 ... temperature=600,
142 ... dc_filename='myrun_sgc.dc',
143 ... chemical_potentials={'Ag': 0, 'Au': 0.8})
144 >>> mc.run(100) # carry out 100 trial swaps
146 """
148 def __init__(self,
149 structure: Atoms,
150 calculator: BaseCalculator,
151 temperature: float,
152 chemical_potentials: Dict[str, float],
153 boltzmann_constant: float = kB,
154 user_tag: str = None,
155 random_seed: int = None,
156 dc_filename: str = None,
157 data_container: str = None,
158 data_container_write_period: float = 600,
159 ensemble_data_write_interval: int = None,
160 trajectory_write_interval: int = None,
161 sublattice_probabilities: List[float] = None) -> None:
163 self._ensemble_parameters = dict(temperature=temperature)
165 # add chemical potentials to ensemble parameters
166 # TODO: add check that chemical symbols in chemical potentials are allowed
167 self._chemical_potentials = get_chemical_potentials(chemical_potentials)
168 for atnum, chempot in self.chemical_potentials.items():
169 mu_sym = 'mu_{}'.format(chemical_symbols[atnum])
170 self._ensemble_parameters[mu_sym] = chempot
172 self._boltzmann_constant = boltzmann_constant
174 super().__init__(
175 structure=structure, calculator=calculator, user_tag=user_tag,
176 random_seed=random_seed,
177 dc_filename=dc_filename,
178 data_container=data_container,
179 data_container_class=DataContainer,
180 data_container_write_period=data_container_write_period,
181 ensemble_data_write_interval=ensemble_data_write_interval,
182 trajectory_write_interval=trajectory_write_interval,
183 boltzmann_constant=boltzmann_constant)
185 if sublattice_probabilities is None: 185 ↛ 188line 185 didn't jump to line 188, because the condition on line 185 was never false
186 self._flip_sublattice_probabilities = self._get_flip_sublattice_probabilities()
187 else:
188 self._flip_sublattice_probabilities = sublattice_probabilities
190 @property
191 def temperature(self) -> float:
192 """ Temperature :math:`T` (see parameters section above). """
193 return self.ensemble_parameters['temperature']
195 def _do_trial_step(self):
196 """ Carries out one Monte Carlo trial step. """
197 sublattice_index = self.get_random_sublattice_index(
198 probability_distribution=self._flip_sublattice_probabilities)
199 return self.do_sgc_flip(self.chemical_potentials, sublattice_index=sublattice_index)
201 @property
202 def chemical_potentials(self) -> Dict[int, float]:
203 r"""
204 Chemical potentials :math:`\mu_i` (see parameters section above).
205 """
206 return self._chemical_potentials
208 def _get_ensemble_data(self) -> Dict:
209 """Returns the data associated with the ensemble. For the SGC
210 ensemble this specifically includes the species counts.
211 """
212 # generic data
213 data = super()._get_ensemble_data()
215 # species counts
216 data.update(self._get_species_counts())
218 return data
221def get_chemical_potentials(chemical_potentials: Union[Dict[str, float], Dict[int, float]]) \
222 -> Dict[int, float]:
223 """ Gets values of chemical potentials. """
224 if not isinstance(chemical_potentials, dict):
225 raise TypeError('chemical_potentials has the wrong type: {}'
226 .format(type(chemical_potentials)))
228 cps = OrderedDict([(key, val) if isinstance(key, int)
229 else (atomic_numbers[key], val)
230 for key, val in chemical_potentials.items()])
232 return cps