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 typing import Dict, List, Any 

2 

3import numpy as np 

4from math import isclose 

5 

6from ase import Atoms 

7from ase.units import kB 

8from ase.data import atomic_numbers, chemical_symbols 

9 

10from .. import DataContainer 

11from ..calculators.base_calculator import BaseCalculator 

12from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble 

13from .vcsgc_ensemble import get_phis 

14from .semi_grand_canonical_ensemble import get_chemical_potentials 

15 

16 

17class HybridEnsemble(ThermodynamicBaseEnsemble): 

18 """ 

19 

20 Instances of this class allows one to combine multiple ensembles. 

21 In particular, a dictionary should be provided for each ensemble, 

22 which must include the type (`ensemble`) as well as the index of 

23 the sublattice (`sublattice_index`). In addition, it is possible 

24 to provide a list of allowed symbols (`allowed_symbols`), which 

25 must represent a subset of the elements that can occupy the sites 

26 on the specified sublattice. Note that additional arguments are 

27 required for the SGC and VCSGC ensembles, namely chemical 

28 potentials (`chemical_potentials`) for the former and constraint 

29 parameters (`phis` and `kappa`) for the latter. For more detailed 

30 information regarding the different ensembles, please see 

31 :class:`CanonicalEnsemble <mchammer.ensembles.CanonicalEnsemble>`, 

32 :class:`SemiGrandCanonicalEnsemble 

33 <mchammer.ensembles.SemiGrandCanonicalEnsemble>`, and 

34 :class:`VCSGCEnsemble <mchammer.ensembles.VCSGCEnsemble>`. 

35 

36 This class is particularly useful for effectively sampling complex 

37 multi-component systems with several active sublattices, in which 

38 case different ensembles can be defined for each of the latter. 

39 The fact that it is possible to set the allowed chemical symbols 

40 means that one can vary the concentrations of a few selected 

41 species, with the help of a VCSGC or semi-grand canonical 

42 ensemble, while still allowing swaps between any two sites, using 

43 a canonical ensemble (see also the below example). It is advisable 

44 to carefully consider how to define the ensemble probabilities. By 

45 default the ensembles are weighted by the sizes of the 

46 corresponding sublattices, which should give suitable 

47 probabilities in most cases. As is shown in the example below, it 

48 might be prudent to provide different values if allowed symbols 

49 are provided as well as for cases where there are several 

50 ensembles that are active on different sublattices. 

51 

52 Parameters 

53 ---------- 

54 

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

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

57 also defines the initial occupation vector 

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

59 calculator to be used for calculating the potential changes 

60 that enter the evaluation of the Metropolis criterion 

61 temperature : float 

62 temperature :math:`T` in appropriate units [commonly Kelvin] 

63 ensemble_specs: List[Dict] 

64 A list of dictionaries, which should contain the following items: 

65 

66 * 'ensemble', which could be either "vcsgc"; "semigrand" 

67 or "canonical", lowercase and uppercase letters or any 

68 combination thereof are accepted (required) 

69 * 'sublattice_index', index for the sublattice of 

70 interest (required) 

71 * 'allowed_symbols', list of allowed chemical symbols 

72 (default: read from ClusterSpace) 

73 * 'chemical_potentials', a dictionary of chemical 

74 potentials for each species 

75 :math:`\\mu_i`; the key denotes the species, the value 

76 :specifies the chemical potential in units that are 

77 :consistent with the underlying cluster expansion (only 

78 :applicable and required for SGC ensembles) 

79 * 'phis ', dictionary with average constraint parameters 

80 ':math:`\\phi_i`; the key denotes the species; for a 

81 N-component sublattice, there should be N - 1 

82 different `\\phi_i` (referred to as 

83 :math:`\\bar{\\phi}` in [SadErh12]_; only applicable 

84 and required for VCSGC ensembles, see also 

85 :class:`VCSGCEnsemble <mchammer.ensembles.VCSGCEnsemble>`) 

86 * 'kappa', parameter that constrains the variance of the 

87 'concentration (referred to as 

88 :math:`\\bar{\\kappa}` in [SadErh12]_; only applicable 

89 :and required for VCSGC ensembles) 

90 

91 probabilities: List[float] 

92 list of floats with the probabilities for choosing a 

93 particular ensemble with the same length as ensemble specs. 

94 If left unspecified the probabilties are weighted by the 

95 sizes of the associated sublattices 

96 boltzmann_constant : float 

97 Boltzmann constant :math:`k_B` in appropriate units, i.e. 

98 units that are consistent with the underlying cluster 

99 expansion and the temperature units [default: eV/K] 

100 user_tag : str 

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

102 random_seed : int 

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

104 simulation 

105 dc_filename : str 

106 name of file the data container associated with the 

107 ensemble will be written to; if the file 

108 exists it will be read, the data container will be appended, 

109 and the file will be updated/overwritten 

110 data_container_write_period : float 

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

112 written to file; writing periodically to file provides both 

113 a way to examine the progress of the simulation and to 

114 back up the data [default: 600 s] 

115 ensemble_data_write_interval : int 

116 interval at which data is written to the data container; 

117 this includes for example the current value of the 

118 calculator (i.e. usually the energy) as well as ensembles 

119 specific fields such as temperature or the number of atoms 

120 of different species 

121 trajectory_write_interval : int 

122 interval at which the current occupation vector of the 

123 atomic configuration is written to the data container. 

124 

125 Example 

126 ------- 

127 The following snippet illustrates how to carry out a simple Monte Carlo 

128 simulation using a combination of one canonical and one VCSGC ensemble. 

129 Specifically, the concentration of one species (Au) is kept constant 

130 while the others (Ag and Pd) are varied, while swaps are still allowed. 

131 Here, the parameters of the cluster expansion are set to emulate a simple 

132 Ising model in order to obtain an example that can be run without 

133 modification. In practice, one should of course use a proper cluster 

134 expansion:: 

135 

136 >>> from ase.build import bulk 

137 >>> from icet import ClusterExpansion, ClusterSpace 

138 >>> from mchammer.calculators import ClusterExpansionCalculator 

139 

140 >>> # prepare cluster expansion 

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

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

143 >>> # pairs are included) 

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

145 >>> cs = ClusterSpace(prim, cutoffs=[4.3], 

146 ... chemical_symbols=['Ag', 'Au', 'Pd']) 

147 >>> ce = ClusterExpansion( 

148 ... cs, [0, 0, 0, 0.1, 0.1, 0.1, -0.02, -0.02, -0.02]) 

149 

150 # define structure object 

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

152 >>> for i, atom in enumerate(structure): 

153 >>> if i % 2 == 0: 

154 >>> atom.symbol = 'Ag' 

155 >>> elif i % 3 == 0: 

156 >>> atom.symbol = 'Pd' 

157 

158 # the default probabilities for this case would be [0.5, 0.5], but 

159 # since the VCSGC ensemble only performs flips on a subset of all 

160 # sites on the sublattice, namely those originally occupied by Ag 

161 # and Pd atoms, specific values will be provided 

162 >>> weights = [len(structure), 

163 ... len([s for s in structure.get_chemical_symbols() if s != 'Au'])] 

164 >>> norm = sum(weights) 

165 >>> probabilities = [w / norm for w in weights] 

166 

167 # set up and run MC simulation 

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

169 >>> ensemble_specs = [ 

170 ... {'ensemble': 'canonical', 'sublattice_index': 0}, 

171 ... {'ensemble': 'vcsgc', 'sublattice_index': 0, 

172 ... 'phis': {'Ag': -0.2}, 'kappa': 200, 

173 ... 'allowed_symbols':['Ag', 'Pd']}] 

174 >>> mc = HybridEnsemble(structure=structure, calculator=calc, 

175 ... ensemble_specs=ensemble_specs, 

176 ... temperature=600, probabilities=probabilities, 

177 ... dc_filename='myrun_hybrid.dc') 

178 >>> mc.run(100) # carry out 100 trial steps 

179 """ 

