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

3import random

5from collections import OrderedDict

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

8import numpy as np

10from ase import Atoms

12from .. import WangLandauDataContainer

13from ..calculators.base_calculator import BaseCalculator

14from .thermodynamic_base_ensemble import BaseEnsemble

15from icet.input_output.logging_tools import logger

17logger = logger.getChild('wang_landau_ensemble')

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

30 The algorithm is initialized as follows.

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.

37 It then proceeds as follows.

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

40 #. Accept or reject the new configuration with probability

42 .. math::

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

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

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:

64 .. math::

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

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.

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

161 >>> from ase import Atoms

162 >>> from icet import ClusterExpansion, ClusterSpace

163 >>> from mchammer.calculators import ClusterExpansionCalculator

164 >>> from mchammer.ensembles import WangLandauEnsemble

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

171 >>> # prepare initial configuration

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

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

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

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

184 """

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:

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

217 # set default values that are system dependent

218 if flatness_check_interval is None:

219 flatness_check_interval = len(structure) * 1000

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

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

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

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

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

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

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)

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

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]

307 @property

308 def fill_factor(self) -> float:

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

310 return self._fill_factor

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

319 @property

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

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

322 return self._converged

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

332 @flatness_limit.setter

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

334 self._flatness_limit = new_value

335 self._converged = None

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

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

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

354 @flatness_check_interval.setter

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

356 self._flatness_check_interval = new_value

358 def run(self, number_of_trial_steps: int):

359 """

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

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.

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)

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

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)

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.

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)

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

430 """Evaluates Metropolis acceptance criterion.

432 Parameters

433 ----------

434 potential_diff

435 change in the thermodynamic potential associated

436 with the trial step

437 """

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

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]

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

497 # update histograms and entropy counters

498 self._update_entropy(bin_cur)

500 return accept

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

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

511 # check flatness of histogram

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

513 self.step > 0 and self._reached_energy_window:

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

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)

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

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

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

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

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)

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.

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

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

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.

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

616 Returns

617 -------

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

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

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

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

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

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.

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

681 None).

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

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, limits[-1] = None, None

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

726 return bounds