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

1from abc import abstractproperty 

2from typing import Dict, List, Type, Any 

3 

4import numpy as np 

5 

6from ase import Atoms 

7from ase.units import kB 

8from ase.data import chemical_symbols 

9 

10from ..calculators.base_calculator import BaseCalculator 

11from ..data_containers.base_data_container import BaseDataContainer 

12from .base_ensemble import BaseEnsemble 

13 

14 

15class ThermodynamicBaseEnsemble(BaseEnsemble): 

16 """ 

17 Thermodynamic base ensemble class. 

18 

19 Parameters 

20 ---------- 

21 structure : :class:`Atoms <ase.Atoms>` 

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

23 also defines the initial occupation vector 

24 calculator : :class:`BaseCalculator <mchammer.calculators.ClusterExpansionCalculator>` 

25 calculator to be used for calculating the potential changes 

26 that enter the evaluation of the Metropolis criterion 

27 boltzmann_constant : float 

28 Boltzmann constant :math:`k_B` in appropriate 

29 units, i.e. units that are consistent 

30 with the underlying cluster expansion 

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

32 user_tag : str 

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

34 random_seed : int 

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

36 simulation 

37 dc_filename : str 

38 name of file the data container associated with the ensemble 

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

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

41 updated/overwritten 

42 data_container_write_period : float 

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

44 written to file; writing periodically to file provides both 

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

46 the data [default: 600 s] 

47 ensemble_data_write_interval : int 

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

49 includes for example the current value of the calculator 

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

51 such as temperature or the number of structure of different species 

52 trajectory_write_interval : int 

53 interval at which the current occupation vector of the atomic 

54 configuration is written to the data container. 

55 """ 

56 

57 def __init__(self, 

58 structure: Atoms, 

59 calculator: BaseCalculator, 

60 user_tag: str = None, 

61 boltzmann_constant: float = kB, 

62 random_seed: int = None, 

63 dc_filename: str = None, 

64 data_container: str = None, 

65 data_container_class: Type[BaseDataContainer] = None, 

66 data_container_write_period: float = 600, 

67 ensemble_data_write_interval: int = None, 

68 trajectory_write_interval: int = None) -> None: 

69 

70 self._boltzmann_constant = boltzmann_constant 

71 

72 super().__init__( 

73 structure=structure, 

74 calculator=calculator, 

75 user_tag=user_tag, 

76 random_seed=random_seed, 

77 dc_filename=dc_filename, 

78 data_container=data_container, 

79 data_container_class=data_container_class, 

80 data_container_write_period=data_container_write_period, 

81 ensemble_data_write_interval=ensemble_data_write_interval, 

82 trajectory_write_interval=trajectory_write_interval) 

83 

84 @abstractproperty 

85 @property 

86 def temperature(self) -> float: 

87 pass 

88 

89 @property 

90 def boltzmann_constant(self) -> float: 

91 """ Boltzmann constant :math:`k_B` (see parameters section above) """ 

92 return self._boltzmann_constant 

93 

94 def _acceptance_condition(self, potential_diff: float) -> bool: 

95 """ 

96 Evaluates Metropolis acceptance criterion. 

97 

98 Parameters 

99 ---------- 

100 potential_diff 

101 change in the thermodynamic potential associated 

102 with the trial step 

103 """ 

104 if potential_diff <= 0: 

105 return True 

106 elif self.temperature <= 1e-16: 

107 return False 

108 else: 

109 p = np.exp(-potential_diff / (self.boltzmann_constant * self.temperature)) 

110 return p > self._next_random_number() 

111 

112 def do_canonical_swap(self, sublattice_index: int, allowed_species: List[int] = None) -> int: 

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

114 

115 Parameters 

116 --------- 

117 sublattice_index 

118 the sublattice the swap will act on 

119 allowed_species 

120 list of atomic numbers for allowed species 

121 

122 Returns 

123 ------- 

124 Returns 1 or 0 depending on if trial move was accepted or rejected 

125 """ 

126 sites, species = self.configuration.get_swapped_state(sublattice_index, allowed_species) 

127 potential_diff = self._get_property_change(sites, species) 

128 

129 if self._acceptance_condition(potential_diff): 

130 self.update_occupations(sites, species) 

131 return 1 

132 return 0 

133 

134 def do_sgc_flip(self, chemical_potentials: Dict[int, float], sublattice_index: int, 

135 allowed_species: List[int] = None) -> int: 

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

137 

138 Parameters 

139 --------- 

140 chemical_potentials 

141 chemical potentials used to calculate the potential 

142 difference 

143 sublattice_index 

144 the sublattice the flip will act on 

145 allowed_species 

146 list of atomic numbers for allowed species 

147 

148 Returns 

149 ------- 

150 Returns 1 or 0 depending on if trial move was accepted or rejected 

151 """ 

152 index, species = self.configuration.get_flip_state(sublattice_index, allowed_species) 

153 potential_diff = self._get_property_change([index], [species]) 

154 

155 # change in chemical potential 