180 

181 def __init__(self, 

182 structure: Atoms, 

183 calculator: BaseCalculator, 

184 temperature: float, 

185 ensemble_specs: List[Dict], 

186 probabilities: List[float] = None, 

187 boltzmann_constant: float = kB, 

188 user_tag: str = None, 

189 random_seed: int = None, 

190 dc_filename: str = None, 

191 data_container: str = None, 

192 data_container_write_period: float = 600, 

193 ensemble_data_write_interval: int = None, 

194 trajectory_write_interval: int = None) -> None: 

195 

196 # define available ensembles 

197 self._ensemble_trial_steps = dict([ 

198 ('canonical', self.do_canonical_swap), 

199 ('semi-grand', self.do_sgc_flip), 

200 ('vcsgc', self.do_vcsgc_flip), 

201 ]) 

202 

203 self._ensemble_parameters = dict(temperature=temperature) 

204 

205 self._trial_steps_per_ensemble = {"ensemble_{}".format(i): 0 for i in 

206 range(len(ensemble_specs))} 

207 

208 # process the list of ensembles and parameters 

209 self._process_ensemble_specs(ensemble_specs) 

210 

211 super().__init__( 

212 structure=structure, 

213 calculator=calculator, 

214 user_tag=user_tag, 

215 random_seed=random_seed, 

216 boltzmann_constant=boltzmann_constant, 

217 dc_filename=dc_filename, 

218 data_container=data_container, 

219 data_container_class=DataContainer, 

220 data_container_write_period=data_container_write_period, 

221 ensemble_data_write_interval=ensemble_data_write_interval, 

222 trajectory_write_interval=trajectory_write_interval) 

