Coverage for mchammer/ensembles/sgc_annealing.py: 92%
62 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1"""Definition of the canonical annealing class."""
3from ase import Atoms
4from ase.data import chemical_symbols
5from ase.units import kB
6from typing import Dict, List, Any
8from .. import DataContainer
9from ..calculators import ClusterExpansionCalculator
10from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble
11from .semi_grand_canonical_ensemble import get_chemical_potentials
12from .canonical_annealing import available_cooling_functions
15class SGCAnnealing(ThermodynamicBaseEnsemble):
16 r"""Instances of this class allow one to carry out simulated annealing
17 in the semi grand canonical ensemble, i.e., the temperature is varied in
18 pre-defined fashion while the chemical potential is kept fixed. See
19 :class:`mchammer.ensembles.SemiGrandCanonicalEnsemble` for more information
20 about the ensemble.
22 The temperature control scheme is selected via the
23 :attr:`cooling_function` keyword argument, while the initial and final
24 temperature are set via the :attr:`T_start` and :attr:`T_stop` arguments.
25 Several pre-defined temperature control schemes are available
26 including ``'linear'`` and ``'exponential'``. In the latter case the
27 temperature varies logarithmatically as a function of the MC step,
28 emulating the exponential temperature dependence of the atomic
29 exchange rate encountered in many materials. It is also possible
30 to provide a user defined cooling function via the keyword
31 argument. This function must comply with the following function
32 header::
34 def cooling_function(step, T_start, T_stop, n_steps):
35 T = ... # compute temperature
36 return T
38 Here :attr:`step` refers to the current MC trial step.
40 Parameters
41 ----------
42 structure
43 Atomic configuration to be used in the Monte Carlo simulation;
44 also defines the initial occupation vector.
45 chemical_potentials
46 Chemical potential for each species :math:`\mu_i`. The key
47 denotes the species, the value specifies the chemical potential in
48 units that are consistent with the underlying cluster expansion.
49 calculator
50 Calculator to be used for calculating the potential changes
51 that enter the evaluation of the Metropolis criterion.
52 T_start
53 Temperature from which the annealing is started.
54 T_stop
55 Final temperature for annealing.
56 n_steps
57 Number of steps to take in the annealing simulation.
58 cooling_function
59 to use the predefined cooling functions provide a string
60 ``'linear'`` or ``'exponential'``, otherwise provide a function.
61 boltzmann_constant
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
67 Human-readable tag for ensemble. Default: ``None``.
68 random_seed
69 Seed for the random number generator used in the Monte Carlo simulation.
70 dc_filename
71 Name of file the data container associated with the ensemble
72 will be written to. If the file exists it will be read, the
73 data container will be appended, and the file will be
74 updated/overwritten.
75 data_container_write_period
76 period in units of seconds at which the data container is
77 written to file. Writing periodically to file provides both
78 a way to examine the progress of the simulation and to back up
79 the data. Default: 600 s.
80 ensemble_data_write_interval
81 Interval at which data is written to the data container. This
82 includes for example the current value of the calculator
83 (i.e., usually the energy) as well as ensembles specific fields
84 such as temperature or the number of atoms of different species.
85 Default: Number of sites in the :attr:`structure`.
86 trajectory_write_interval
87 Interval at which the current occupation vector of the atomic
88 configuration is written to the data container.
89 Default: Number of sites in the :attr:`structure`.
90 sublattice_probabilities
91 Probability for picking a sublattice when doing a random swap.
92 This should be as long as the number of sublattices and should
93 sum up to 1.
95 """
97 def __init__(self,
98 structure: Atoms,
99 calculator: ClusterExpansionCalculator,
100 T_start: float,
101 T_stop: float,
102 n_steps: int,
103 chemical_potentials: Dict[str, float],
104 cooling_function: str = 'exponential',
105 boltzmann_constant: float = kB,
106 user_tag: str = None,
107 random_seed: int = None,
108 dc_filename: str = None,
109 data_container: str = None,
110 data_container_write_period: float = 600,
111 ensemble_data_write_interval: int = None,
112 trajectory_write_interval: int = None,
113 sublattice_probabilities: List[float] = None) -> None:
115 self._ensemble_parameters: Dict[str, Any] = dict(n_steps=n_steps)
117 self._chemical_potentials = get_chemical_potentials(chemical_potentials)
119 # add chemical potentials to ensemble parameters
120 self._chemical_potentials = get_chemical_potentials(chemical_potentials)
121 for atnum, chempot in self.chemical_potentials.items():
122 mu_sym = 'mu_{}'.format(chemical_symbols[atnum])
123 self._ensemble_parameters[mu_sym] = chempot
125 super().__init__(
126 structure=structure,
127 calculator=calculator,
128 user_tag=user_tag,
129 random_seed=random_seed,
130 dc_filename=dc_filename,
131 data_container=data_container,
132 data_container_class=DataContainer,
133 data_container_write_period=data_container_write_period,
134 ensemble_data_write_interval=ensemble_data_write_interval,
135 trajectory_write_interval=trajectory_write_interval,
136 boltzmann_constant=boltzmann_constant)
138 self._temperature = T_start
139 self._T_start = T_start
140 self._T_stop = T_stop
141 self._n_steps = n_steps
143 self._ground_state_candidate = self.configuration.structure
144 self._ground_state_candidate_potential = calculator.calculate_total(
145 occupations=self.configuration.occupations)
147 # setup cooling function
148 if isinstance(cooling_function, str):
149 available = sorted(available_cooling_functions.keys())
150 if cooling_function not in available: 150 ↛ 151line 150 didn't jump to line 151, because the condition on line 150 was never true
151 raise ValueError(
152 'Select from the available cooling_functions {}'.format(available))
153 self._cooling_function = available_cooling_functions[cooling_function]
154 elif callable(cooling_function): 154 ↛ 157line 154 didn't jump to line 157, because the condition on line 154 was never false
155 self._cooling_function = cooling_function
156 else:
157 raise TypeError('cooling_function must be either str or a function')
159 if sublattice_probabilities is None: 159 ↛ 162line 159 didn't jump to line 162, because the condition on line 159 was never false
160 self._flip_sublattice_probabilities = self._get_flip_sublattice_probabilities()
161 else:
162 self._flip_sublattice_probabilities = sublattice_probabilities
164 @property
165 def chemical_potentials(self) -> Dict[int, float]:
166 r"""
167 Chemical potentials :math:`\mu_i` (see parameters section above).
168 """
169 return self._chemical_potentials
171 @property
172 def temperature(self) -> float:
173 """ Current temperature. """
174 return self._temperature
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. """
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._flip_sublattice_probabilities)
197 return self.do_sgc_flip(
198 sublattice_index=sublattice_index, chemical_potentials=self.chemical_potentials)
200 def _get_ensemble_data(self) -> Dict:
201 """Returns the data associated with the ensemble. For
202 CanonicalAnnealing this specifically includes the temperature.
203 """
204 data = super()._get_ensemble_data()
206 # species counts
207 data.update(self._get_species_counts())
209 data['temperature'] = self.temperature
210 if data['potential'] < self._ground_state_candidate_potential:
211 self._ground_state_candidate_potential = data['potential']
212 self._ground_state_candidate = self.configuration.structure
213 return data