Coverage for mchammer/ensembles/canonical_annealing.py: 93%
69 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-09-14 04:08 +0000
1"""Definition of the canonical annealing class."""
3import numpy as np
5from ase import Atoms
6from ase.units import kB
8from .. import DataContainer
9from ..calculators import ClusterExpansionCalculator
10from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble
13class CanonicalAnnealing(ThermodynamicBaseEnsemble):
14 """Instances of this class allow one to carry out simulated annealing
15 in the canonical ensemble, i.e. the temperature is varied in
16 pre-defined fashion while the composition is kept fixed. See
17 :class:`CanonicalEnsemble <mchammer.ensembles.CanonicalEnsemble>`
18 for more information about the standard canonical ensemble.
20 The canonical annealing ensemble can be useful, for example, for
21 finding ground states or generating low energy configurations.
23 The temperature control scheme is selected via the
24 :attr:`cooling_function` keyword argument, while the initial and final
25 temperature are set via the :attr:`T_start` and :attr:`T_stop` arguments.
26 Several pre-defined temperature control schemes are available
27 including ``'linear'`` and ``'exponential'``. In the latter case the
28 temperature varies logarithmatically as a function of the MC step,
29 emulating the exponential temperature dependence of the atomic
30 exchange rate encountered in many materials. It is also possible
31 to provide a user defined cooling function via the keyword
32 argument. This function must comply with the following function
33 header::
35 def cooling_function(step, T_start, T_stop, n_steps):
36 T = ... # compute temperature
37 return T
39 Here :attr:`step` refers to the current MC trial step.
41 Parameters
42 ----------
43 structure
44 Atomic configuration to be used in the Monte Carlo simulation;
45 also defines the initial occupation vector.
46 calculator
47 Calculator to be used for calculating the potential changes
48 that enter the evaluation of the Metropolis criterion.
49 T_start
50 Temperature from which the annealing is started.
51 T_stop
52 Final temperature for annealing.
53 n_steps
54 Number of steps to take in the annealing simulation.
55 cooling_function
56 to use the predefined cooling functions provide a string
57 ``'linear'`` or ``'exponential'``, otherwise provide a function.
58 boltzmann_constant
59 Boltzmann constant :math:`k_B` in appropriate
60 units, i.e., units that are consistent
61 with the underlying cluster expansion
62 and the temperature units. Default: eV/K.
63 user_tag
64 Human-readable tag for ensemble. Default: ``None``.
65 random_seed
66 Seed for the random number generator used in the Monte Carlo simulation.
67 dc_filename
68 Name of file the data container associated with the ensemble
69 will be written to. If the file exists it will be read, the
70 data container will be appended, and the file will be
71 updated/overwritten.
72 data_container_write_period
73 Period in units of seconds at which the data container is
74 written to file. Writing periodically to file provides both
75 a way to examine the progress of the simulation and to back up
76 the data. Default: 600 s.
77 ensemble_data_write_interval
78 interval at which data is written to the data container. This
79 includes for example the current value of the calculator
80 (i.e., usually the energy) as well as ensembles specific fields
81 such as temperature or the number of atoms of different species.
82 Default: Number of sites in the :attr:`structure`.
83 trajectory_write_interval
84 interval at which the current occupation vector of the atomic
85 configuration is written to the data container.
86 Default: Number of sites in the :attr:`structure`.
87 sublattice_probabilities
88 Probability for picking a sublattice when doing a random swap.
89 This should be as long as the number of sublattices and should
90 sum up to 1.
92 """
94 def __init__(self,
95 structure: Atoms,
96 calculator: ClusterExpansionCalculator,
97 T_start: float,
98 T_stop: float,
99 n_steps: int,
100 cooling_function: str = 'exponential',
101 user_tag: str = None,
102 boltzmann_constant: float = kB,
103 random_seed: int = None,
104 dc_filename: str = None,
105 data_container_write_period: float = 600,
106 ensemble_data_write_interval: int = None,
107 trajectory_write_interval: int = None,
108 sublattice_probabilities: list[float] = None) -> None:
110 self._ensemble_parameters = dict(n_steps=n_steps)
112 # add species count to ensemble parameters
113 for sl in calculator.sublattices:
114 for symbol in sl.chemical_symbols:
115 key = 'n_atoms_{}'.format(symbol)
116 count = structure.get_chemical_symbols().count(symbol)
117 self._ensemble_parameters[key] = count
119 super().__init__(
120 structure=structure, calculator=calculator, user_tag=user_tag,
121 random_seed=random_seed,
122 dc_filename=dc_filename,
123 data_container_class=DataContainer,
124 data_container_write_period=data_container_write_period,
125 ensemble_data_write_interval=ensemble_data_write_interval,
126 trajectory_write_interval=trajectory_write_interval,
127 boltzmann_constant=boltzmann_constant)
129 self._temperature = T_start
130 self._T_start = T_start
131 self._T_stop = T_stop
132 self._n_steps = n_steps
134 self._ground_state_candidate = self.configuration.structure
135 self._ground_state_candidate_potential = calculator.calculate_total(
136 occupations=self.configuration.occupations)
138 # setup cooling function
139 if isinstance(cooling_function, str):
140 available = sorted(available_cooling_functions.keys())
141 if cooling_function not in available: 141 ↛ 142line 141 didn't jump to line 142 because the condition on line 141 was never true
142 raise ValueError(
143 'Select from the available cooling_functions {}'.format(available))
144 self._cooling_function = available_cooling_functions[cooling_function]
145 elif callable(cooling_function): 145 ↛ 148line 145 didn't jump to line 148 because the condition on line 145 was always true
146 self._cooling_function = cooling_function
147 else:
148 raise TypeError('cooling_function must be either str or a function')
150 if sublattice_probabilities is None: 150 ↛ 153line 150 didn't jump to line 153 because the condition on line 150 was always true
151 self._swap_sublattice_probabilities = self._get_swap_sublattice_probabilities()
152 else:
153 self._swap_sublattice_probabilities = sublattice_probabilities
155 @property
156 def temperature(self) -> float:
157 """ Current temperature. """
158 return self._temperature
160 @property
161 def T_start(self) -> float:
162 """ Starting temperature. """
163 return self._T_start
165 @property
166 def T_stop(self) -> float:
167 """ Starting temperature """
168 return self._T_stop
170 @property
171 def n_steps(self) -> int:
172 """ Number of steps to carry out. """
173 return self._n_steps
175 @property
176 def estimated_ground_state(self):
177 """ Structure with lowest observed potential during run. """
178 return self._ground_state_candidate.copy()
180 @property
181 def estimated_ground_state_potential(self):
182 """ Lowest observed potential during run. """
183 return self._ground_state_candidate_potential
185 def run(self):
186 """ Runs the annealing simulation. """
187 if self.step >= self.n_steps:
188 raise Exception('Annealing has already finished')
189 super().run(self.n_steps - self.step)
191 def _do_trial_step(self) -> int:
192 """ Carries out one Monte Carlo trial step. """
193 self._temperature = self._cooling_function(
194 self.step, self.T_start, self.T_stop, self.n_steps)
195 sublattice_index = self.get_random_sublattice_index(self._swap_sublattice_probabilities)
196 return self.do_canonical_swap(sublattice_index=sublattice_index)
198 def _get_ensemble_data(self) -> dict:
199 """ Returns the data associated with the ensemble. For the
200 CanonicalAnnealing this specifically includes the temperature.
201 """
202 data = super()._get_ensemble_data()
203 data['temperature'] = self.temperature
204 if data['potential'] < self._ground_state_candidate_potential:
205 self._ground_state_candidate_potential = data['potential']
206 self._ground_state_candidate = self.configuration.structure
207 return data
210def _cooling_linear(step, T_start, T_stop, n_steps):
211 return T_start + (T_stop-T_start) * (step + 1) / n_steps
214def _cooling_exponential(step, T_start, T_stop, n_steps):
215 return T_start - (T_start - T_stop) * np.log(step+1) / np.log(n_steps)
218available_cooling_functions = dict(linear=_cooling_linear, exponential=_cooling_exponential)