223 

224 # postprocess the list of ensembles and parameters 

225 self._postprocess_ensemble_args() 

226 

227 # set the probabilities 

228 self._process_probabilities(probabilities) 

229 

230 @property 

231 def temperature(self) -> float: 

232 """ Current temperature """ 

233 return self._ensemble_parameters['temperature'] 

234 

235 @property 

236 def probabilities(self) -> List[float]: 

237 """ Ensemble propabilities """ 

238 return self._probabilities 

239 

240 @property 

241 def trial_steps_per_ensemble(self) -> Dict[str, int]: 

242 """ Number of Monte Carlo trial steps for each ensemble """ 

243 return self._trial_steps_per_ensemble 

244 

245 def _process_ensemble_specs(self, ensemble_specs: List[Dict]) -> None: 

246 """Process the list of ensembles and parameters 

247 

248 Parameters 

249 ---------- 

250 ensemble_specs: List[Dict] 

251 A list of dictionaries, which should contain the following items: 

252 * 'ensemble', which could be either "vcsgc"; "semigrand" or "canonical", lowercase and 

253 upercase letters or any combination thereof are accepted 

254 * 'sublattice_index', index for the sublattice of interest 

255 * 'allowed_symbols', list of allowed chemical symbols 

256 * 'chemical_potentials', a dictionary of chemical potentials for each species 

257 :math:`\\mu_i`; the key denotes the species, the value specifies the chemical potential 

258 in units that are consistent with the underlying cluster expansion 

259 * 'phis ', dictionary with average constraint parameters :math:`\\phi_i`; the key 

260 denotes the species; for a N-component sublattice, there should be N - 1 different 

261 `\\phi_i` 

262 * 'kappa', parameter that constrains the variance of the concentration 

263 """ 

264 

265 ensemble_args = [] 

266 

267 for ind, ensemble_spec in enumerate(ensemble_specs): 

268 

269 ensemble_arg = {} # type: Dict[str, Any] 

270 tag = "ensemble_{}".format(ind) 

271 ensemble_arg['tag'] = tag 

272 

273 # check the ensemble name 

274 if 'ensemble' not in ensemble_spec: 

275 raise ValueError("The dictionary {} lacks the required key" 

276 " 'ensemble'".format(ensemble_spec)) 

277 ensemble = ensemble_spec['ensemble'].lower() 

278 if ensemble not in self._ensemble_trial_steps.keys(): 

279 msg = ['Ensemble not available'] 

280 msg += ['Please choose one of the following:'] 

281 for key in self._ensemble_trial_steps.keys(): 

282 msg += [' * ' + key] 

283 raise ValueError('\n'.join(msg)) 

284 ensemble_arg['ensemble'] = ensemble 

285 self._ensemble_parameters[tag] = ensemble 

286 

287 # check that all required keys, and no unknown keys, are present 

288 keys = ['ensemble', 'sublattice_index', 'allowed_symbols'] 

289 if ensemble == 'semi-grand': 

290 keys = ['chemical_potentials'] + keys 

291 elif ensemble == 'vcsgc': 

292 keys = ['phis', 'kappa'] + keys 

293 for key in keys[:-1]: 

294 if key not in ensemble_spec: 

