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

2 

3from ase import Atoms 

4from ase.units import kB 

5from typing import List 

6 

7from .. import DataContainer 

8from ..calculators.base_calculator import BaseCalculator 

9from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble 

10 

11 

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

17 

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

19 proportional to the well-known Boltzmann factor, 

20 

21 .. math:: 

22 

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

24 

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 

30 

31 .. math:: 

32 

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

34 

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

36 swap. 

37 

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. 

50 

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. 

93 

94 

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

102 

103 >>> from ase.build import bulk 

104 >>> from icet import ClusterExpansion, ClusterSpace 

105 >>> from mchammer.calculators import ClusterExpansionCalculator 

106 

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

114 

115 >>> # prepare initial configuration 

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

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

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

119 

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

127 

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: 

141 

142 self._ensemble_parameters = dict(temperature=temperature) 

143 

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 

151 

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) 

164 

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 

169 

170 @property 

171 def temperature(self) -> float: 

172 """ Current temperature """ 

173 return self._ensemble_parameters['temperature'] 

174 

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)