Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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

10 

11 

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:`ase.Atoms` objects at the same time. 

19 

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

27 object 

28 T_start 

29 artificial temperature at which annealing is started 

30 T_stop : float 

31 artificial temperature at which annealing is stopped 

32 random_seed : int 

33 seed for random number generator used in the Monte Carlo 

34 simulation 

35 """ 

36 

37 def __init__(self, structure: List[Atoms], 

38 calculators: List[TargetVectorCalculator], 

39 T_start: float = 5.0, T_stop: float = 0.001, 

40 random_seed: int = None) -> None: 

41 

42 if isinstance(structure, Atoms): 

43 raise ValueError( 

44 'A list of ASE Atoms (supercells) must be provided') 

45 if len(structure) != len(calculators): 

46 raise ValueError('There must be as many supercells as there ' 

47 'are calculators ({} != {})'.format(len(structure), 

48 len(calculators))) 

49 

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

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

52 

53 # random number generator 

54 if random_seed is None: 

55 self._random_seed = random.randint(0, int(1e16)) 

56 else: 

57 self._random_seed = random_seed 

58 random.seed(a=self._random_seed) 

59 

60 # Initialize an ensemble for each supercell 

61 sub_ensembles = [] 

62 for ens_id, (supercell, calculator) in enumerate(zip(structure, calculators)): 

63 sub_ensembles.append(CanonicalEnsemble(structure=supercell, 

64 calculator=calculator, 

65 random_seed=random.randint( 

66 0, int(1e16)), 

67 user_tag='ensemble_{}'.format( 

68 ens_id), 

69 temperature=T_start, 

70 dc_filename=None)) 

71 self._sub_ensembles = sub_ensembles 

72 self._current_score = self._sub_ensembles[0].calculator.calculate_total( 

73 self._sub_ensembles[0].configuration.occupations) 

74 self._best_score = self._current_score 

75 self._best_structure = structure[0].copy() 

76 self._temperature = T_start 

77 self._T_start = T_start 

78 self._T_stop = T_stop 

79 self._total_trials = 0 

80 self._accepted_trials = 0 

81 self._n_steps = 42 

82 

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

84 """ 

85 Run a structure annealing simulation. 

86 

87 Parameters 

88 ---------- 

89 number_of_trial_steps 

90 Total number of trial steps to perform. If None, 

91 run (on average) 3000 steps per supercell 

92 """ 

93 if number_of_trial_steps is None: 

94 self._n_steps = 3000 * len(self._sub_ensembles) 

95 else: 

96 self._n_steps = number_of_trial_steps 

97 

98 self._temperature = self._T_start 

99 self._total_trials = 0 

100 self._accepted_trials = 0 

101 while self.total_trials < self.n_steps: 

102 if self._total_trials % 1000 == 0: 

103 logger.info('MC step {}/{} ({} accepted trials, ' 

104 'temperature {:.3f}), ' 

105 'best score: {:.3f}'.format(self.total_trials, 

106 self.n_steps, 

107 self.accepted_trials, 

108 self.temperature, 

109 self.best_score)) 

110 self._do_trial_step() 

111 return self.best_structure 

112 

113 def _do_trial_step(self): 

114 """ Carries out one Monte Carlo trial step. """ 

115 self._temperature = _cooling_exponential( 

116 self.total_trials, self.T_start, self.T_stop, self.n_steps) 

117 self._total_trials += 1 

118 

119 # Choose a supercell 

120 ensemble = random.choice(self._sub_ensembles) 

121 

122 # Choose two sites and swap 

123 sublattice_index = ensemble.get_random_sublattice_index( 

124 ensemble._swap_sublattice_probabilities) 

125 sites, species = ensemble.configuration.get_swapped_state( 

126 sublattice_index) 

127 

128 # Update occupations so that the cluster vector (and its score) 

129 # can be calculated 

130 ensemble.configuration.update_occupations(sites, species) 

131 new_score = ensemble.calculator.calculate_total( 

132 ensemble.configuration.occupations) 

133 

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

135 self._current_score = new_score 

136 self._accepted_trials += 1 

137 

138 # Since we are looking for the best structures we want to 

139 # keep track of the best one we have found as yet (the 

140 # current one may have a worse score) 

141 if self._current_score < self._best_score: 

142 self._best_structure = ensemble.structure 

143 self._best_score = self._current_score 

144 else: 

145 ensemble.configuration.update_occupations( 

146 sites, list(reversed(species))) 

147 

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

149 """ 

150 Evaluates Metropolis acceptance criterion. 

151 

152 Parameters 

153 ---------- 

154 potential_diff 

155 change in the thermodynamic potential associated 

156 with the trial step 

157 """ 

158 if potential_diff < 0: 

159 return True 

160 elif abs(self.temperature) < 1e-6: # temperature is numerically zero 160 ↛ 161line 160 didn't jump to line 161, because the condition on line 160 was never true

161 return False 

162 else: 

163 p = np.exp(-potential_diff / self.temperature) 

164 return p > random.random() 

165 

166 @property 

167 def temperature(self) -> float: 

168 """ Current temperature """ 

169 return self._temperature 

170 

171 @property 

172 def T_start(self) -> float: 

173 """ Starting temperature """ 

174 return self._T_start 

175 

176 @property 

177 def T_stop(self) -> float: 

178 """ Stop temperature """ 

179 return self._T_stop 

180 

181 @property 

182 def n_steps(self) -> int: 

183 """ Number of steps to carry out """ 

184 return self._n_steps 

185 

186 @property 

187 def total_trials(self) -> int: 

188 """ Number of steps carried out so far """ 

189 return self._total_trials 

190 

191 @property 

192 def accepted_trials(self) -> int: 

193 """ Number of accepted trials carried out so far """ 

194 return self._accepted_trials 

195 

196 @property 

197 def current_score(self) -> float: 

198 """ Current target vector score """ 

199 return self._current_score 

200 

201 @property 

202 def best_score(self) -> float: 

203 """ Best target vector score found so far """ 

204 return self._best_score 

205 

206 @property 

207 def best_structure(self) -> float: 

208 """ Structure most closely matching target vector so far """ 

209 return self._best_structure