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"""Definition of the Wang-Landau algorithm class.""" 

2 

3import random 

4 

5from collections import OrderedDict 

6from typing import Any, Dict, List, Optional, Tuple, Union 

7 

8import numpy as np 

9 

10from ase import Atoms 

11 

12from .. import WangLandauDataContainer 

13from ..calculators.base_calculator import BaseCalculator 

14from .thermodynamic_base_ensemble import BaseEnsemble 

15from icet.input_output.logging_tools import logger 

16 

17logger = logger.getChild('wang_landau_ensemble') 

18 

19 

20class WangLandauEnsemble(BaseEnsemble): 

21 """Instances of this class allow one to sample a system using the 

22 Wang-Landau (WL) algorithm, see Phys. Rev. Lett. **86**, 2050 

23 (2001) [WanLan01a]_. The WL algorithm enables one to acquire the 

24 density of states (DOS) as a function of energy, from which one 

25 can readily calculate many thermodynamic observables as a function 

26 of temperature. To this end, the WL algorithm accumulates both the 

27 microcanonical entropy :math:`S(E)` and a histogram :math:`H(E)` 

28 on an energy grid with a predefined spacing (``energy_spacing``). 

29 

30 The algorithm is initialized as follows. 

31 

32 #. Generate an initial configuration. 

33 #. Initialize counters for the microcanonical entropy 

34 :math:`S(E)` and the histogram :math:`H(E)` to zero. 

35 #. Set the fill factor :math:`f=1`. 

36 

37 It then proceeds as follows. 

38 

39 #. Propose a new configuration (see ``trial_move``). 

40 #. Accept or reject the new configuration with probability 

41 

42 .. math:: 

43 

44 P = \\min \\{ 1, \\, \\exp [ S(E_\\mathrm{new}) - S(E_\\mathrm{cur}) ] \\}, 

45 

46 where :math:`E_\\mathrm{cur}` and :math:`E_\\mathrm{new}` are the 

47 energies of the current and new configurations, respectively. 

48 #. Update the microcanonical entropy :math:`S(E)\\leftarrow S(E) + f` 

49 and histogram :math:`H(E) \\leftarrow H(E) + 1` where 

50 :math:`E` is the energy of the system at the end of the move. 

51 #. Check the flatness of the histogram :math:`H(E)`. If 

52 :math:`H(E) > \\chi \\langle H(E)\\rangle\\,\\forall E` reset the histogram 

53 :math:`H(E) = 0` and reduce the fill factor :math:`f \\leftarrow f / 2`. 

54 The parameter :math:`\\chi` is set via ``flatness_limit``. 

55 #. If :math:`f` is smaller than ``fill_factor_limit`` terminate 

56 the loop, otherwise return to 1. 

57 

58 The microcanonical entropy :math:`S(E)` and the histogram along 

59 with related information are written to the data container every 

60 time :math:`f` is updated. Using the density :math:`\\rho(E) = \\exp S(E)` 

61 one can then readily compute various thermodynamic quantities, 

62 including, e.g., the average energy: 

63 

64 .. math:: 

65 

66 \\left<E\\right> = \\frac{\\sum_E E \\rho(E) \\exp(-E / k_B T)}{ 

67 \\sum_E \\rho(E) \\exp(-E / k_B T)} 

68 

69 Parameters 

70 ---------- 

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

72 atomic configuration to be used in the Wang-Landau simulation; 

73 also defines the initial occupation vector 

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

75 calculator to be used for calculating potential changes 

76 trial_move : str 

77 One can choose between two different trial moves for 

78 generating new configurations. In a 'swap' move two sites are 

79 selected and their occupations are swapped; in a 'flip' move 

80 one site is selected and its occupation is flipped to a 

81 different species. While 'swap' moves conserve the 

82 concentrations of the species in the system, 'flip' moves 

83 allow one in principle to sample the full composition space. 

84 energy_spacing : float 

85 defines the bin size of the energy grid on which the microcanonical 

86 entropy :math:`S(E)`, and thus the density :math:`\\exp S(E)`, is 

87 evaluated; the spacing should be small enough to capture the features 

88 of the density of states; too small values will, however, render the 

89 convergence very tedious if not impossible 

90 energy_limit_left : float 

91 defines the lower limit of the energy range within which the 

92 microcanonical entropy :math:`S(E)` will be sampled. By default 

93 (`None`) no limit is imposed. Setting limits can be useful if only a 

94 part of the density of states is required. 

95 energy_limit_right : float 

96 defines the upper limit of the energy range within which the 

97 microcanonical entropy :math:`S(E)` will be sampled. By default 

98 (`None`) no limit is imposed. Setting limits can be useful if only a 

99 part of the density of states is required. 

100 fill_factor_limit : float 

101 If the fill_factor :math:`f` falls below this value, the 

102 WL sampling loop is terminated. 

103 flatness_check_interval : int 

104 For computational efficiency the flatness condition is only 

105 evaluated every ``flatness_check_interval``-th trial step. By 

106 default (``None``) ``flatness_check_interval`` is set to 1000 

107 times the number of sites in ``structure``, i.e. 1000 Monte 

108 Carlo sweeps. 

109 flatness_limit : float 

110 The histogram :math:`H(E)` is deemed sufficiently flat if 

111 :math:`H(E) > \\chi \\left<H(E)\\right>\\,\\forall 

112 E`. ``flatness_limit`` sets the parameter :math:`\\chi`. 

113 window_search_penalty : float 

114 If `energy_limit_left` and/or `energy_limit_right` have been 

115 provided, a modified acceptance probability, 

116 :math:`P=\\min\\{1,\\,\\exp[C_\\mathrm{WSP}(d_\\mathrm{new}- 

117 d_\\mathrm{cur})]\\}`, will be used until a configuration is 

118 found within the interval of interest. This parameter, 

119 specifically, corresponds to :math:`C_\\mathrm{WSP}`, which 

120 controls how strongly moves that lead to an increase in the 

121 distance, i.e. difference in energy divided by the energy 

122 spacing, to the energy window (:math:`d_\\mathrm{new}> 

123 d_\\mathrm{cur}`) should be penalized. A higher value leads 

124 to a lower acceptance probability for such moves. 

125 user_tag : str 

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

127 dc_filename : str 

128 name of file the data container associated with the ensemble 

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

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

131 updated/overwritten 

132 random_seed : int 

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

134 simulation 

135 ensemble_data_write_interval : int 

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

137 includes for example the current value of the calculator 

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

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

140 data_container_write_period : float 

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

142 written to file; writing periodically to file provides both 

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

144 the data [default: 600 s] 

145 trajectory_write_interval : int 

146 interval at which the current occupation vector of the atomic 

147 configuration is written to the data container. 

148 sublattice_probabilities : List[float] 

149 probability for picking a sublattice when doing a random swap. 

150 The list must contain as many elements as there are sublattices 

151 and it needs to sum up to 1. 

152 

153 Example 

154 ------- 

155 The following snippet illustrates how to carry out a Wang-Landau 

156 simulation. For the purpose of demonstration, the parameters of 

157 the cluster expansion are set to obtain a two-dimensional square 

158 Ising model, one of the systems studied in the original work by 

159 Wang and Landau:: 

160 

161 >>> from ase import Atoms 

162 >>> from icet import ClusterExpansion, ClusterSpace 

163 >>> from mchammer.calculators import ClusterExpansionCalculator 

164 >>> from mchammer.ensembles import WangLandauEnsemble 

165 

166 >>> # prepare cluster expansion 

167 >>> prim = Atoms('Au', positions=[[0, 0, 0]], cell=[1, 1, 10], pbc=True) 

168 >>> cs = ClusterSpace(prim, cutoffs=[1.1], chemical_symbols=['Ag', 'Au']) 

169 >>> ce = ClusterExpansion(cs, [0, 0, 2]) 

170 

171 >>> # prepare initial configuration 

172 >>> structure = prim.repeat((4, 4, 1)) 

173 >>> for k in range(8): 

174 ... structure[k].symbol = 'Ag' 

175 

176 >>> # set up and run Wang-Landau simulation 

177 >>> calculator = ClusterExpansionCalculator(structure, ce) 

178 >>> mc = WangLandauEnsemble(structure=structure, 

179 ... calculator=calculator, 

180 ... energy_spacing=1, 

181 ... dc_filename='ising_2d_run.dc') 

182 >>> mc.run(number_of_trial_steps=len(structure)*1000) # in practice one requires more steps 

183 

184 """ 