295 raise ValueError("The dictionary {} lacks the key '{}', which is required for" 

296 " {} ensembles".format(ensemble_spec, key, ensemble)) 

297 for key in ensemble_spec.keys(): 

298 if key not in keys: 

299 raise ValueError("Unknown key '{}', for a {} ensemble," 

300 " in the dictionary {}".format(key, ensemble, ensemble_spec)) 

301 

302 # record the sublattice index 

303 ensemble_arg['sublattice_index'] = ensemble_spec['sublattice_index'] 

304 

305 # process chemical potentials 

306 if 'chemical_potentials' in ensemble_spec: 

307 chemical_potentials = get_chemical_potentials(ensemble_spec['chemical_potentials']) 

308 ensemble_arg['chemical_potentials'] = chemical_potentials 

309 for atnum, chempot in chemical_potentials.items(): 

310 mu_sym = '{}_mu_{}'.format(tag, chemical_symbols[atnum]) 

311 self._ensemble_parameters[mu_sym] = chempot 

312 

313 # process phis 

314 if 'phis' in ensemble_spec: 

315 phis = get_phis(ensemble_spec['phis']) 

316 ensemble_arg['phis'] = phis 

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

318 if isinstance(sym, str): 318 ↛ 319line 318 didn't jump to line 319, because the condition on line 318 was never true

319 chemical_symbol = sym 

320 else: 

321 chemical_symbol = chemical_symbols[sym] 

322 phi_sym = '{}_phi_{}'.format(tag, chemical_symbol) 

323 self._ensemble_parameters[phi_sym] = phi 

324 

325 # process kappa 

326 if 'kappa' in ensemble_spec: 

327 ensemble_arg['kappa'] = ensemble_spec['kappa'] 

328 self._ensemble_parameters['{}_kappa'.format(tag)] = ensemble_spec['kappa'] 

329 

330 # record the allowed chemical symbols 

331 if 'allowed_symbols' in ensemble_spec: 

332 ensemble_arg['allowed_symbols'] = ensemble_spec['allowed_symbols'] 

333 

334 ensemble_args.append(ensemble_arg) 

335 

336 self._ensemble_args = ensemble_args 

337 

338 def _postprocess_ensemble_args(self): 

339 """Process the list of dictionaries with ensemble specific parameters 

340 """ 

341 

342 for i in range(len(self._ensemble_args)): 

343 

344 # check the sublattice index 

345 self._check_sublattice_index(self._ensemble_args[i]['sublattice_index']) 

346 

347 # extract the allowed species 

348 if 'allowed_symbols' in self._ensemble_args[i]: 

349 self._ensemble_args[i]['allowed_species'] =\ 

350 self._extract_allowed_species(self._ensemble_args[i]['allowed_symbols'], 

351 self._ensemble_args[i]['sublattice_index']) 

352 del self._ensemble_args[i]['allowed_symbols'] 

353 else: 

354 self._ensemble_args[i]['allowed_species'] = None 

355 

356 if self._ensemble_args[i]['ensemble'] == 'vcsgc': 

357 # Check that each sublattice has N - 1 phis 

358 count_specified_elements = 0 

359 if self._ensemble_args[i]['allowed_species'] is None: 359 ↛ 363line 359 didn't jump to line 363, because the condition on line 359 was never false

360 allowed_species =\ 

361 self.sublattices[self._ensemble_args[i]['sublattice_index']].atomic_numbers 

362 else: 

363 allowed_species = self._ensemble_args[i]['allowed_species'] 

364 for number in allowed_species: 

365 if number in self._ensemble_args[i]['phis'].keys(): 

366 count_specified_elements += 1 

367 if count_specified_elements != len(allowed_species) - 1: 

368 raise ValueError("phis must be set for N - 1 elements on a sublattice with" 

369 " N elements") 

370 

371 def _check_sublattice_index(self, sublattice_index: int): 

372 """Check the 'sublattice_index' item in the 'ensemble_spec' dictionary 

373 

374 Parameters 

375 ---------- 

376 sublattice_index: 

377 Specific sublattice to consider provided as as an index or a symbol 

378 """ 

379 

380 if not isinstance(sublattice_index, int): 

381 raise TypeError("'sublattice_index' must be an integer, not" 

382 " {}".format(type(sublattice_index))) 

383 

384 # check that the sublattice exists 

