Coverage for mchammer/ensembles/target_cluster_vector_annealing.py: 98%
93 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-09-16 04:14 +0000
1from ase import Atoms
2from mchammer.calculators.target_vector_calculator import TargetVectorCalculator
3from .canonical_ensemble import CanonicalEnsemble
4from .canonical_annealing import _cooling_exponential
5import numpy as np
6import random
7from icet.input_output.logging_tools import logger
8logger = logger.getChild('target_cluster_vector_annealing')
11class TargetClusterVectorAnnealing:
12 """
13 Instances of this class allow one to carry out simulated annealing
14 towards a target cluster vector. Because it is impossible
15 to know *a priori* which supercell shape accomodates the best
16 match, this ensemble allows the annealing to be done for multiple
17 :class:`Atoms <ase.Atoms>` objects at the same time.
19 Parameters
20 ----------
21 structure
22 Atomic configurations to be used in the Monte Carlo simulation;
23 also defines the initial occupation vectors.
24 calculators
25 Calculators corresponding to each :class:`Atoms <ase.Atoms>` object.
26 T_start
27 Artificial temperature at which annealing is started.
28 T_stop
29 Artificial temperature at which annealing is stopped.
30 random_seed
31 Seed for random number generator used in the Monte Carlo simulation.
32 """
34 def __init__(self, structure: list[Atoms],
35 calculators: list[TargetVectorCalculator],
36 T_start: float = 5.0, T_stop: float = 0.001,
37 random_seed: int = None) -> None:
39 if isinstance(structure, Atoms):
40 raise ValueError(
41 'A list of ASE Atoms (supercells) must be provided')
42 if len(structure) != len(calculators):
43 raise ValueError('There must be as many supercells as there '
44 'are calculators ({} != {})'.format(len(structure),
45 len(calculators)))
47 logger.info('Initializing target cluster vector annealing '
48 'with {} supercells'.format(len(structure)))
50 # random number generator
51 if random_seed is None:
52 self._random_seed = random.randint(0, int(1e16))
53 else:
54 self._random_seed = random_seed
55 random.seed(a=self._random_seed)
57 # Initialize an ensemble for each supercell
58 sub_ensembles = []
59 for ens_id, (supercell, calculator) in enumerate(zip(structure, calculators)):
60 sub_ensembles.append(CanonicalEnsemble(structure=supercell,
61 calculator=calculator,
62 random_seed=random.randint(
63 0, int(1e16)),
64 user_tag='ensemble_{}'.format(
65 ens_id),
66 temperature=T_start,
67 dc_filename=None))
68 self._sub_ensembles = sub_ensembles
69 self._current_score = self._sub_ensembles[0].calculator.calculate_total(
70 self._sub_ensembles[0].configuration.occupations)
71 self._best_score = self._current_score
72 self._best_structure = structure[0].copy()
73 self._temperature = T_start
74 self._T_start = T_start
75 self._T_stop = T_stop
76 self._total_trials = 0
77 self._accepted_trials = 0
78 self._n_steps = 42
80 def generate_structure(self, number_of_trial_steps: int = None) -> Atoms:
81 """
82 Runs a structure annealing simulation.
84 Parameters
85 ----------
86 number_of_trial_steps
87 Total number of trial steps to perform. If ``None``
88 run (on average) 3000 steps per supercell.
89 """
90 if number_of_trial_steps is None:
91 self._n_steps = 3000 * len(self._sub_ensembles)
92 else:
93 self._n_steps = number_of_trial_steps
95 self._temperature = self._T_start
96 self._total_trials = 0
97 self._accepted_trials = 0
98 while self.total_trials < self.n_steps:
99 if self._total_trials % 1000 == 0:
100 logger.info('MC step {}/{} ({} accepted trials, '
101 'temperature {:.3f}), '
102 'best score: {:.3f}'.format(self.total_trials,
103 self.n_steps,
104 self.accepted_trials,
105 self.temperature,
106 self.best_score))
107 self._do_trial_step()
108 return self.best_structure
110 def _do_trial_step(self):
111 """ Carries out one Monte Carlo trial step. """
112 self._temperature = _cooling_exponential(
113 self.total_trials, self.T_start, self.T_stop, self.n_steps)
114 self._total_trials += 1
116 # Choose a supercell
117 ensemble = random.choice(self._sub_ensembles)
119 # Choose two sites and swap
120 sublattice_index = ensemble.get_random_sublattice_index(
121 ensemble._swap_sublattice_probabilities)
122 sites, species = ensemble.configuration.get_swapped_state(
123 sublattice_index)
125 # Update occupations so that the cluster vector (and its score)
126 # can be calculated
127 ensemble.configuration.update_occupations(sites, species)
128 new_score = ensemble.calculator.calculate_total(
129 ensemble.configuration.occupations)
131 if self._acceptance_condition(new_score - self.current_score):
132 self._current_score = new_score
133 self._accepted_trials += 1
135 # Since we are looking for the best structures we want to
136 # keep track of the best one we have found as yet (the
137 # current one may have a worse score)
138 if self._current_score < self._best_score:
139 self._best_structure = ensemble.structure
140 self._best_score = self._current_score
141 else:
142 ensemble.configuration.update_occupations(
143 sites, list(reversed(species)))
145 def _acceptance_condition(self, potential_diff: float) -> bool:
146 """
147 Evaluates Metropolis acceptance criterion.
149 Parameters
150 ----------
151 potential_diff
152 Change in the thermodynamic potential associated with the trial step.
153 """
154 if potential_diff < 0:
155 return True
156 elif abs(self.temperature) < 1e-6: # temperature is numerically zero 156 ↛ 157line 156 didn't jump to line 157 because the condition on line 156 was never true
157 return False
158 else:
159 p = np.exp(-potential_diff / self.temperature)
160 return p > random.random()
162 @property
163 def temperature(self) -> float:
164 """ Current temperature """
165 return self._temperature
167 @property
168 def T_start(self) -> float:
169 """ Starting temperature. """
170 return self._T_start
172 @property
173 def T_stop(self) -> float:
174 """ Stop temperature. """
175 return self._T_stop
177 @property
178 def n_steps(self) -> int:
179 """ Number of steps to carry out. """
180 return self._n_steps
182 @property
183 def total_trials(self) -> int:
184 """ Number of steps carried out so far. """
185 return self._total_trials
187 @property
188 def accepted_trials(self) -> int:
189 """ Number of accepted trials carried out so far. """
190 return self._accepted_trials
192 @property
193 def current_score(self) -> float:
194 """ Current target vector score. """
195 return self._current_score
197 @property
198 def best_score(self) -> float:
199 """ Best target vector score found so far. """
200 return self._best_score
202 @property
203 def best_structure(self) -> float:
204 """ Structure most closely matching target vector so far. """
205 return self._best_structure