Coverage for mchammer/ensembles/canonical_annealing.py: 93%

70 statements  

« prev     ^ index     » next       coverage.py v7.5.0, created at 2024-12-26 04:12 +0000

1"""Definition of the canonical annealing class.""" 

2 

3import numpy as np 

4 

5from ase import Atoms 

6from ase.units import kB 

7from typing import Dict, List 

8 

9from .. import DataContainer 

10from ..calculators import ClusterExpansionCalculator 

11from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble 

12 

13 

14class CanonicalAnnealing(ThermodynamicBaseEnsemble): 

15 """Instances of this class allow one to carry out simulated annealing 

16 in the canonical ensemble, i.e. the temperature is varied in 

17 pre-defined fashion while the composition is kept fixed. See 

18 :class:`CanonicalEnsemble <mchammer.ensembles.CanonicalEnsemble>` 

19 for more information about the standard canonical ensemble. 

20 

21 The canonical annealing ensemble can be useful, for example, for 

22 finding ground states or generating low energy configurations. 

23 

24 The temperature control scheme is selected via the 

25 :attr:`cooling_function` keyword argument, while the initial and final 

26 temperature are set via the :attr:`T_start` and :attr:`T_stop` arguments. 

27 Several pre-defined temperature control schemes are available 

28 including ``'linear'`` and ``'exponential'``. In the latter case the 

29 temperature varies logarithmatically as a function of the MC step, 

30 emulating the exponential temperature dependence of the atomic 

31 exchange rate encountered in many materials. It is also possible 

32 to provide a user defined cooling function via the keyword 

33 argument. This function must comply with the following function 

34 header:: 

35 

36 def cooling_function(step, T_start, T_stop, n_steps): 

37 T = ... # compute temperature 

38 return T 

39 

40 Here :attr:`step` refers to the current MC trial step. 

41 

42 Parameters 

43 ---------- 

44 structure 

45 Atomic configuration to be used in the Monte Carlo simulation; 

46 also defines the initial occupation vector. 

47 calculator 

48 Calculator to be used for calculating the potential changes 

49 that enter the evaluation of the Metropolis criterion. 

50 T_start 

51 Temperature from which the annealing is started. 

52 T_stop 

53 Final temperature for annealing. 

54 n_steps 

55 Number of steps to take in the annealing simulation. 

56 cooling_function 

57 to use the predefined cooling functions provide a string 

58 ``'linear'`` or ``'exponential'``, otherwise provide a function. 

59 boltzmann_constant 

60 Boltzmann constant :math:`k_B` in appropriate 

61 units, i.e., units that are consistent 

62 with the underlying cluster expansion 

63 and the temperature units. Default: eV/K. 

64 user_tag 

65 Human-readable tag for ensemble. Default: ``None``. 

66 random_seed 

67 Seed for the random number generator used in the Monte Carlo simulation. 

68 dc_filename 

69 Name of file the data container associated with the ensemble 

70 will be written to. If the file exists it will be read, the 

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

72 updated/overwritten. 

73 data_container_write_period 

74 Period in units of seconds at which the data container is 

75 written to file. Writing periodically to file provides both 

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

77 the data. Default: 600 s. 

78 ensemble_data_write_interval 

79 interval at which data is written to the data container. This 

80 includes for example the current value of the calculator 

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

82 such as temperature or the number of atoms of different species. 

83 Default: Number of sites in the :attr:`structure`. 

84 trajectory_write_interval 

85 interval at which the current occupation vector of the atomic 

86 configuration is written to the data container. 

87 Default: Number of sites in the :attr:`structure`. 

88 sublattice_probabilities 

89 Probability for picking a sublattice when doing a random swap. 

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

91 sum up to 1. 

92 

93 """ 

94 

