 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 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 """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 <vcsgc_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 Sadigh, B. and Erhart, P., Phys.

79 Rev. B **86**, 134204 (2012) [SadErh12]_.

81 Parameters

82 ----------

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

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

85 also defines the initial occupation vector

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

87 calculator to be used for calculating the potential changes

88 that enter the evaluation of the Metropolis criterion

89 temperature : float

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

91 phis : Dict[str, float]

92 average constraint parameters :math:\\phi_i; the key denotes the

93 species; for a N-component sublattice, there should be N - 1

94 different :math:\\phi_i (referred to as :math:\\bar{\\phi}

96 kappa : float

97 parameter that constrains the variance of the concentration

98 (referred to as :math:\\bar{\\kappa} in [SadErh12]_)

99 boltzmann_constant : float

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

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

106 random_seed : int

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

108 simulation

109 dc_filename : str

110 name of file the data container associated with the ensemble

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

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

113 updated/overwritten

114 data_container_write_period : float

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

116 written to file; writing periodically to file provides both

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

118 the data [default: 600 s]

119 ensemble_data_write_interval : int

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

121 includes for example the current value of the calculator

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

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

124 trajectory_write_interval : int

125 interval at which the current occupation vector of the atomic

126 configuration is written to the data container.

127 sublattice_probabilities : List[float]

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

129 The list should be as long as the number of sublattices and should

130 sum up to 1.

132 Example

133 -------

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

135 simulation in the variance-constrained semi-canonical ensemble. Here, the

136 parameters of the cluster expansion are set to emulate a simple Ising model

137 in order to obtain an example that can be run without modification. In

138 practice, one should of course use a proper cluster expansion::

140 >>> from ase.build import bulk

141 >>> from icet import ClusterExpansion, ClusterSpace

142 >>> from mchammer.calculators import ClusterExpansionCalculator

144 >>> # prepare cluster expansion

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

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

147 >>> # pairs are included)

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

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

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

152 >>> # set up and run MC simulation

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

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

155 >>> phi = 0.6

156 >>> mc = VCSGCEnsemble(structure=structure, calculator=calc,

157 ... temperature=600,

158 ... dc_filename='myrun_vcsgc.dc',

159 ... phis={'Au': phi},

160 ... kappa=200)

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

162 """

164 def __init__(self, structure: Atoms,

165 calculator: BaseCalculator,

166 temperature: float,

167 phis: Dict[str, float],

168 kappa: float,

169 boltzmann_constant: float = kB,

170 user_tag: str = None,

171 random_seed: int = None,

172 dc_filename: str = None,

173 data_container: str = None,

174 data_container_write_period: float = 600,

175 ensemble_data_write_interval: int = None,

176 trajectory_write_interval: int = None,

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

179 self._ensemble_parameters = dict(temperature=temperature,

180 kappa=kappa)

181 self._boltzmann_constant = boltzmann_constant

183 # Save ensemble parameters

184 for sym, phi in phis.items():

185 if isinstance(sym, str):

186 chemical_symbol = sym

187 else:

188 chemical_symbol = chemical_symbols[sym]

189 phi_sym = 'phi_{}'.format(chemical_symbol)

190 self._ensemble_parameters[phi_sym] = phi

192 super().__init__(

193 structure=structure,

194 calculator=calculator,

195 user_tag=user_tag,

196 random_seed=random_seed,

197 dc_filename=dc_filename,

198 data_container=data_container,

199 data_container_class=DataContainer,

200 data_container_write_period=data_container_write_period,

201 ensemble_data_write_interval=ensemble_data_write_interval,

202 trajectory_write_interval=trajectory_write_interval,

203 boltzmann_constant=boltzmann_constant

204 )

206 # Save phis (need self.configuration to check sublattices so

207 # we do it last)

208 self._phis = get_phis(phis)

210 # Check that each sublattice has N - 1 phis

211 for sl in self.sublattices.active_sublattices:

212 count_specified_elements = 0

213 for number in sl.atomic_numbers:

214 if number in self._phis.keys():

215 count_specified_elements += 1

216 if count_specified_elements != len(sl.atomic_numbers) - 1:

217 raise ValueError('phis must be set for N - 1 elements on a '

218 'sublattice with N species')

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

221 self._flip_sublattice_probabilities = self._get_flip_sublattice_probabilities()

222 else:

223 self._flip_sublattice_probabilities = sublattice_probabilities

225 def _do_trial_step(self):

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

227 sublattice_index = self.get_random_sublattice_index(

228 probability_distribution=self._flip_sublattice_probabilities)

229 return self.do_vcsgc_flip(

230 phis=self.phis, kappa=self.kappa, sublattice_index=sublattice_index)

232 @property

233 def temperature(self) -> float:

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

235 return self.ensemble_parameters['temperature']

237 @property

238 def phis(self) -> Dict[int, float]:

239 """

240 phis :math:\\phi_i for all species but one

241 (referred to as :math:\\bar{\\phi} in [SadErh12]_)

242 """

243 return self._phis

245 @property

246 def kappa(self) -> float:

247 """

248 kappa :math:\\bar{\\kappa} constrain parameter

249 (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 the phis that will be transformed to the format

276 used by the ensemble

277 """

278 if not isinstance(phis, dict):

279 raise TypeError('phis has the wrong type: {}'.format(type(phis)))

281 # Translate to atomic numbers if necessary

282 phis_ret = {}

283 for key, phi in phis.items():

284 if isinstance(key, str):

285 atomic_number = atomic_numbers[key]

286 phis_ret[atomic_number] = phi

287 elif isinstance(key, int): 287 ↛ 283line 287 didn't jump to line 283, because the condition on line 287 was never false

288 phis_ret[key] = phi

290 return phis_ret