 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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 """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

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 average recorded

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 <sgc_ensemble>.

71 Parameters

72 ----------

73 structure : :class:Atoms <ase.Atoms>

74 atomic configuration to be used in the Monte Carlo simulation;

75 also defines the initial occupation vector

76 calculator : :class:BaseCalculator <mchammer.calculators.ClusterExpansionCalculator>

77 calculator to be used for calculating the potential changes

78 that enter the evaluation of the Metropolis criterion

79 temperature : float

80 temperature :math:T in appropriate units [commonly Kelvin]

81 chemical_potentials : Dict[str, float]

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 : float

86 Boltzmann constant :math:k_B in appropriate

87 units, i.e. units that are consistent

88 with the underlying cluster expansion

89 and the temperature units [default: eV/K]

90 user_tag : str

91 human-readable tag for ensemble [default: None]

92 random_seed : int

93 seed for the random number generator used in the Monte Carlo

94 simulation

95 dc_filename : str

96 name of file the data container associated with the ensemble

97 will be written to; if the file exists it will be read, the

98 data container will be appended, and the file will be

99 updated/overwritten

100 data_container_write_period : float

101 period in units of seconds at which the data container is

102 written to file; writing periodically to file provides both

103 a way to examine the progress of the simulation and to back up

104 the data [default: 600 s]

105 ensemble_data_write_interval : int

106 interval at which data is written to the data container; this

107 includes for example the current value of the calculator

108 (i.e. usually the energy) as well as ensembles specific fields

109 such as temperature or the number of atoms of different species

110 trajectory_write_interval : int

111 interval at which the current occupation vector of the atomic

112 configuration is written to the data container.

113 sublattice_probabilities : List[float]

114 probability for picking a sublattice when doing a random flip.

115 This should be as long as the number of sublattices and should

116 sum up to 1.

118 Example

119 -------

120 The following snippet illustrate how to carry out a simple Monte Carlo

121 simulation in the semi-canonical ensemble. Here, the parameters of the

122 cluster expansion are set to emulate a simple Ising model in order to

123 obtain an example that can be run without modification. In practice, one

124 should of course use a proper cluster expansion::

126 >>> from ase.build import bulk

127 >>> from icet import ClusterExpansion, ClusterSpace

128 >>> from mchammer.calculators import ClusterExpansionCalculator

130 >>> # prepare cluster expansion

131 >>> # the setup emulates a second nearest-neighbor (NN) Ising model

132 >>> # (zerolet and singlet ECIs are zero; only first and second neighbor

133 >>> # pairs are included)

134 >>> prim = bulk('Au')

135 >>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])

136 >>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])

138 >>> # set up and run MC simulation (T=600 K, delta_mu=0.8 eV/atom)

139 >>> structure = prim.repeat(3)

140 >>> calc = ClusterExpansionCalculator(structure, ce)

141 >>> mc = SemiGrandCanonicalEnsemble(structure=structure, calculator=calc,

142 ... temperature=600,

143 ... dc_filename='myrun_sgc.dc',

144 ... chemical_potentials={'Ag': 0, 'Au': 0.8})

145 >>> mc.run(100) # carry out 100 trial swaps

147 """

149 def __init__(self,

150 structure: Atoms,

151 calculator: BaseCalculator,

152 temperature: float,

153 chemical_potentials: Dict[str, float],

154 boltzmann_constant: float = kB,

155 user_tag: str = None,

156 random_seed: int = None,

157 dc_filename: str = None,

158 data_container: str = None,

159 data_container_write_period: float = 600,

160 ensemble_data_write_interval: int = None,

161 trajectory_write_interval: int = None,

162 sublattice_probabilities: List[float] = None) -> None:

164 self._ensemble_parameters = dict(temperature=temperature)

166 # add chemical potentials to ensemble parameters

167 # TODO: add check that chemical symbols in chemical potentials are allowed

168 self._chemical_potentials = get_chemical_potentials(chemical_potentials)

169 for atnum, chempot in self.chemical_potentials.items():

170 mu_sym = 'mu_{}'.format(chemical_symbols[atnum])

171 self._ensemble_parameters[mu_sym] = chempot

173 self._boltzmann_constant = boltzmann_constant

175 super().__init__(

176 structure=structure, calculator=calculator, user_tag=user_tag,

177 random_seed=random_seed,

178 dc_filename=dc_filename,

179 data_container=data_container,

180 data_container_class=DataContainer,

181 data_container_write_period=data_container_write_period,

182 ensemble_data_write_interval=ensemble_data_write_interval,

183 trajectory_write_interval=trajectory_write_interval,

184 boltzmann_constant=boltzmann_constant)

186 if sublattice_probabilities is None: 186 ↛ 189line 186 didn't jump to line 189, because the condition on line 186 was never false

187 self._flip_sublattice_probabilities = self._get_flip_sublattice_probabilities()

188 else:

189 self._flip_sublattice_probabilities = sublattice_probabilities

191 @property

192 def temperature(self) -> float:

193 """ temperature :math:T (see parameters section above) """

194 return self.ensemble_parameters['temperature']

196 def _do_trial_step(self):

197 """ Carries out one Monte Carlo trial step. """

198 sublattice_index = self.get_random_sublattice_index(

199 probability_distribution=self._flip_sublattice_probabilities)

200 return self.do_sgc_flip(self.chemical_potentials, sublattice_index=sublattice_index)

202 @property

203 def chemical_potentials(self) -> Dict[int, float]:

204 """

205 chemical potentials :math:\\mu_i (see parameters section above)

206 """

207 return self._chemical_potentials

209 def _get_ensemble_data(self) -> Dict:

210 """Returns the data associated with the ensemble. For the SGC

211 ensemble this specifically includes the species counts.

212 """

213 # generic data

214 data = super()._get_ensemble_data()

216 # species counts

217 data.update(self._get_species_counts())

219 return data

222def get_chemical_potentials(chemical_potentials: Union[Dict[str, float], Dict[int, float]]) \

223 -> Dict[int, float]:

224 """ Gets values of chemical potentials."""

225 if not isinstance(chemical_potentials, dict):

226 raise TypeError('chemical_potentials has the wrong type: {}'

227 .format(type(chemical_potentials)))

229 cps = OrderedDict([(key, val) if isinstance(key, int)

230 else (atomic_numbers[key], val)

231 for key, val in chemical_potentials.items()])

233 return cps