Hide keyboard shortcuts

Hot-keys on this page

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

4 

5from ase import Atoms 

6from ase.data import atomic_numbers, chemical_symbols 

7from ase.units import kB 

8from typing import Dict, Union, List 

9 

10from .. import DataContainer 

11from ..calculators.base_calculator import BaseCalculator 

12from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble 

13 

14 

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. 

22 

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 

27 

28 .. math:: 

29 

30 \\rho_{\\text{VCSGC}} \\propto \\exp\\Big[ - E / k_B T 

31 + \\kappa N ( c_1 + \\phi_1 / 2 )^2 \\Big], 

32 

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. 

42 

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 

47 

48 .. math:: 

49 

50 P = \\min \\{ 1, \\, \\exp [ - \\Delta E / k_B T 

51 + \\kappa N \\Delta c_1 (\\phi_1 + \\Delta c_1 + 2 c_1 ) ] \\}. 

52 

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, 

65 

66 .. math:: 

67 

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 ), 

70 

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. 

77 

78 When using the VCSGC ensemble, please cite Sadigh, B. and Erhart, P., Phys. 

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

80 

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}` 

95 in [SadErh12]_) 

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. 

131 

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

139 

140 >>> from ase.build import bulk 

141 >>> from icet import ClusterExpansion, ClusterSpace 

142 >>> from mchammer.calculators import ClusterExpansionCalculator 

143 

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]) 

151 

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

163 

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: 

178 

179 self._ensemble_parameters = dict(temperature=temperature, 

180 kappa=kappa) 

181 self._boltzmann_constant = boltzmann_constant 

182 

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 

191 

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 ) 

205 

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

207 # we do it last) 

208 self._phis = get_phis(phis) 

209 

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') 

219 

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 

224 

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) 

231 

232 @property 

233 def temperature(self) -> float: 

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

235 return self.ensemble_parameters['temperature'] 

236 

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 

244 

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'] 

252 

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() 

259 

260 # free energy derivative 

261 data.update(self._get_vcsgc_free_energy_derivatives(self.phis, self.kappa)) 

262 

263 # species counts 

264 data.update(self._get_species_counts()) 

265 

266 return data 

267 

268 

269def get_phis(phis: Union[Dict[int, float], Dict[str, float]]) -> Dict[int, float]: 

270 """Get phis as used in the vcsgc ensemble. 

271 

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

280 

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 

289 

290 return phis_ret