Coverage for /builds/kinetik161/ase/ase/geometry/dimensionality/interval_analysis.py: 100.00%

1"""Implements the dimensionality scoring parameter. 

2 

3Method is described in: 

4Definition of a scoring parameter to identify low-dimensional materials 

5components 

6P.M. Larsen, M. Pandey, M. Strange, and K. W. Jacobsen 

7Phys. Rev. Materials 3 034003, 2019 

8https://doi.org/10.1103/PhysRevMaterials.3.034003 

9""" 

10from collections import namedtuple 

11 

12import numpy as np 

13 

14from ase.geometry.dimensionality import rank_determination, topology_scaling 

15from ase.geometry.dimensionality.bond_generator import next_bond 

16 

17KInterval = namedtuple('KInterval', 'dimtype score a b h components cdim') 

18 

19 

20def f(x): 

21 if x == float("inf"): 

22 return 1 

23 k = 1 / 0.15**2 

24 return k * max(0, x - 1)**2 / (1. + k * max(0, x - 1)**2) 

25 

26 

27def calculate_score(a, b): 

28 return f(b) - f(a) 

29 

30 

31def reduced_histogram(h): 

32 h = [int(e > 0) for e in h] 

33 return tuple(h) 

34 

35 

36def build_dimtype(h): 

37 h = reduced_histogram(h) 

38 return ''.join([str(i) for i, e in enumerate(h) if e > 0]) + 'D' 

39 

40 

41def build_kinterval(a, b, h, components, cdim, score=None): 

42 if score is None: 

43 score = calculate_score(a, b) 

44 

45 return KInterval(dimtype=build_dimtype(h), score=score, 

46 a=a, b=b, h=h, components=components, cdim=cdim) 

47 

48 

49def merge_intervals(intervals): 

50 """Merges intervals of the same dimensionality type. 

51 

52 For example, two histograms with component histograms [10, 4, 0, 0] and 

53 [6, 2, 0, 0] are both 01D structures so they will be merged. 

54 

55 Intervals are merged by summing the scores, and taking the minimum and 

56 maximum k-values. Component IDs in the merged interval are taken from the 

57 interval with the highest score. 

58 

59 On rare occasions, intervals to be merged are not adjacent. In this case, 

60 the score of the merged interval is not equal to the score which would be 

61 calculated from its k-interval. This is necessary to maintain the property 

62 that the scores sum to 1. 

63 """ 

64 dimtypes = {e.dimtype for e in intervals} 

65 

66 merged_intervals = [] 

67 for dimtype in dimtypes: 

68 relevant = [e for e in intervals if e.dimtype == dimtype] 

69 combined_score = sum([e.score for e in relevant]) 

70 amin = min([e.a for e in relevant]) 

71 bmax = max([e.b for e in relevant]) 

72 best = max(relevant, key=lambda x: x.score) 

73 merged = build_kinterval(amin, bmax, best.h, best.components, 

74 best.cdim, score=combined_score) 

75 merged_intervals.append(merged) 

76 return merged_intervals 

77 

78 

79def build_kintervals(atoms, method_name): 

80 """The interval analysis is performed by inserting bonds one at a time 

81 until the component analysis finds a single component.""" 

82 method = {'RDA': rank_determination.RDA, 

83 'TSA': topology_scaling.TSA}[method_name] 

84 

85 assert all(e in [0, 1] for e in atoms.pbc) 

86 num_atoms = len(atoms) 

87 intervals = [] 

88 kprev = 0 

89 calc = method(num_atoms) 

90 hprev = calc.check() 

91 components_prev, cdim_prev = calc.get_components() 

92 

93 """ 

94 The end state is a single component, whose dimensionality depends on 

95 the periodic boundary conditions: 

96 """ 

97 end_state = np.zeros(4) 

98 end_dim = sum(atoms.pbc) 

99 end_state[end_dim] = 1 

100 end_state = tuple(end_state) 

101 

102 # Insert each new bond into the component graph. 

103 for (k, i, j, offset) in next_bond(atoms): 

104 calc.insert_bond(i, j, offset) 

105 h = calc.check() 

106 if h == hprev: # Test if any components were merged 

107 continue 

108 

109 components, cdim = calc.get_components() 

110 

111 # If any components were merged, create a new interval 

112 if k != kprev: 

113 # Only keep intervals of non-zero width 

114 intervals.append(build_kinterval(kprev, k, hprev, 

115 components_prev, cdim_prev)) 

116 kprev = k 

117 hprev = h 

118 components_prev = components 

119 cdim_prev = cdim 

120 

121 # Stop once all components are merged 

122 if h == end_state: 

123 intervals.append(build_kinterval(k, float("inf"), h, 

124 components, cdim)) 

125 return intervals 

126 

127 

128def analyze_kintervals(atoms, method='RDA', merge=True): 

129 """Performs a k-interval analysis. 

130 

131 In each k-interval the components (connected clusters) are identified. 

132 The intervals are sorted according to the scoring parameter, from high 

133 to low. 

134 

135 Parameters: 

136 

137 atoms: ASE atoms object 

138 The system to analyze. The periodic boundary conditions determine 

139 the maximum achievable component dimensionality, i.e. pbc=[1,1,0] sets 

140 a maximum dimensionality of 2. 

141 method: string 

142 Analysis method to use, either 'RDA' (default option) or 'TSA'. 

143 These correspond to the Rank Determination Algorithm of Mounet et al. 

144 and the Topological Scaling Algorithm (TSA) of Ashton et al. 

145 merge: boolean 

146 Decides if k-intervals of the same type (e.g. 01D or 3D) should be 

147 merged. Default: true 

148 

149 Returns: 

150 

151 intervals: list 

152 List of KIntervals for each interval identified. A KInterval is a 

153 namedtuple with the following field names: 

154 

155 score: float 

156 Dimensionality score in the range [0, 1] 

157 a: float 

158 The start of the k-interval 

159 b: float 

160 The end of the k-interval 

161 dimtype: str 

162 The dimensionality type 

163 h: tuple 

164 The histogram of the number of components of each dimensionality. 

165 For example, (8, 0, 3, 0) means eight 0D and three 2D components. 

166 components: array 

167 The component ID of each atom. 

168 cdim: dict 

169 The component dimensionalities 

170 """ 

171 intervals = build_kintervals(atoms, method) 

172 if merge: 

173 intervals = merge_intervals(intervals) 

174 

175 # Sort intervals by score. Interval width resolves ambiguity when score=0. 

176 return sorted(intervals, reverse=True, key=lambda x: (x.score, x.b - x.a))