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