Coverage for /builds/debichem-team/python-ase/ase/geometry/dimensionality/interval_analysis.py : 100.00%
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
1"""Implements the dimensionality scoring parameter.
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"""
10import numpy as np
11from collections import namedtuple
12from ase.geometry.dimensionality import rank_determination
13from ase.geometry.dimensionality import topology_scaling
14from ase.geometry.dimensionality.bond_generator import next_bond
17Kinterval = namedtuple('KInterval', 'dimtype score a b h components cdim')
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)
27def calculate_score(a, b):
28 return f(b) - f(a)
31def reduced_histogram(h):
32 h = [int(e > 0) for e in h]
33 return tuple(h)
36def build_dimtype(h):
37 h = reduced_histogram(h)
38 return ''.join([str(i) for i, e in enumerate(h) if e > 0]) + 'D'
41def build_kinterval(a, b, h, components, cdim, score=None):
42 if score is None:
43 score = calculate_score(a, b)
45 return Kinterval(dimtype=build_dimtype(h), score=score,
46 a=a, b=b, h=h, components=components, cdim=cdim)
49def merge_intervals(intervals):
50 """Merges intervals of the same dimensionality type.
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.
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.
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 = set([e.dimtype for e in intervals])
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
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]
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()
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)
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
109 components, cdim = calc.get_components()
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
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
128def analyze_kintervals(atoms, method='RDA', merge=True):
129 """Performs a k-interval analysis.
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.
135 Parameters:
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
149 Returns:
151 intervals: list
152 List of KIntervals for each interval identified. A KInterval is a
153 namedtuple with the following field names:
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)
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))