Coverage for /builds/ase/ase/ase/geometry/dimensionality/isolation.py: 99.34%
152 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""
4Implements functions for extracting ('isolating') a low-dimensional material
5component in its own unit cell.
7This uses the rank-determination method described in:
8Definition of a scoring parameter to identify low-dimensional materials
9components
10P.M. Larsen, M. Pandey, M. Strange, and K. W. Jacobsen
11Phys. Rev. Materials 3 034003, 2019
12https://doi.org/10.1103/PhysRevMaterials.3.034003
13"""
14import collections
15import itertools
17import numpy as np
19from ase import Atoms
20from ase.geometry.cell import complete_cell
21from ase.geometry.dimensionality import (
22 analyze_dimensionality,
23 rank_determination,
24)
25from ase.geometry.dimensionality.bond_generator import next_bond
26from ase.geometry.dimensionality.interval_analysis import merge_intervals
29def orthogonal_basis(X, Y=None):
31 is_1d = Y is None
32 b = np.zeros((3, 3))
33 b[0] = X
34 if not is_1d:
35 b[1] = Y
36 b = complete_cell(b)
38 Q = np.linalg.qr(b.T)[0].T
39 if np.dot(b[0], Q[0]) < 0:
40 Q[0] = -Q[0]
42 if np.dot(b[2], Q[1]) < 0:
43 Q[1] = -Q[1]
45 if np.linalg.det(Q) < 0:
46 Q[2] = -Q[2]
48 if is_1d:
49 Q = Q[[1, 2, 0]]
51 return Q
54def select_cutoff(atoms):
55 intervals = analyze_dimensionality(atoms, method='RDA', merge=False)
56 dimtype = max(merge_intervals(intervals), key=lambda x: x.score).dimtype
57 m = next(e for e in intervals if e.dimtype == dimtype)
58 if m.b == float("inf"):
59 return m.a + 0.1
60 else:
61 return (m.a + m.b) / 2
64def traverse_graph(atoms, kcutoff):
65 if kcutoff is None:
66 kcutoff = select_cutoff(atoms)
68 rda = rank_determination.RDA(len(atoms))
69 for (k, i, j, offset) in next_bond(atoms):
70 if k > kcutoff:
71 break
72 rda.insert_bond(i, j, offset)
74 rda.check()
75 return rda.graph.find_all(), rda.all_visited, rda.ranks
78def build_supercomponent(atoms, components, k, v, anchor=True):
79 # build supercomponent by mapping components into visited cells
80 positions = []
81 numbers = []
83 for c, offset in dict(v[::-1]).items():
84 indices = np.where(components == c)[0]
85 ps = atoms.positions + np.dot(offset, atoms.get_cell())
86 positions.extend(ps[indices])
87 numbers.extend(atoms.numbers[indices])
88 positions = np.array(positions)
89 numbers = np.array(numbers)
91 # select an 'anchor' atom, which will lie at the origin
92 anchor_index = next(i for i in range(len(atoms)) if components[i] == k)
93 if anchor:
94 positions -= atoms.positions[anchor_index]
95 return positions, numbers
98def select_chain_rotation(scaled):
100 best = (-1, [1, 0, 0])
101 for s in scaled:
102 vhat = np.array([s[0], s[1], 0])
103 norm = np.linalg.norm(vhat)
104 if norm < 1E-6:
105 continue
106 vhat /= norm
107 obj = np.sum(np.dot(scaled, vhat)**2)
108 best = max(best, (obj, vhat), key=lambda x: x[0])
109 _, vhat = best
110 cost, sint, _ = vhat
111 rot = np.array([[cost, -sint, 0], [sint, cost, 0], [0, 0, 1]])
112 return np.dot(scaled, rot)
115def isolate_chain(atoms, components, k, v):
117 # identify the vector along the chain; this is the new cell vector
118 basis_points = np.array([offset for c, offset in v if c == k])
119 assert len(basis_points) >= 2
120 assert (0, 0, 0) in [tuple(e) for e in basis_points]
122 sizes = np.linalg.norm(basis_points, axis=1)
123 index = np.argsort(sizes)[1]
124 basis = basis_points[index]
125 vector = np.dot(basis, atoms.get_cell())
126 norm = np.linalg.norm(vector)
127 vhat = vector / norm
129 # project atoms into new basis
130 positions, numbers = build_supercomponent(atoms, components, k, v)
131 scaled = np.dot(positions, orthogonal_basis(vhat).T / norm)
133 # move atoms into new cell
134 scaled[:, 2] %= 1.0
136 # subtract barycentre in x and y directions
137 scaled[:, :2] -= np.mean(scaled, axis=0)[:2]
139 # pick a good chain rotation (i.e. non-random)
140 scaled = select_chain_rotation(scaled)
142 # make cell large enough in x and y directions
143 init_cell = norm * np.eye(3)
144 pos = np.dot(scaled, init_cell)
145 rmax = np.max(np.linalg.norm(pos[:, :2], axis=1))
146 rmax = max(1, rmax)
147 cell = np.diag([4 * rmax, 4 * rmax, norm])
149 # construct a new atoms object containing the isolated chain
150 return Atoms(numbers=numbers, positions=pos, cell=cell, pbc=[0, 0, 1])
153def construct_inplane_basis(atoms, k, v):
155 basis_points = np.array([offset for c, offset in v if c == k])
156 assert len(basis_points) >= 3
157 assert (0, 0, 0) in [tuple(e) for e in basis_points]
159 sizes = np.linalg.norm(basis_points, axis=1)
160 indices = np.argsort(sizes)
161 basis_points = basis_points[indices]
163 # identify primitive basis
164 best = (float("inf"), None)
165 for u, v in itertools.combinations(basis_points, 2):
167 basis = np.array([[0, 0, 0], u, v])
168 if np.linalg.matrix_rank(basis) < 2:
169 continue
171 a = np.dot(u, atoms.get_cell())
172 b = np.dot(v, atoms.get_cell())
173 norm = np.linalg.norm(np.cross(a, b))
174 best = min(best, (norm, a, b), key=lambda x: x[0])
175 _, a, b = best
176 return a, b, orthogonal_basis(a, b)
179def isolate_monolayer(atoms, components, k, v):
181 a, b, basis = construct_inplane_basis(atoms, k, v)
183 # project atoms into new basis
184 c = np.cross(a, b)
185 c /= np.linalg.norm(c)
186 init_cell = np.dot(np.array([a, b, c]), basis.T)
188 positions, numbers = build_supercomponent(atoms, components, k, v)
189 scaled = np.linalg.solve(init_cell.T, np.dot(positions, basis.T).T).T
191 # move atoms into new cell
192 scaled[:, :2] %= 1.0
194 # subtract barycentre in z-direction
195 scaled[:, 2] -= np.mean(scaled, axis=0)[2]
197 # make cell large enough in z-direction
198 pos = np.dot(scaled, init_cell)
199 zmax = np.max(np.abs(pos[:, 2]))
200 cell = np.copy(init_cell)
201 cell[2] *= 4 * zmax
203 # construct a new atoms object containing the isolated chain
204 return Atoms(numbers=numbers, positions=pos, cell=cell, pbc=[1, 1, 0])
207def isolate_bulk(atoms, components, k, v):
208 positions, numbers = build_supercomponent(atoms, components, k, v,
209 anchor=False)
210 atoms = Atoms(numbers=numbers, positions=positions, cell=atoms.cell,
211 pbc=[1, 1, 1])
212 atoms.wrap(eps=0)
213 return atoms
216def isolate_cluster(atoms, components, k, v):
217 positions, numbers = build_supercomponent(atoms, components, k, v)
218 positions -= np.min(positions, axis=0)
219 cell = np.diag(np.max(positions, axis=0))
220 return Atoms(numbers=numbers, positions=positions, cell=cell, pbc=False)
223def isolate_components(atoms, kcutoff=None):
224 """Isolates components by dimensionality type.
226 Given a k-value cutoff the components (connected clusters) are
227 identified. For each component an Atoms object is created, which contains
228 that component only. The geometry of the resulting Atoms object depends
229 on the component dimensionality type:
231 0D: The cell is a tight box around the atoms. pbc=[0, 0, 0].
232 The cell has no physical meaning.
234 1D: The chain is aligned along the z-axis. pbc=[0, 0, 1].
235 The x and y cell directions have no physical meaning.
237 2D: The layer is aligned in the x-y plane. pbc=[1, 1, 0].
238 The z cell direction has no physical meaning.
240 3D: The original cell is used. pbc=[1, 1, 1].
242 Parameters:
244 atoms: ASE atoms object
245 The system to analyze.
246 kcutoff: float
247 The k-value cutoff to use. Default=None, in which case the
248 dimensionality scoring parameter is used to select the cutoff.
250 Returns:
252 components: dict
253 key: the component dimenionalities.
254 values: a list of Atoms objects for each dimensionality type.
255 """
256 data = {}
257 components, all_visited, ranks = traverse_graph(atoms, kcutoff)
259 for k, v in all_visited.items():
260 v = sorted(list(v))
262 # identify the components which constitute the component
263 key = tuple(np.unique([c for c, offset in v]))
264 dim = ranks[k]
266 if dim == 0:
267 data[('0D', key)] = isolate_cluster(atoms, components, k, v)
268 elif dim == 1:
269 data[('1D', key)] = isolate_chain(atoms, components, k, v)
270 elif dim == 2:
271 data[('2D', key)] = isolate_monolayer(atoms, components, k, v)
272 elif dim == 3:
273 data[('3D', key)] = isolate_bulk(atoms, components, k, v)
275 result = collections.defaultdict(list)
276 for (dim, _), atoms in data.items():
277 result[dim].append(atoms)
278 return result