385 if sublattice_index not in range(len(self.sublattices)): 

386 raise ValueError("There is no sublattice with index {}".format(sublattice_index)) 

387 

388 # check that the sublattice is active 

389 if len(self.sublattices[sublattice_index].chemical_symbols) == 1: 

390 raise ValueError("The sublattice {} is inactive".format(sublattice_index)) 

391 

392 def _extract_allowed_species(self, allowed_symbols: List[str], sublattice_index: int 

393 ) -> List[int]: 

394 """Check and extract the allowed species from the 'allowed_symbols' in the 'ensemble_spec' 

395 dictionary 

396 

397 Parameters 

398 ---------- 

399 allowed_symbols: 

400 list of allowed chemical symbols 

401 sublattice_index: 

402 Index for the relevant sublattice 

403 """ 

404 

405 if not isinstance(allowed_symbols, list) or not all( 

406 [isinstance(i, str) for i in allowed_symbols]): 

407 raise TypeError( 

408 "'allowed_symbols' must be a List[str], not {}".format(type(allowed_symbols))) 

409 for symbol in allowed_symbols: 

410 if symbol not in self.sublattices[sublattice_index].chemical_symbols: 

411 raise ValueError("The species {} is not allowed on sublattice" 

412 " {}".format(symbol, sublattice_index)) 

413 

414 return [atomic_numbers[s] for s in allowed_symbols] 

415 

416 def _process_probabilities(self, probabilities: List[float]): 

417 """Process the list of probabilities 

418 

419 Parameters 

420 ---------- 

421 probabilities: 

422 list of floats with the probabilities for choosing a particular ensemble with the same 

423 length as self._ensemble_args. 

424 """ 

425 

426 if probabilities is None: 

427 # use the sizes of the associated sublattices when calculating the ensemble 

428 # probabilities 

429 weights = [len(self.sublattices[ensemble_arg['sublattice_index']].indices) for 

430 ensemble_arg in self._ensemble_args] 

431 norm = sum(weights) 

432 probabilities = [w / norm for w in weights] 

433 else: 

434 if len(probabilities) != len(self._ensemble_args): 

435 raise ValueError("The number of probabilities must be match the number of" 

436 " ensembles") 

437 

438 if not isclose(sum(probabilities), 1.0): 

439 raise ValueError("The sum of all probabilities must be equal to 1") 

440 

441 self._probabilities = probabilities 

442 

443 def _do_trial_step(self): 

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

445 

446 # randomly pick an ensemble 

447 ensemble_arg = np.random.choice(self._ensemble_args, p=self._probabilities) 

448 

449 # count number of trial steps for each ensemble 

450 self._trial_steps_per_ensemble[ensemble_arg['tag']] += 1 

451 

452 if ensemble_arg['ensemble'] == 'canonical' and not self.configuration.is_swap_possible( 452 ↛ 454line 452 didn't jump to line 454, because the condition on line 452 was never true

453 ensemble_arg['sublattice_index'], ensemble_arg['allowed_species']): 

454 return 0 

455 else: 

456 arguments = {key: val for key, val in ensemble_arg.items() if key not in 

457 ['ensemble', 'tag']} 

458 return self._ensemble_trial_steps[ensemble_arg['ensemble']](**arguments) 

459 

460 def _get_ensemble_data(self) -> Dict: 

461 """ 

462 Returns a dict with the default data of the ensemble. This includes 

463 atom counts and free energy derivative. 

464 """ 

465 data = super()._get_ensemble_data() 

466 

467 ensemble_types = [e['ensemble'] for e in self._ensemble_args] 

468 

469 # free energy derivative 

470 if 'vcsgc' in ensemble_types: 470 ↛ 478line 470 didn't jump to line 478, because the condition on line 470 was never false

471 for ensemble_arg in self._ensemble_args: 

472 if 'vcsgc' == ensemble_arg['ensemble']: 

473 data.update(self._get_vcsgc_free_energy_derivatives( 

474 ensemble_arg['phis'], ensemble_arg['kappa'], 

475 ensemble_arg['sublattice_index'])) 

476 

477 # species counts 

478 if any([e in ensemble_types for e in ['vcsgc', 'semi-grand']]): 478 ↛ 481line 478 didn't jump to line 481, because the condition on line 478 was never false

479 data.update(self._get_species_counts()) 

480 

481 return data