95 def __init__(self, 

96 structure: Atoms, 

97 calculator: ClusterExpansionCalculator, 

98 T_start: float, 

99 T_stop: float, 

100 n_steps: int, 

101 cooling_function: str = 'exponential', 

102 user_tag: str = None, 

103 boltzmann_constant: float = kB, 

104 random_seed: int = None, 

105 dc_filename: str = None, 

106 data_container_write_period: float = 600, 

107 ensemble_data_write_interval: int = None, 

108 trajectory_write_interval: int = None, 

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

110 

111 self._ensemble_parameters = dict(n_steps=n_steps) 

112 

113 # add species count to ensemble parameters 

114 for sl in calculator.sublattices: 

115 for symbol in sl.chemical_symbols: 

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

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

118 self._ensemble_parameters[key] = count 

119 

120 super().__init__( 

121 structure=structure, calculator=calculator, user_tag=user_tag, 

122 random_seed=random_seed, 

123 dc_filename=dc_filename, 

124 data_container_class=DataContainer, 

125 data_container_write_period=data_container_write_period, 

126 ensemble_data_write_interval=ensemble_data_write_interval, 

127 trajectory_write_interval=trajectory_write_interval, 

128 boltzmann_constant=boltzmann_constant) 

129 

130 self._temperature = T_start 

131 self._T_start = T_start 

132 self._T_stop = T_stop 

133 self._n_steps = n_steps 

134 

135 self._ground_state_candidate = self.configuration.structure 

136 self._ground_state_candidate_potential = calculator.calculate_total( 

137 occupations=self.configuration.occupations) 

138 

139 # setup cooling function 

140 if isinstance(cooling_function, str): 

141 available = sorted(available_cooling_functions.keys()) 

142 if cooling_function not in available: 142 ↛ 143line 142 didn't jump to line 143, because the condition on line 142 was never true

143 raise ValueError( 

144 'Select from the available cooling_functions {}'.format(available)) 

145 self._cooling_function = available_cooling_functions[cooling_function] 

146 elif callable(cooling_function): 146 ↛ 149line 146 didn't jump to line 149, because the condition on line 146 was never false

147 self._cooling_function = cooling_function 

148 else: 

149 raise TypeError('cooling_function must be either str or a function') 

150 

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

152 self._swap_sublattice_probabilities = self._get_swap_sublattice_probabilities() 

153 else: 

154 self._swap_sublattice_probabilities = sublattice_probabilities 

155 

156 @property 

157 def temperature(self) -> float: 

158 """ Current temperature. """ 

159 return self._temperature 

160 

161 @property 

162 def T_start(self) -> float: 

163 """ Starting temperature. """ 

164 return self._T_start 

165 

166 @property 

167 def T_stop(self) -> float: 

168 """ Starting temperature """ 

169 return self._T_stop 

170 

171 @property 

172 def n_steps(self) -> int: 

173 """ Number of steps to carry out. """ 

174 return self._n_steps 

175 

176 @property 

177 def estimated_ground_state(self): 

178 """ Structure with lowest observed potential during run. """ 

179 return self._ground_state_candidate.copy() 

180 

181 @property 

182 def estimated_ground_state_potential(self): 

183 """ Lowest observed potential during run. """ 

184 return self._ground_state_candidate_potential 

185 

186 def run(self): 

187 """ Runs the annealing simulation. """ 

188 if self.step >= self.n_steps: 

189 raise Exception('Annealing has already finished') 

190 super().run(self.n_steps - self.step) 

191 

192 def _do_trial_step(self) -> int: 

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

194 self._temperature = self._cooling_function( 

195 self.step, self.T_start, self.T_stop, self.n_steps) 

196 sublattice_index = self.get_random_sublattice_index(self._swap_sublattice_probabilities) 

197 return self.do_canonical_swap(sublattice_index=sublattice_index) 

198 

199 def _get_ensemble_data(self) -> Dict: 

200 """ Returns the data associated with the ensemble. For the 

201 CanonicalAnnealing this specifically includes the temperature. 

202 """ 

203 data = super()._get_ensemble_data() 

204 data['temperature'] = self.temperature 

205 if data['potential'] < self._ground_state_candidate_potential: 

206 self._ground_state_candidate_potential = data['potential'] 

207 self._ground_state_candidate = self.configuration.structure 

208 return data 

209 

210 

211def _cooling_linear(step, T_start, T_stop, n_steps): 

212 return T_start + (T_stop-T_start) * (step + 1) / n_steps 

213 

214 

215def _cooling_exponential(step, T_start, T_stop, n_steps): 

216 return T_start - (T_start - T_stop) * np.log(step+1) / np.log(n_steps) 

217 

218 

219available_cooling_functions = dict(linear=_cooling_linear, exponential=_cooling_exponential)