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 typing import List, Union 

2 

3import numpy as np 

4 

5from _icet import _ClusterExpansionCalculator 

6from icet.input_output.logging_tools import logger 

7from ase import Atoms 

8from icet import ClusterExpansion 

9from icet.core.structure import Structure 

10from icet.core.sublattices import Sublattices 

11from mchammer.calculators.base_calculator import BaseCalculator 

12 

13 

14class ClusterExpansionCalculator(BaseCalculator): 

15 """A ClusterExpansionCalculator object enables the efficient 

16 calculation of properties described by a cluster expansion. It is 

17 specific for a particular (supercell) structure and commonly 

18 employed when setting up a Monte Carlo simulation, see 

19 :ref:`ensembles`. 

20 

21 Cluster expansions, e.g., of the energy, typically yield property 

22 values *per site*. When running a Monte Carlo simulation one, 

23 however, considers changes in the *total* energy of the 

24 system. The default behavior is therefore to multiply the output 

25 of the cluster expansion by the number of sites. This behavior can 

26 be changed via the ``scaling`` keyword parameter. 

27 

28 Parameters 

29 ---------- 

30 structure : ase.Atoms 

31 structure for which to set up the calculator 

32 cluster_expansion : ClusterExpansion 

33 cluster expansion from which to build calculator 

34 name 

35 human-readable identifier for this calculator 

36 scaling 

37 scaling factor applied to the property value predicted by the 

38 cluster expansion 

39 use_local_energy_calculator 

40 evaluate energy changes using only the local environment; this method 

41 is generally *much* faster; unless you know what you are doing do *not* 

42 set this option to `False` 

43 """ 

44 

45 def __init__(self, 

46 structure: Atoms, cluster_expansion: ClusterExpansion, 

47 name: str = 'Cluster Expansion Calculator', 

48 scaling: Union[float, int] = None, 

49 use_local_energy_calculator: bool = True) -> None: 

50 super().__init__(name=name) 

51 

52 structure_cpy = structure.copy() 

53 cluster_expansion.prune() 

54 

55 if cluster_expansion._cluster_space.is_supercell_self_interacting(structure): 

56 logger.warning('The ClusterExpansionCalculator self-interacts, ' 

57 'which may lead to erroneous results. To avoid ' 

58 'self-interaction, use a larger supercell or a ' 

59 'cluster space with shorter cutoffs.') 

60 

61 self.use_local_energy_calculator = use_local_energy_calculator 

62 if self.use_local_energy_calculator: 62 ↛ 68line 62 didn't jump to line 68, because the condition on line 62 was never false

63 self.cpp_calc = _ClusterExpansionCalculator( 

64 cluster_space=cluster_expansion.get_cluster_space_copy(), 

65 structure=Structure.from_atoms(structure_cpy), 

66 fractional_position_tolerance=cluster_expansion.fractional_position_tolerance) 

67 

68 self._cluster_expansion = cluster_expansion 

69 if scaling is None: 69 ↛ 72line 69 didn't jump to line 72, because the condition on line 69 was never false

70 self._property_scaling = len(structure) 

71 else: 

72 self._property_scaling = scaling 

73 

74 self._sublattices = self.cluster_expansion._cluster_space.get_sublattices(structure) 

75 

76 @property 

77 def cluster_expansion(self) -> ClusterExpansion: 

78 """ cluster expansion from which calculator was constructed """ 

79 return self._cluster_expansion 

80 

81 def calculate_total(self, *, occupations: List[int]) -> float: 

82 """ 

83 Calculates and returns the total property value of the current 

84 configuration. 

85 

86 Parameters 

87 ---------- 

88 occupations 

89 the entire occupation vector (i.e. list of atomic species) 

90 """ 

91 

92 cv = self.cpp_calc.get_full_cluster_vector(occupations) 

93 return np.dot(cv, self.cluster_expansion.parameters) * self._property_scaling 

94 

95 def calculate_change(self, *, sites: List[int], 

96 current_occupations: List[int], 

97 new_site_occupations: List[int]) -> float: 

98 """ 

99 Calculates and returns the sum of the contributions to the property 

100 due to the sites specified in `local_indices` 

101 

102 Parameters 

103 ---------- 

104 sites 

105 index of sites at which occupations will be changed 

106 current_occupations 

107 entire occupation vector (atomic numbers) before change 

108 new_site_occupations 

109 atomic numbers after change at the sites defined by `sites` 

110 """ 

111 occupations = np.array(current_occupations) 

112 new_site_occupations = np.array(new_site_occupations) 

113 

114 if not self.use_local_energy_calculator: 

115 e_before = self.calculate_total(occupations=occupations) 

116 occupations[sites] = np.array(new_site_occupations) 

117 e_after = self.calculate_total(occupations=occupations) 

118 return e_after - e_before 

119 

120 local_contribution = 0 

121 try: 

122 exclude_indices = [] # type: List[int] 

123 for index in sites: 

124 local_contribution -= self._calculate_local_contribution( 

125 occupations=occupations, index=index, 

126 exclude_indices=exclude_indices) 

127 exclude_indices.append(index) 

128 

129 occupations[sites] = np.array(new_site_occupations) 

130 exclude_indices = [] # type: List[int] 

131 for index in sites: 

132 local_contribution += self._calculate_local_contribution( 

133 occupations=occupations, index=index, 

134 exclude_indices=exclude_indices) 

135 exclude_indices.append(index) 

136 except Exception as e: 

137 msg = 'Caught exception {}. Try setting parameter '.format(e) 

138 msg += 'use_local_energy_calculator to False in init' 

139 raise RuntimeError(msg) 

140 

141 return local_contribution * self._property_scaling 

142 

143 def _calculate_local_contribution(self, occupations: List[int], index: int, 

144 exclude_indices: List[int] = []): 

145 """ 

146 Internal method to calculate the local contribution for one 

147 index. 

148 

149 Parameters 

150 ---------- 

151 occupations 

152 entire occupation vector 

153 index : int 

154 lattice index 

155 exclude_indices 

156 previously calculated indices, these indices will 

157 be ignored in order to avoid double counting bonds 

158 

159 """ 

160 local_cv = self.cpp_calc.get_local_cluster_vector( 

161 occupations, index, exclude_indices) 

162 return np.dot(local_cv, self.cluster_expansion.parameters) 

163 

164 @property 

165 def sublattices(self) -> Sublattices: 

166 """Sublattices of the calculators structure.""" 

167 return self._sublattices