Coverage for mchammer/calculators/cluster_expansion_calculator.py: 93%

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 :class:`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 :attr:`scaling` keyword parameter. 

27 

28 Parameters 

29 ---------- 

30 structure 

31 Structure for which to set up the calculator. 

32 cluster_expansion 

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 self.cpp_calc = _ClusterExpansionCalculator( 

63 cluster_space=cluster_expansion.get_cluster_space_copy(), 

64 structure=Structure.from_atoms(structure_cpy), 

65 fractional_position_tolerance=cluster_expansion.fractional_position_tolerance) 

66 

67 self._cluster_expansion = cluster_expansion 

68 if scaling is None: 

69 self._property_scaling = len(structure) 

70 else: 

71 self._property_scaling = scaling 

72 

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

74 

75 @property 

76 def cluster_expansion(self) -> ClusterExpansion: 

77 """ Cluster expansion from which calculator was set up. """ 

78 return self._cluster_expansion 

79 

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

81 """ 

82 Calculates and returns the total property value of the current 

83 configuration. 

84 

85 Parameters 

86 ---------- 

87 occupations 

88 The entire occupation vector (i.e., list of atomic species). 

89 """ 

90 

91 cv = self.cpp_calc.get_cluster_vector(occupations) 

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

93 

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

95 current_occupations: List[int], 

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

97 """ 

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

99 due to the sites specified in :attr:`sites`. 

100 

101 Parameters 

102 ---------- 

103 sites 

104 Indices of sites at which occupations will be changed. 

105 current_occupations 

106 Entire occupation vector (atomic numbers) before change. 

107 new_site_occupations 

108 Atomic numbers after change at the sites defined by :attr:`sites`. 

109 """ 

110 occupations = np.array(current_occupations) 

111 

112 if not self.use_local_energy_calculator: 

113 e_before = self.calculate_total(occupations=occupations) 

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

115 e_after = self.calculate_total(occupations=occupations) 

116 return e_after - e_before 

117 

118 change = 0.0 

119 try: 

120 for index, new_occupation in zip(sites, new_site_occupations): 

121 change += self._calculate_partial_change(occupations=occupations, 

122 flip_index=index, 

123 new_occupation=new_occupation) 

124 occupations[index] = new_occupation # Safe because we work with a copy 

125 except Exception as e: 

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

127 msg += 'use_local_energy_calculator to False in init' 

128 raise RuntimeError(msg) 

129 return change * self._property_scaling 

130 

131 def _calculate_partial_change(self, occupations: List[int], flip_index: int, 

132 new_occupation: int): 

133 """ 

134 Internal method to calculate the local contribution for one site index. 

135 

136 Parameters 

137 ---------- 

138 occupations 

139 Entire occupation vector. 

140 flip_index 

141 Lattice index for site where a flipped has occurred. 

142 new_occupation 

143 Atomic number of new occupation at site :attr:`flip_index`. 

144 """ 

145 cv_change = self.cpp_calc.get_cluster_vector_change(occupations=occupations, 

146 flip_index=flip_index, 

147 new_occupation=new_occupation) 

148 

149 return np.dot(cv_change, self.cluster_expansion.parameters) 

150 

151 @property 

152 def sublattices(self) -> Sublattices: 

153 """ Sublattices of the calculator structure. """ 

154 return self._sublattices