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 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 collections import OrderedDict 

9from typing import Dict, Union, List 

10 

11from .. import DataContainer 

12from ..calculators.base_calculator import BaseCalculator 

13from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble 

14 

15 

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

22 

23 The probability for a particular state in the SGC ensemble for a 

24 :math:`m`-component system can be written 

25 

26 .. math:: 

27 

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] 

30 

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 

37 

38 .. math:: 

39 

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\\}, 

43 

44 where :math:`\\Delta E` is the change in potential energy caused by the 

45 swap. 

46 

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 

50 

51 .. math:: \\Delta \\mu = - \\frac{1}{N} \\frac{\\partial F}{\\partial c} ( 

52 N, V, T, \\langle c \\rangle). 

53 

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

70 

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. 

117 

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

125 

126 >>> from ase.build import bulk 

127 >>> from icet import ClusterExpansion, ClusterSpace 

128 >>> from mchammer.calculators import ClusterExpansionCalculator 

129 

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

137 

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 

146 

147 """ 

148 

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: 

163 

164 self._ensemble_parameters = dict(temperature=temperature) 

165 

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 

172 

173 self._boltzmann_constant = boltzmann_constant 

174 

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) 

185 

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 

190 

191 @property 

192 def temperature(self) -> float: 

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

194 return self.ensemble_parameters['temperature'] 

195 

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) 

201 

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 

208 

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

215 

216 # species counts 

217 data.update(self._get_species_counts()) 

218 

219 return data 

220 

221 

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

228 

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

230 else (atomic_numbers[key], val) 

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

232 

233 return cps