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 canonical ensemble class."""

3from ase import Atoms

4from ase.units import kB

5from typing import List

7from .. import DataContainer

8from ..calculators.base_calculator import BaseCalculator

9from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble

12class CanonicalEnsemble(ThermodynamicBaseEnsemble):

13 """Instances of this class allow one to simulate systems in the

14 canonical ensemble (:math:N_iVT), i.e. at constant temperature

15 (:math:T), number of atoms of each species (:math:N_i), and

16 volume (:math:V).

18 The probability for a particular state in the canonical ensemble is

19 proportional to the well-known Boltzmann factor,

21 .. math::

23 \\rho_{\\text{C}} \\propto \\exp [ - E / k_B T ].

25 Since the concentrations or equivalently the number of atoms of each

26 species is held fixed in the canonical ensemble, a trial step must

27 conserve the concentrations. This is accomplished by randomly picking two

28 unlike atoms and swapping their identities. The swap is accepted with

29 probability

31 .. math::

33 P = \\min \\{ 1, \\, \\exp [ - \\Delta E / k_B T ] \\},

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

36 swap.

38 The canonical ensemble provides an ideal framework for studying the

39 properties of a system at a specific concentration. Properties such as

40 potential energy or phenomena such as chemical ordering at a specific

41 temperature can conveniently be studied by simulating at that temperature.

42 The canonical ensemble is also a convenient tool for "optimizing" a

43 system, i.e., finding its lowest energy chemical ordering. In practice,

44 this is usually achieved by simulated annealing, i.e. the system is

45 equilibrated at a high temperature, after which the temperature is

46 continuously lowered until the acceptance probability is almost zero. In a

47 well-behaved system, the chemical ordering at that point corresponds to a

48 low-energy structure, possibly the global minimum at that particular

49 concentration.

51 Parameters

52 ----------

53 structure : :class:Atoms <ase.Atoms>

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

55 also defines the initial occupation vector

56 calculator : :class:BaseCalculator <mchammer.calculators.ClusterExpansionCalculator>

57 calculator to be used for calculating the potential changes

58 that enter the evaluation of the Metropolis criterion

59 temperature : float

60 temperature :math:T in appropriate units [commonly Kelvin]

61 boltzmann_constant : float

62 Boltzmann constant :math:k_B in appropriate

63 units, i.e. units that are consistent

64 with the underlying cluster expansion

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

66 user_tag : str

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

68 random_seed : int

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

70 simulation

71 dc_filename : str

72 name of file the data container associated with the ensemble

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

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

75 updated/overwritten

76 data_container_write_period : float

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

78 written to file; writing periodically to file provides both

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

80 the data [default: 600 s]

81 ensemble_data_write_interval : int

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

83 includes for example the current value of the calculator

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

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

86 trajectory_write_interval : int

87 interval at which the current occupation vector of the atomic

88 configuration is written to the data container.

89 sublattice_probabilities : List[float]

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

91 This should be as long as the number of sublattices and should

92 sum up to 1.

95 Example

96 -------

97 The following snippet illustrate how to carry out a simple Monte Carlo

98 simulation in the canonical ensemble. Here, the parameters of the cluster

99 expansion are set to emulate a simple Ising model in order to obtain an

100 example that can be run without modification. In practice, one should of

101 course use a proper cluster expansion::

103 >>> from ase.build import bulk

104 >>> from icet import ClusterExpansion, ClusterSpace

105 >>> from mchammer.calculators import ClusterExpansionCalculator

107 >>> # prepare cluster expansion

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

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

110 >>> # pairs are included)

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

112 >>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])

113 >>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])

115 >>> # prepare initial configuration

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

117 >>> for k in range(5):

118 >>> structure[k].symbol = 'Ag'

120 >>> # set up and run MC simulation

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

122 >>> mc = CanonicalEnsemble(structure=structure, calculator=calc,

123 ... temperature=600,

124 ... dc_filename='myrun_canonical.dc')

125 >>> mc.run(100) # carry out 100 trial swaps

126 """

128 def __init__(self,

129 structure: Atoms,

130 calculator: BaseCalculator,

131 temperature: float,

132 user_tag: str = None,

133 boltzmann_constant: float = kB,

134 random_seed: int = None,

135 dc_filename: str = None,

136 data_container: str = None,

137 data_container_write_period: float = 600,

138 ensemble_data_write_interval: int = None,

139 trajectory_write_interval: int = None,

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

142 self._ensemble_parameters = dict(temperature=temperature)

144 # add species count to ensemble parameters

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

146 for symbol in sub.chemical_symbols])

147 for symbol in symbols:

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

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

150 self._ensemble_parameters[key] = count

152 super().__init__(

153 structure=structure,

154 calculator=calculator,

155 user_tag=user_tag,

156 random_seed=random_seed,

157 data_container=data_container,

158 dc_filename=dc_filename,

159 data_container_class=DataContainer,

160 data_container_write_period=data_container_write_period,

161 ensemble_data_write_interval=ensemble_data_write_interval,

162 trajectory_write_interval=trajectory_write_interval,

163 boltzmann_constant=boltzmann_constant)

165 if sublattice_probabilities is None: 165 ↛ 168line 165 didn't jump to line 168, because the condition on line 165 was never false

166 self._swap_sublattice_probabilities = self._get_swap_sublattice_probabilities()

167 else:

168 self._swap_sublattice_probabilities = sublattice_probabilities

170 @property

171 def temperature(self) -> float:

172 """ Current temperature """

173 return self._ensemble_parameters['temperature']

175 def _do_trial_step(self):

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

177 sublattice_index = self.get_random_sublattice_index(self._swap_sublattice_probabilities)

178 return self.do_canonical_swap(sublattice_index=sublattice_index)