Coverage for mchammer/ensembles/vcsgc_ensemble.py: 96%
57 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1"""
2Definition of the variance-constrained semi-grand canonical ensemble class.
3"""
5from ase import Atoms
6from ase.data import atomic_numbers, chemical_symbols
7from ase.units import kB
8from typing import Dict, Union, List
10from .. import DataContainer
11from ..calculators.base_calculator import BaseCalculator
12from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble
15class VCSGCEnsemble(ThermodynamicBaseEnsemble):
16 r"""Instances of this class allow one to simulate systems in the
17 variance-constrained semi-grand canonical (VCSGC) ensemble
18 (:math:`N\phi\kappa VT`), i.e. at constant temperature (:math:`T`), total
19 number of sites (:math:`N=\sum_i N_i`), and two additional dimensionless
20 parameters :math:`\phi` and :math:`\kappa`, which constrain average and
21 variance of the concentration, respectively.
23 The below examples treat the binary case, but the generalization of
24 to ternaries and higher-order systems is straight-forward.
25 The probability for a particular state in the VCSGC ensemble for a
26 :math:`2`-component system can be written
28 .. math::
30 \rho_{\text{VCSGC}} \propto \exp\Big[ - E / k_B T
31 + \kappa N ( c_1 + \phi_1 / 2 )^2 \Big],
33 where :math:`c_1` represents the concentration of species 1, i.e.,
34 :math:`c_1=N_1/N`. (Please note that the quantities :math:`\kappa` and
35 :math:`\phi` correspond, respectively, to :math:`\bar{\kappa}` and
36 :math:`\bar{\phi}` in [SadErh12]_.) The :math:`\phi` may refer to any
37 of the two species. If :math:`\phi` is specified for species A, an
38 equivalent simulation can be carried out by specifying :math:`\phi_B` as
39 :math:`-2-\phi_A`. In general, simulations of :math:`N`-component
40 systems requires the specification of :math:`\phi` for :math:`N-1`
41 elements.
43 Just like the :ref:`semi-grand canonical ensemble <canonical_ensemble>`,
44 the VCSGC ensemble allows concentrations to change. A trial step consists
45 of changing the identity of a randomly chosen atom and accepting the change
46 with probability
48 .. math::
50 P = \min \{ 1, \, \exp [ - \Delta E / k_B T
51 + \kappa N \Delta c_1 (\phi_1 + \Delta c_1 + 2 c_1 ) ] \}.
53 Note that for a sufficiently large value of :math:`\kappa`, say 200, the
54 probability density :math:`\rho_{\text{VCSGC}}` is sharply peaked around
55 :math:`c_1=-\phi_1 / 2`. In practice, this means that we can gradually
56 change :math:`\phi_1` from (using some margins) :math:`-2.1` to
57 :math:`0.1` and take the system continuously from :math:`c_1 = 0` to
58 :math:`1`. The parameter :math:`\kappa` constrains the fluctuations (or
59 the variance) of the concentration at each value of :math:`\phi_1`, with
60 higher values of :math:`\kappa` meaning less fluctuations. Unlike the
61 :ref:`semi-grand canonical ensemble <sgc_ensemble>`, one value of
62 :math:`\phi_1` maps to one and only one concentration also in multiphase
63 regions. Since the derivative of the canonical free energy can be expressed
64 in terms of parameters and observables of the VCSGC ensemble,
66 .. math::
68 k_B T \kappa ( \phi_1 + 2 \langle c_1 \rangle ) = - \frac{1}{N}
69 \frac{\partial F}{\partial c_1} (N, V, T, \langle c_1 \rangle ),
71 this ensemble allows for thermodynamic integration across multiphase
72 regions. This means that we can construct phase diagrams by directly
73 comparing the free energies of the different phases. This often makes the
74 VCSGC ensemble more convenient than the :ref:`semi-grand canonical ensemble
75 <sgc_ensemble>` when simulating materials with multiphase regions, such as
76 alloys with miscibility gaps.
78 When using the VCSGC ensemble, please cite
79 Phys. Rev. B **86**, 134204 (2012) [SadErh12]_.
81 Parameters
82 ----------
83 structure
84 Atomic configuration to be used in the Monte Carlo simulation;
85 also defines the initial occupation vector.
86 calculator
87 Calculator to be used for calculating the potential changes
88 that enter the evaluation of the Metropolis criterion.
89 temperature
90 Temperature :math:`T` in appropriate units, commonly Kelvin.
91 phis
92 Average constraint parameters :math:`\phi_i`. The key denotes the
93 species. For a :math:`N`-component sublattice, there should be :math:`N - 1`
94 different :math:`\phi_i` (referred to as :math:`\bar{\phi}`
95 in [SadErh12]_).
96 kappa
97 parameter that constrains the variance of the concentration
98 (referred to as :math:`\bar{\kappa}` in [SadErh12]_)
99 boltzmann_constant
100 Boltzmann constant :math:`k_B` in appropriate
101 units, i.e., units that are consistent
102 with the underlying cluster expansion
103 and the temperature units. Default: eV/K.
104 user_tag
105 Human-readable tag for ensemble. Default: ``None``.
106 random_seed
107 Seed for the random number generator used in the Monte Carlo simulation.
108 dc_filename
109 Name of file the data container associated with the ensemble
110 will be written to. If the file exists it will be read, the
111 data container will be appended, and the file will be
112 updated/overwritten.
113 data_container_write_period
114 Period in units of seconds at which the data container is
115 written to file. Writing periodically to file provides both
116 a way to examine the progress of the simulation and to back up
117 the data. Default: 600 s.
118 ensemble_data_write_interval
119 Interval at which data is written to the data container. This
120 includes for example the current value of the calculator
121 (i.e., usually the energy) as well as ensembles specific fields
122 such as temperature or the number of atoms of different species.
123 Default: Number of sites in the :attr:`structure`.
124 trajectory_write_interval
125 Interval at which the current occupation vector of the atomic
126 configuration is written to the data container.
127 Default: Number of sites in the :attr:`structure`.
128 sublattice_probabilities
129 Probability for picking a sublattice when doing a random flip.
130 The list should be as long as the number of sublattices and should
131 sum up to 1.
133 Example
134 -------
135 The following snippet illustrate how to carry out a simple Monte Carlo
136 simulation in the variance-constrained semi-canonical ensemble. Here, the
137 parameters of the cluster expansion are set to emulate a simple Ising model
138 in order to obtain an example that can be run without modification. In
139 practice, one should of course use a proper cluster expansion::
141 >>> from ase.build import bulk
142 >>> from icet import ClusterExpansion, ClusterSpace
143 >>> from mchammer.calculators import ClusterExpansionCalculator
145 >>> # prepare cluster expansion
146 >>> # the setup emulates a second nearest-neighbor (NN) Ising model
147 >>> # (zerolet and singlet ECIs are zero; only first and second neighbor
148 >>> # pairs are included)
149 >>> prim = bulk('Au')
150 >>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
151 >>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])
153 >>> # set up and run MC simulation
154 >>> structure = prim.repeat(3)
155 >>> calc = ClusterExpansionCalculator(structure, ce)
156 >>> phi = 0.6
157 >>> mc = VCSGCEnsemble(structure=structure, calculator=calc,
158 ... temperature=600,
159 ... dc_filename='myrun_vcsgc.dc',
160 ... phis={'Au': phi},
161 ... kappa=200)
162 >>> mc.run(100) # carry out 100 trial swaps
163 """
165 def __init__(self, structure: Atoms,
166 calculator: BaseCalculator,
167 temperature: float,
168 phis: Dict[str, float],
169 kappa: float,
170 boltzmann_constant: float = kB,
171 user_tag: str = None,
172 random_seed: int = None,
173 dc_filename: str = None,
174 data_container: str = None,
175 data_container_write_period: float = 600,
176 ensemble_data_write_interval: int = None,
177 trajectory_write_interval: int = None,
178 sublattice_probabilities: List[float] = None) -> None:
180 self._ensemble_parameters = dict(temperature=temperature,
181 kappa=kappa)
182 self._boltzmann_constant = boltzmann_constant
184 # Save ensemble parameters
185 for sym, phi in phis.items():
186 if isinstance(sym, str):
187 chemical_symbol = sym
188 else:
189 chemical_symbol = chemical_symbols[sym]
190 phi_sym = 'phi_{}'.format(chemical_symbol)
191 self._ensemble_parameters[phi_sym] = phi
193 super().__init__(
194 structure=structure,
195 calculator=calculator,
196 user_tag=user_tag,
197 random_seed=random_seed,
198 dc_filename=dc_filename,
199 data_container=data_container,
200 data_container_class=DataContainer,
201 data_container_write_period=data_container_write_period,
202 ensemble_data_write_interval=ensemble_data_write_interval,
203 trajectory_write_interval=trajectory_write_interval,
204 boltzmann_constant=boltzmann_constant
205 )
207 # Save phis (need self.configuration to check sublattices so
208 # we do it last)
209 self._phis = get_phis(phis)
211 # Check that each sublattice has N - 1 phis
212 for sl in self.sublattices.active_sublattices:
213 count_specified_elements = 0
214 for number in sl.atomic_numbers:
215 if number in self._phis.keys():
216 count_specified_elements += 1
217 if count_specified_elements != len(sl.atomic_numbers) - 1:
218 raise ValueError('phis must be set for N - 1 elements on a '
219 'sublattice with N species')
221 if sublattice_probabilities is None: 221 ↛ 224line 221 didn't jump to line 224, because the condition on line 221 was never false
222 self._flip_sublattice_probabilities = self._get_flip_sublattice_probabilities()
223 else:
224 self._flip_sublattice_probabilities = sublattice_probabilities
226 def _do_trial_step(self):
227 """ Carries out one Monte Carlo trial step. """
228 sublattice_index = self.get_random_sublattice_index(
229 probability_distribution=self._flip_sublattice_probabilities)
230 return self.do_vcsgc_flip(
231 phis=self.phis, kappa=self.kappa, sublattice_index=sublattice_index)
233 @property
234 def temperature(self) -> float:
235 """ Temperature :math:`T` (see parameters section above). """
236 return self.ensemble_parameters['temperature']
238 @property
239 def phis(self) -> Dict[int, float]:
240 r"""
241 Average constraint parameters :math:`\phi_i` for all species but one
242 (referred to as :math:`\bar{\phi}` in [SadErh12]_).
243 """
244 return self._phis
246 @property
247 def kappa(self) -> float:
248 r"""
249 Variance constraint parameter :math:`\bar{\kappa}` (see parameters section above).
250 """
251 return self.ensemble_parameters['kappa']
253 def _get_ensemble_data(self) -> Dict:
254 """
255 Returns a dict with the default data of the ensemble. This includes
256 atom counts and free energy derivative.
257 """
258 data = super()._get_ensemble_data()
260 # free energy derivative
261 data.update(self._get_vcsgc_free_energy_derivatives(self.phis, self.kappa))
263 # species counts
264 data.update(self._get_species_counts())
266 return data
269def get_phis(phis: Union[Dict[int, float], Dict[str, float]]) -> Dict[int, float]:
270 """Get phis as used in the VCSGC ensemble.
272 Parameters
273 ----------
274 phis
275 phis that will be transformed to the format used by the ensemble.
276 """
277 if not isinstance(phis, dict):
278 raise TypeError(f'phis has the wrong type: {type(phis)}')
280 # Translate to atomic numbers if necessary
281 phis_ret = {}
282 for key, phi in phis.items():
283 if isinstance(key, str):
284 atomic_number = atomic_numbers[key]
285 phis_ret[atomic_number] = phi
286 elif isinstance(key, int): 286 ↛ 282line 286 didn't jump to line 282, because the condition on line 286 was never false
287 phis_ret[key] = phi
289 return phis_ret