156 old_species = self.configuration.occupations[index] 

157 chemical_potential_diff = chemical_potentials[old_species] - chemical_potentials[species] 

158 potential_diff += chemical_potential_diff 

159 

160 if self._acceptance_condition(potential_diff): 

161 self.update_occupations([index], [species]) 

162 return 1 

163 return 0 

164 

165 def do_vcsgc_flip(self, phis: Dict[int, float], kappa: float, sublattice_index: int, 

166 allowed_species: List[int] = None) -> int: 

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

168 

169 Parameters 

170 ---------- 

171 phis 

172 average constraint parameters 

173 kappa 

174 parameter that constrains the variance of the concentration 

175 sublattice_index 

176 the sublattice the flip will act on 

177 allowed_species 

178 list of atomic numbers for allowed species 

179 

180 Returns 

181 ------- 

182 Returns 1 or 0 depending on if trial move was accepted or rejected 

183 """ 

184 index, new_species = self.configuration.get_flip_state( 

185 sublattice_index, allowed_species) 

186 old_species = self.configuration.occupations[index] 

187 

188 # Calculate difference in VCSGC thermodynamic potential. 

189 # Note that this assumes that only one atom was flipped. 

190 sl_occupations = self.configuration.get_occupations_on_sublattice(sublattice_index) 

191 N = len(sl_occupations) 

192 potential_diff = 1.0 # dN 

193 for species in phis: 

194 if species == old_species: 

195 factor = -1 

196 elif species == new_species: 196 ↛ 199line 196 didn't jump to line 199, because the condition on line 196 was never false

197 factor = 1 

198 else: 

199 continue 

200 potential_diff += factor * (N * phis[species] + 2 * sl_occupations.count(species)) 

201 potential_diff *= kappa * self.boltzmann_constant * self.temperature / N 

202 potential_diff += self._get_property_change([index], [new_species]) 

203 

204 if self._acceptance_condition(potential_diff): 

205 self.update_occupations([index], [new_species]) 

206 return 1 

207 return 0 

208 

209 def _get_swap_sublattice_probabilities(self) -> List[float]: 

210 """ Returns sublattice probabilities suitable for swaps.""" 

211 sublattice_probabilities = [] # type: List[Any] 

212 for i, sl in enumerate(self.sublattices): 

213 if self.configuration.is_swap_possible(i): 

214 sublattice_probabilities.append(len(sl.indices)) 

215 else: 

216 sublattice_probabilities.append(0) 

217 norm = sum(sublattice_probabilities) 

218 if norm == 0: 

219 raise ValueError('No canonical swaps are possible on any of the active sublattices.') 

220 sublattice_probabilities = [p / norm for p in sublattice_probabilities] 

221 return sublattice_probabilities 

222 

223 def _get_flip_sublattice_probabilities(self) -> List[float]: 

224 """Returns the default sublattice probability which is based on 

225 the sizes of a sublattice. 

226 """ 

227 probability_distribution = [] # type: List[Any] 

228 for i, sl in enumerate(self.sublattices): 

229 if len(sl.chemical_symbols) > 1: 

230 probability_distribution.append(len(sl.indices)) 

231 else: 

232 probability_distribution.append(0) 

233 

234 norm = sum(probability_distribution) 

235 probability_distribution = [p / norm for p in probability_distribution] 

236 return probability_distribution 

237 

238 def _get_vcsgc_free_energy_derivatives(self, phis: Dict[int, float], kappa: float, 

239 sublattice_index: int = None) -> Dict: 

240 """ 

241 Returns a dict with the free energy derivatives. 

242 

243 Parameters 

244 ---------- 

245 phis 

246 average constraint parameters 

247 kappa 

248 parameter that constrains the variance of the concentration 

249 sublattice_index 

250 sublattice index 

251 """ 

252 data = {} 

253 

254 for atnum in phis: 

255 for i, sublattice in enumerate(self.sublattices): 

256 if sublattice_index is not None and i != sublattice_index: 

257 continue 

258 if len(sublattice.chemical_symbols) > 0 and atnum in sublattice.atomic_numbers: 

259 N = len(sublattice.indices) 

260 sl_occupations = self.configuration.get_occupations_on_sublattice(i) 

261 concentration = sl_occupations.count(atnum) / N 

262 data['free_energy_derivative_{}'.format(chemical_symbols[atnum])] \ 

263 = kappa * self.boltzmann_constant * self.temperature * \ 

264 (- 2 * concentration - phis[atnum]) 

265 

266 return data 

267 

268 def _get_species_counts(self) -> Dict: 

269 """ 

270 Returns a dict with the species counts. 

271 """ 

272 data = {} 

273 structure = self.configuration.structure 

274 unique, counts = np.unique(structure.numbers, return_counts=True) 

275 for sl in self.sublattices: 

276 for symbol in sl.chemical_symbols: 

277 data['{}_count'.format(symbol)] = 0 

278 for atnum, count in zip(unique, counts): 

279 data['{}_count'.format(chemical_symbols[atnum])] = count 

280 

281 return data