Coverage for mchammer/ensembles/canonical_annealing.py: 93%
70 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1"""Definition of the canonical annealing class."""
3import numpy as np
5from ase import Atoms
6from ase.units import kB
7from typing import Dict, List
9from .. import DataContainer
10from ..calculators import ClusterExpansionCalculator
11from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble
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.
21 The canonical annealing ensemble can be useful, for example, for
22 finding ground states or generating low energy configurations.
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::
36 def cooling_function(step, T_start, T_stop, n_steps):
37 T = ... # compute temperature
38 return T
40 Here :attr:`step` refers to the current MC trial step.
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.
93 """
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:
111 self._ensemble_parameters = dict(n_steps=n_steps)
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
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)
130 self._temperature = T_start
131 self._T_start = T_start
132 self._T_stop = T_stop
133 self._n_steps = n_steps
135 self._ground_state_candidate = self.configuration.structure
136 self._ground_state_candidate_potential = calculator.calculate_total(
137 occupations=self.configuration.occupations)
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')
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
156 @property
157 def temperature(self) -> float:
158 """ Current temperature. """
159 return self._temperature
161 @property
162 def T_start(self) -> float:
163 """ Starting temperature. """
164 return self._T_start
166 @property
167 def T_stop(self) -> float:
168 """ Starting temperature """
169 return self._T_stop
171 @property
172 def n_steps(self) -> int:
173 """ Number of steps to carry out. """
174 return self._n_steps
176 @property
177 def estimated_ground_state(self):
178 """ Structure with lowest observed potential during run. """
179 return self._ground_state_candidate.copy()
181 @property
182 def estimated_ground_state_potential(self):
183 """ Lowest observed potential during run. """
184 return self._ground_state_candidate_potential
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)
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)
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
211def _cooling_linear(step, T_start, T_stop, n_steps):
212 return T_start + (T_stop-T_start) * (step + 1) / n_steps
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)
219available_cooling_functions = dict(linear=_cooling_linear, exponential=_cooling_exponential)