185 

186 def __init__(self, 

187 structure: Atoms, 

188 calculator: BaseCalculator, 

189 energy_spacing: float, 

190 energy_limit_left: float = None, 

191 energy_limit_right: float = None, 

192 trial_move: str = 'swap', 

193 fill_factor_limit: float = 1e-6, 

194 flatness_check_interval: int = None, 

195 flatness_limit: float = 0.8, 

196 window_search_penalty: float = 2.0, 

197 user_tag: str = None, 

198 dc_filename: str = None, 

199 data_container: str = None, 

200 random_seed: int = None, 

201 data_container_write_period: float = 600, 

202 ensemble_data_write_interval: int = None, 

203 trajectory_write_interval: int = None, 

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

205 

206 # set trial move 

207 if trial_move == 'swap': 

208 self.do_move = self._do_swap 

209 self._get_sublattice_probabilities = self._get_swap_sublattice_probabilities 

210 elif trial_move == 'flip': 

211 self.do_move = self._do_flip 

212 self._get_sublattice_probabilities = self._get_flip_sublattice_probabilities 

213 else: 

214 raise ValueError('Invalid value for trial_move: {}.' 

215 ' Must be either "swap" or "flip".'.format(trial_move)) 

216 

217 # set default values that are system dependent 

218 if flatness_check_interval is None: 

219 flatness_check_interval = len(structure) * 1000 

220 

221 # parameters pertaining to construction of entropy and histogram 

222 self._energy_spacing = energy_spacing 

223 self._fill_factor_limit = fill_factor_limit 

224 self._flatness_check_interval = flatness_check_interval 

225 self._flatness_limit = flatness_limit 

226 

227 # energy window 

228 self._window_search_penalty = window_search_penalty 

229 self._bin_left = self._get_bin_index(energy_limit_left) 

230 self._bin_right = self._get_bin_index(energy_limit_right) 

231 if self._bin_left is not None and \ 

232 self._bin_right is not None and self._bin_left >= self._bin_right: 

233 raise ValueError('Invalid energy window: left boundary ({}, {}) must be' 

234 ' smaller than right boundary ({}, {})' 

235 .format(energy_limit_left, self._bin_left, 

236 energy_limit_right, self._bin_right)) 

237 

238 # ensemble parameters 

239 self._ensemble_parameters = {} 

240 self._ensemble_parameters['energy_spacing'] = energy_spacing 

241 self._ensemble_parameters['trial_move'] = trial_move 

242 self._ensemble_parameters['energy_limit_left'] = energy_limit_left 

243 self._ensemble_parameters['energy_limit_right'] = energy_limit_right 

244 # The following parameters are _intentionally excluded_ from 

245 # the ensemble_parameters dict as it would prevent users from 

246 # changing their values between restarts. The latter is advantageous 

247 # as these runs can require restarts and possibly parameter adjustments 

248 # to achieve convergence. 

249 # * fill_factor_limit 

250 # * flatness_check_interval 

251 # * flatness_limit 

252 # * entropy_write_frequency 

253 # * window_search_penalty 

254 

255 # add species count to ensemble parameters 

256 symbols = set([symbol for sub in calculator.sublattices 

257 for symbol in sub.chemical_symbols]) 

258 for symbol in symbols: 

259 key = 'n_atoms_{}'.format(symbol) 

260 count = structure.get_chemical_symbols().count(symbol) 

261 self._ensemble_parameters[key] = count 

262 

263 # set the convergence, which may be updated in case of a restart 

264 self._converged = None # type: Optional[bool] 

265 

266 # the constructor of the parent classes must be called *after* 

267 # the ensemble_parameters dict has been populated 

268 super().__init__( 

269 structure=structure, 

270 calculator=calculator, 

271 user_tag=user_tag, 

272 random_seed=random_seed, 

273 dc_filename=dc_filename, 

274 data_container=data_container, 

275 data_container_class=WangLandauDataContainer, 

276 data_container_write_period=data_container_write_period, 

277 ensemble_data_write_interval=ensemble_data_write_interval, 

278 trajectory_write_interval=trajectory_write_interval) 

279 

280 # handle probabilities for swaps on different sublattices 

281 if sublattice_probabilities is None: 

282 self._sublattice_probabilities = self._get_sublattice_probabilities() 

283 else: 

284 self._sublattice_probabilities = sublattice_probabilities 

285 

286 # initialize Wang-Landau algorithm; in the case of a restart 

287 # these quantities are read from the data container file; the 

288 # if-conditions prevent these values from being overwritten 

289 self._potential = self.calculator.calculate_total( 

290 occupations=self.configuration.occupations) 

291 self._reached_energy_window = self._inside_energy_window( 

292 self._get_bin_index(self._potential)) 

293 if not hasattr(self, '_fill_factor'): 

294 self._fill_factor = 1.0 

295 if not hasattr(self, '_fill_factor_history'): 

296 if self._reached_energy_window: 

297 self._fill_factor_history = {self.step: self._fill_factor} 

298 else: 

299 self._fill_factor_history = {} 

300 if not hasattr(self, '_entropy_history'): 

301 self._entropy_history = {} 

302 if not hasattr(self, '_histogram'): 

303 self._histogram = {} # type: Dict[int, int] 

304 if not hasattr(self, '_entropy'): 

305 self._entropy = {} # type: Dict[int, float] 

306 

307 @property 

308 def fill_factor(self) -> float: 

309 """ current value of fill factor """ 

310 return self._fill_factor 

311 

312 @property 

313 def fill_factor_history(self) -> Dict[int, float]: 

314 """evolution of the fill factor in the Wang-Landau algorithm (key=MC 

315 trial step, value=fill factor) 

316 """ 

317 return self._fill_factor_history 

318 

319 @property 

320 def converged(self) -> Optional[bool]: 

321 """ True if convergence has been achieved """ 

322 return self._converged 

323 

324 @property 

325 def flatness_limit(self) -> float: 

326 """The histogram :math:`H(E)` is deemed sufficiently flat if 

327 :math:`H(E) > \\chi \\left<H(E)\\right>\\,\\forall 

328 E` where ``flatness_limit`` sets the parameter :math:`\\chi`. 

329 """ 

330 return self._flatness_limit 

331 

332 @flatness_limit.setter 

333 def flatness_limit(self, new_value) -> None: 

334 self._flatness_limit = new_value 

335 self._converged = None 

336 

337 @property 

338 def fill_factor_limit(self) -> float: 

339 """ If the fill_factor :math:`f` falls below this value, the 

340 Wang-Landau sampling loop is terminated. """ 

341 return self._fill_factor_limit 

342 

343 @fill_factor_limit.setter 

344 def fill_factor_limit(self, new_value) -> None: 

345 self._fill_factor_limit = new_value 

346 self._converged = None 

347 

348 @property 

349 def flatness_check_interval(self) -> int: 

350 """ number of MC trial steps between checking the flatness 

351 condition """ 

352 return self._flatness_check_interval 

353 

354 @flatness_check_interval.setter 

355 def flatness_check_interval(self, new_value: int) -> None: 

356 self._flatness_check_interval = new_value 

357 

358 def run(self, number_of_trial_steps: int): 

359 """ 

360 Samples the ensemble for the given number of trial steps. 

361 

362 Parameters 

363 ---------- 

364 number_of_trial_steps 

365 maximum number of MC trial steps to run in total (the 

366 run will terminate earlier if `fill_factor_limit` is reached) 

367 reset_step 

368 if True the MC trial step counter and the data container will 

369 be reset to zero and empty, respectively. 

370 

371 Raises 

372 ------ 

373 TypeError 

374 if `number_of_trial_steps` is not an int 

375 """ 

376 if self.converged: 

377 logger.warning('Convergence has already been reached.') 

378 else: 

379 super().run(number_of_trial_steps) 

380 

381 def _terminate_sampling(self) -> bool: 

382 """Returns True if the Wang-Landau algorithm has converged. This is 

383 used in the run method implemented of BaseEnsemble to 

384 evaluate whether the sampling loop should be terminated. 

385 """ 

386 # N.B.: self._converged can be None 

387 if self._converged is not None: 

388 return self._converged 

389 else: 

390 return False 

391 

392 def _restart_ensemble(self): 

393 """Restarts ensemble using the last state saved in the data container 

394 file. Note that this method does _not_ use the last_state property of 

395 the data container but rather uses the last data written the data frame. 

396 """ 

397 super()._restart_ensemble() 

398 self._fill_factor = self.data_container._last_state['fill_factor'] 

399 self._fill_factor_history = self.data_container._last_state['fill_factor_history'] 

400 self._entropy_history = self.data_container._last_state['entropy_history'] 

401 self._histogram = self.data_container._last_state['histogram'] 

402 self._entropy = self.data_container._last_state['entropy'] 

403 histogram = np.array(list(self._histogram.values())) 

404 limit = self._flatness_limit * np.average(histogram) 

405 self._converged = (self._fill_factor <= self._fill_factor_limit 

406 ) & np.all(histogram >= limit) 

407 

408 def write_data_container(self, outfile: Union[str, bytes]): 

409 """Updates last state of the Wang-Landau simulation and 

410 writes DataContainer to file. 

411 

412 Parameters 

413 ---------- 

414 outfile 

415 file to which to write 

416 """ 

417 self._data_container._update_last_state( 

418 last_step=self.step, 

419 occupations=self.configuration.occupations.tolist(), 

420 accepted_trials=self._accepted_trials, 

421 random_state=random.getstate(), 

422 fill_factor=self._fill_factor, 

423 fill_factor_history=self._fill_factor_history, 

424 entropy_history=self._entropy_history, 

425 histogram=OrderedDict(sorted(self._histogram.items())), 

426 entropy=OrderedDict(sorted(self._entropy.items()))) 

427 self.data_container.write(outfile) 

428 

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

430 """Evaluates Metropolis acceptance criterion. 

431 

432 Parameters 

433 ---------- 

434 potential_diff 

435 change in the thermodynamic potential associated 

436 with the trial step 

437 """ 

438 

439 # acceptance/rejection step 

440 bin_old = self._get_bin_index(self._potential) 

441 bin_new = self._get_bin_index(self._potential + potential_diff) 

442 bin_cur = bin_old 

443 if self._allow_move(bin_cur, bin_new): 

444 S_cur = self._entropy.get(bin_cur, 0) 

445 S_new = self._entropy.get(bin_new, 0) 

446 delta = np.exp(S_cur - S_new) 

447 if delta >= 1 or delta >= self._next_random_number(): 

448 accept = True 

449 self._potential += potential_diff 

450 bin_cur = bin_new 

451 else: 

452 accept = False 

453 else: 

454 accept = False 

455 

456 if not self._reached_energy_window: 

457 # check whether the target energy window has been reached 

458 self._reached_energy_window = self._inside_energy_window(bin_cur) 

459 # if the target window has been reached remove unused bins 

460 # from histogram and entropy counters 

461 if self._reached_energy_window: 

462 self._fill_factor_history[self.step] = self._fill_factor 

463 # flush data from data container except for initial step 

464 self._data_container._data_list = [self._data_container._data_list[0]] 

465 self._entropy = {k: self._entropy[k] 

466 for k in self._entropy if self._inside_energy_window(k)} 

467 self._histogram = {k: self._histogram[k] 

468 for k in self._histogram if self._inside_energy_window(k)} 

469 else: 

470 # then reconsider accept/reject based on whether we 

471 # approached the window or not 

472 dist_new = np.inf 

473 dist_old = np.inf 

474 if self._bin_left is not None: 474 ↛ 477line 474 didn't jump to line 477, because the condition on line 474 was never false

475 dist_new = min(dist_new, abs(bin_new - self._bin_left)) 

476 dist_old = min(dist_old, abs(bin_old - self._bin_left)) 

477 if self._bin_right is not None: 

478 dist_new = min(dist_new, abs(bin_new - self._bin_right)) 

479 dist_old = min(dist_old, abs(bin_old - self._bin_right)) 

480 assert dist_new < np.inf and dist_old < np.inf 

481 exp_dist = np.exp((dist_old - dist_new) * self._window_search_penalty) 

482 if exp_dist >= 1 or exp_dist >= self._next_random_number(): 

483 # should be accepted 

484 if not accept: 

485 # reset potential 

486 self._potential += potential_diff 

487 bin_cur = bin_new 

488 accept = True 

489 else: 

490 # should be rejected 

491 if accept: 491 ↛ 494line 491 didn't jump to line 494, because the condition on line 491 was never false

492 # reset potential 

493 self._potential -= potential_diff 

494 bin_cur = bin_old 

495 accept = False 

496 

497 # update histograms and entropy counters 

498 self._update_entropy(bin_cur) 

499 

500 return accept 

501 

502 def _update_entropy(self, bin_cur: int) -> None: 

503 """Updates counters for histogram and entropy, checks histogram 

504 flatness, and updates fill factor if indicated. 

505 """ 

506 

507 # update histogram and entropy 

508 self._entropy[bin_cur] = self._entropy.get(bin_cur, 0) + self._fill_factor 

509 self._histogram[bin_cur] = self._histogram.get(bin_cur, 0) + 1 

510 

511 # check flatness of histogram 

512 if self.step % self._flatness_check_interval == 0 and \ 

513 self.step > 0 and self._reached_energy_window: 

514 

515 # shift entropy counter in order to avoid overflow 

516 entropy_ref = np.min(list(self._entropy.values())) 

517 for k in self._entropy: 

518 self._entropy[k] -= entropy_ref 

519 

520 # check whether the Wang-Landau algorithm has converged 

521 histogram = np.array(list(self._histogram.values())) 

522 limit = self._flatness_limit * np.average(histogram) 

523 is_flat = np.all(histogram >= limit) 

524 self._converged = (self._fill_factor <= self._fill_factor_limit) & is_flat 

525 if is_flat and not self._converged: 

526 # update fill factor 

527 self._fill_factor /= 2 

528 self._fill_factor_history[self.step] = self._fill_factor 

529 # update entropy history 

530 self._entropy_history[self.step] = OrderedDict( 

531 sorted(self._entropy.items())) 

532 # reset histogram 

533 self._histogram = dict.fromkeys(self._histogram, 0) 

534 

535 def _get_bin_index(self, energy: float) -> Optional[int]: 

536 """ Returns bin index for histogram and entropy dictionaries. """ 

537 if energy is None or np.isnan(energy): 

538 return None 

539 return int(np.around(energy / self._energy_spacing)) 

540 

541 def _allow_move(self, bin_cur: Optional[int], bin_new: int) -> bool: 

542 """Returns True if the current move is to be included in the 

543 accumulation of histogram and entropy. This logic has been 

544 moved into a separate function in order to enhance 

545 readability. 

546 """ 

547 if self._bin_left is None and self._bin_right is None: 

548 # no limits on energy window 

549 return True 

550 if self._bin_left is not None: 

551 if bin_cur < self._bin_left: 

552 # not yet in window (left limit) 

553 return True 

554 if bin_new < self._bin_left: 

555 # imposing left limit 

556 return False 

557 if self._bin_right is not None: 

558 if bin_cur > self._bin_right: 

559 # not yet in window (right limit) 

560 return True 

561 if bin_new > self._bin_right: 

562 # imposing right limit 

563 return False 

564 return True 

565 

566 def _inside_energy_window(self, bin_k: int) -> bool: 

567 """Returns True if bin_k is inside the energy window specified for 

568 this simulation. 

569 """ 

570 if self._bin_left is not None and bin_k < self._bin_left: 

571 return False 

572 if self._bin_right is not None and bin_k > self._bin_right: 

573 return False 

574 return True 

575 

576 def _do_trial_step(self): 

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

578 sublattice_index = self.get_random_sublattice_index(self._sublattice_probabilities) 

579 return self.do_move(sublattice_index=sublattice_index) 

580 

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

582 """Carries out a Monte Carlo trial that involves swapping the species 

583 on two sites. This method has been copied from 

584 ThermodynamicBaseEnsemble. 

585 

586 Parameters 

587 --------- 

588 sublattice_index 

589 the sublattice the swap will act on 

590 allowed_species 

591 list of atomic numbers for allowed species 

592 

593 Returns 

594 ------- 

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

596 """ 

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

598 potential_diff = self._get_property_change(sites, species) 

599 if self._acceptance_condition(potential_diff): 

600 self.update_occupations(sites, species) 

601 return 1 

602 return 0 

603 

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

605 """Carries out one Monte Carlo trial step that involves flipping the 

606 species on one site. This method has been adapted from 

607 ThermodynamicBaseEnsemble. 

608 

609 Parameters 

610 --------- 

611 sublattice_index 

612 the sublattice the flip will act on 

613 allowed_species 

614 list of atomic numbers for allowed species 

615 

616 Returns 

617 ------- 

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

619 

620 """ 

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

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

623 if self._acceptance_condition(potential_diff): 

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

625 return 1 

626 return 0 

627 

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

629 """Returns sublattice probabilities suitable for swaps. This method 

630 has been copied without modification from ThermodynamicBaseEnsemble. 

631 """ 

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

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

634 if self.configuration.is_swap_possible(i): 

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

636 else: 

637 sublattice_probabilities.append(0) 

638 norm = sum(sublattice_probabilities) 

639 if norm == 0: 

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

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

642 return sublattice_probabilities 

643 

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

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

646 sizes of a sublattice. This method has been copied without 

647 modification from ThermodynamicBaseEnsemble. 

648 """ 

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

650 for _, sl in enumerate(self.sublattices): 

651 if len(sl.chemical_symbols) > 1: 

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

653 else: 

654 sublattice_probabilities.append(0) 

655 norm = sum(sublattice_probabilities) 

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

657 return sublattice_probabilities 

658 

659 

660def get_bins_for_parallel_simulations(n_bins: int, 

661 energy_spacing: float, 

662 minimum_energy: float, 

663 maximum_energy: float, 

664 overlap: int = 4, 

665 bin_size_exponent: float = 1.0) -> List[Tuple[float, float]]: 

666 """Generates a list of energy bins (lower and upper bound) suitable for 

667 parallel Wang-Landau simulations. For the latter, the energy range is 

668 split up into a several bins (``n_bins``). Each bin is then sampled in a 

669 separate Wang-Landau simulation. Once the density of states in the 

670 individual bins has been converged the total density of states can be 

671 constructed by patching the segments back together. To this end, one 

672 requires some over overlap between the segments (``overlap``). 

673 

674 The function returns a list of tuples. Each tuple provides the lower 

675 (``energy_limit_left``) and upper (``energy_limit_right``) bound of one 

676 bin, which are then to be used to set ``energy_limit_left`` and 

677 ``energy_limit_right`` when initializing a :class:`WangLandauEnsemble` 

678 instance. 

679 

680 N.B.: The left-most/right-most bin has no lower/upper bound (set to 

681 ``None``). 

682 

683 Parameters 

684 ---------- 

685 n_bins 

686 number of bins 

687 energy_spacing 

688 defines the bin size of the energy grid used by the Wang-Landau 

689 simulation, see :class:`WangLandauEnsemble` for details 

690 minimum_energy 

691 an estimate for the lowest energy to be encountered in this system 

692 maximum_energy 

693 an estimate for the highest energy to be encountered in this system 

694 overlap 

695 amount of overlap between bins in units of ``energy_spacing`` 

696 bin_size_exponent 

697 *Expert option*: This parameter allows one to generate a non-uniform 

698 distribution of bin sizes. If ``bin_size_exponent`` is smaller than 

699 one bins at the lower and upper end of the energy range (specified via 

700 ``minimum_energy`` and ``maximum_energy``) will be shrunk relative to 

701 the bins in the middle of the energy range. In principle this can be 

702 used one to achieve a more even distribution of computational load 

703 between the individual Wang-Landau simulations. 

704 """ 

705 

706 limits = np.linspace(-1, 1, n_bins + 1) 

707 limits = np.sign(limits) * np.abs(limits) ** bin_size_exponent 

708 limits *= 0.5 * (maximum_energy - minimum_energy) 

709 limits += 0.5 * (maximum_energy + minimum_energy) 

710 limits[0], limits[-1] = None, None 

711 

712 bounds = [] 

713 for k, (energy_limit_left, energy_limit_right) in enumerate(zip(limits[:-1], limits[1:])): 

714 if energy_limit_left is not None and energy_limit_right is not None and \ 

715 (energy_limit_right - energy_limit_left) / energy_spacing < 2 * overlap: 

716 raise ValueError('Energy window too small. min/max: {}/{}' 

717 .format(energy_limit_right, energy_limit_left) + 

718 ' Try decreasing n_bins ({}) and/or overlap ({}).' 

719 .format(n_bins, overlap)) 

720 if energy_limit_left is not None: 

721 energy_limit_left -= overlap * energy_spacing 

722 if energy_limit_right is not None: 

723 energy_limit_right += overlap * energy_spacing 

724 bounds.append((energy_limit_left, energy_limit_right)) 

725 

726 return bounds