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

1# fmt: off 

2 

3"""Implements the dimensionality scoring parameter. 

4 

5Method is described in: 

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

7components 

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

9Phys. Rev. Materials 3 034003, 2019 

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

11""" 

12from collections import namedtuple 

13 

14import numpy as np 

15 

16from ase.geometry.dimensionality import rank_determination, topology_scaling 

17from ase.geometry.dimensionality.bond_generator import next_bond 

18 

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

20 

21 

22def f(x): 

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

24 return 1 

25 k = 1 / 0.15**2 

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

27 

28 

29def calculate_score(a, b): 

30 return f(b) - f(a) 

31 

32 

33def reduced_histogram(h): 

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

35 return tuple(h) 

36 

37 

38def build_dimtype(h): 

39 h = reduced_histogram(h) 

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

41 

42 

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

44 if score is None: 

45 score = calculate_score(a, b) 

46 

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

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

49 

50 

51def merge_intervals(intervals): 

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

53 

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

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

56 

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

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

59 interval with the highest score. 

60 

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

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

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

64 that the scores sum to 1. 

65 """ 

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

67 

68 merged_intervals = [] 

69 for dimtype in dimtypes: 

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

71 combined_score = sum(e.score for e in relevant) 

72 amin = min(e.a for e in relevant) 

73 bmax = max(e.b for e in relevant) 

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

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

76 best.cdim, score=combined_score) 

77 merged_intervals.append(merged) 

78 return merged_intervals 

79 

80 

81def build_kintervals(atoms, method_name): 

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

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

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

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

86 

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

88 num_atoms = len(atoms) 

89 intervals = [] 

90 kprev = 0 

91 calc = method(num_atoms) 

92 hprev = calc.check() 

93 components_prev, cdim_prev = calc.get_components() 

94 

95 """ 

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

97 the periodic boundary conditions: 

98 """ 

99 end_state = np.zeros(4) 

100 end_dim = sum(atoms.pbc) 

101 end_state[end_dim] = 1 

102 end_state = tuple(end_state) 

103 

104 # Insert each new bond into the component graph. 

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

106 calc.insert_bond(i, j, offset) 

107 h = calc.check() 

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

109 continue 

110 

111 components, cdim = calc.get_components() 

112 

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

114 if k != kprev: 

115 # Only keep intervals of non-zero width 

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

117 components_prev, cdim_prev)) 

118 kprev = k 

119 hprev = h 

120 components_prev = components 

121 cdim_prev = cdim 

122 

123 # Stop once all components are merged 

124 if h == end_state: 

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

126 components, cdim)) 

127 return intervals 

128 

129 

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

131 """Performs a k-interval analysis. 

132 

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

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

135 to low. 

136 

137 Parameters: 

138 

139 atoms: ASE atoms object 

140 The system to analyze. The periodic boundary conditions determine 

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

142 a maximum dimensionality of 2. 

143 method: string 

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

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

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

147 merge: boolean 

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

149 merged. Default: true 

150 

151 Returns: 

152 

153 intervals: list 

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

155 namedtuple with the following field names: 

156 

157 score: float 

158 Dimensionality score in the range [0, 1] 

159 a: float 

160 The start of the k-interval 

161 b: float 

162 The end of the k-interval 

163 dimtype: str 

164 The dimensionality type 

165 h: tuple 

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

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

168 components: array 

169 The component ID of each atom. 

170 cdim: dict 

171 The component dimensionalities 

172 """ 

173 intervals = build_kintervals(atoms, method) 

174 if merge: 

175 intervals = merge_intervals(intervals) 

176 

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

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