Coverage for mchammer/ensembles/target_cluster_vector_annealing.py: 98%

93 statements  

« prev     ^ index     » next       coverage.py v7.10.1, created at 2025-09-14 04:08 +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') 

9 

10 

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. 

18 

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 """ 

33 

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: 

38 

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))) 

46 

47 logger.info('Initializing target cluster vector annealing ' 

48 'with {} supercells'.format(len(structure))) 

49 

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) 

56 

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 

79 

80 def generate_structure(self, number_of_trial_steps: int = None) -> Atoms: 

81 """ 

82 Runs a structure annealing simulation. 

83 

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 

94 

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 

109 

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 

115 

116 # Choose a supercell 

117 ensemble = random.choice(self._sub_ensembles) 

118 

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) 

124 

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) 

130 

131 if self._acceptance_condition(new_score - self.current_score): 

132 self._current_score = new_score 

133 self._accepted_trials += 1 

134 

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))) 

144 

145 def _acceptance_condition(self, potential_diff: float) -> bool: 

146 """ 

147 Evaluates Metropolis acceptance criterion. 

148 

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() 

161 

162 @property 

163 def temperature(self) -> float: 

164 """ Current temperature """ 

165 return self._temperature 

166 

167 @property 

168 def T_start(self) -> float: 

169 """ Starting temperature. """ 

170 return self._T_start 

171 

172 @property 

173 def T_stop(self) -> float: 

174 """ Stop temperature. """ 

175 return self._T_stop 

176 

177 @property 

178 def n_steps(self) -> int: 

179 """ Number of steps to carry out. """ 

180 return self._n_steps 

181 

182 @property 

183 def total_trials(self) -> int: 

184 """ Number of steps carried out so far. """ 

185 return self._total_trials 

186 

187 @property 

188 def accepted_trials(self) -> int: 

189 """ Number of accepted trials carried out so far. """ 

190 return self._accepted_trials 

191 

192 @property 

193 def current_score(self) -> float: 

194 """ Current target vector score. """ 

195 return self._current_score 

196 

197 @property 

198 def best_score(self) -> float: 

199 """ Best target vector score found so far. """ 

200 return self._best_score 

201 

202 @property 

203 def best_structure(self) -> float: 

204 """ Structure most closely matching target vector so far. """ 

205 return self._best_structure