Coverage for /builds/debichem-team/python-ase/ase/geometry/dimensionality/isolation.py: 99.34%
152 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1"""
2Implements functions for extracting ('isolating') a low-dimensional material
3component in its own unit cell.
5This uses the rank-determination method 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"""
12import collections
13import itertools
15import numpy as np
17from ase import Atoms
18from ase.geometry.cell import complete_cell
19from ase.geometry.dimensionality import (
20 analyze_dimensionality,
21 rank_determination,
22)
23from ase.geometry.dimensionality.bond_generator import next_bond
24from ase.geometry.dimensionality.interval_analysis import merge_intervals
27def orthogonal_basis(X, Y=None):
29 is_1d = Y is None
30 b = np.zeros((3, 3))
31 b[0] = X
32 if not is_1d:
33 b[1] = Y
34 b = complete_cell(b)
36 Q = np.linalg.qr(b.T)[0].T
37 if np.dot(b[0], Q[0]) < 0:
38 Q[0] = -Q[0]
40 if np.dot(b[2], Q[1]) < 0:
41 Q[1] = -Q[1]
43 if np.linalg.det(Q) < 0:
44 Q[2] = -Q[2]
46 if is_1d:
47 Q = Q[[1, 2, 0]]
49 return Q
52def select_cutoff(atoms):
53 intervals = analyze_dimensionality(atoms, method='RDA', merge=False)
54 dimtype = max(merge_intervals(intervals), key=lambda x: x.score).dimtype
55 m = next(e for e in intervals if e.dimtype == dimtype)
56 if m.b == float("inf"):
57 return m.a + 0.1
58 else:
59 return (m.a + m.b) / 2
62def traverse_graph(atoms, kcutoff):
63 if kcutoff is None:
64 kcutoff = select_cutoff(atoms)
66 rda = rank_determination.RDA(len(atoms))
67 for (k, i, j, offset) in next_bond(atoms):
68 if k > kcutoff:
69 break
70 rda.insert_bond(i, j, offset)
72 rda.check()
73 return rda.graph.find_all(), rda.all_visited, rda.ranks
76def build_supercomponent(atoms, components, k, v, anchor=True):
77 # build supercomponent by mapping components into visited cells
78 positions = []
79 numbers = []
81 for c, offset in dict(v[::-1]).items():
82 indices = np.where(components == c)[0]
83 ps = atoms.positions + np.dot(offset, atoms.get_cell())
84 positions.extend(ps[indices])
85 numbers.extend(atoms.numbers[indices])
86 positions = np.array(positions)
87 numbers = np.array(numbers)
89 # select an 'anchor' atom, which will lie at the origin
90 anchor_index = next(i for i in range(len(atoms)) if components[i] == k)
91 if anchor:
92 positions -= atoms.positions[anchor_index]
93 return positions, numbers
96def select_chain_rotation(scaled):
98 best = (-1, [1, 0, 0])
99 for s in scaled:
100 vhat = np.array([s[0], s[1], 0])
101 norm = np.linalg.norm(vhat)
102 if norm < 1E-6:
103 continue
104 vhat /= norm
105 obj = np.sum(np.dot(scaled, vhat)**2)
106 best = max(best, (obj, vhat), key=lambda x: x[0])
107 _, vhat = best
108 cost, sint, _ = vhat
109 rot = np.array([[cost, -sint, 0], [sint, cost, 0], [0, 0, 1]])
110 return np.dot(scaled, rot)
113def isolate_chain(atoms, components, k, v):
115 # identify the vector along the chain; this is the new cell vector
116 basis_points = np.array([offset for c, offset in v if c == k])
117 assert len(basis_points) >= 2
118 assert (0, 0, 0) in [tuple(e) for e in basis_points]
120 sizes = np.linalg.norm(basis_points, axis=1)
121 index = np.argsort(sizes)[1]
122 basis = basis_points[index]
123 vector = np.dot(basis, atoms.get_cell())
124 norm = np.linalg.norm(vector)
125 vhat = vector / norm
127 # project atoms into new basis
128 positions, numbers = build_supercomponent(atoms, components, k, v)
129 scaled = np.dot(positions, orthogonal_basis(vhat).T / norm)
131 # move atoms into new cell
132 scaled[:, 2] %= 1.0
134 # subtract barycentre in x and y directions
135 scaled[:, :2] -= np.mean(scaled, axis=0)[:2]
137 # pick a good chain rotation (i.e. non-random)
138 scaled = select_chain_rotation(scaled)
140 # make cell large enough in x and y directions
141 init_cell = norm * np.eye(3)
142 pos = np.dot(scaled, init_cell)
143 rmax = np.max(np.linalg.norm(pos[:, :2], axis=1))
144 rmax = max(1, rmax)
145 cell = np.diag([4 * rmax, 4 * rmax, norm])
147 # construct a new atoms object containing the isolated chain
148 return Atoms(numbers=numbers, positions=pos, cell=cell, pbc=[0, 0, 1])
151def construct_inplane_basis(atoms, k, v):
153 basis_points = np.array([offset for c, offset in v if c == k])
154 assert len(basis_points) >= 3
155 assert (0, 0, 0) in [tuple(e) for e in basis_points]
157 sizes = np.linalg.norm(basis_points, axis=1)
158 indices = np.argsort(sizes)
159 basis_points = basis_points[indices]
161 # identify primitive basis
162 best = (float("inf"), None)
163 for u, v in itertools.combinations(basis_points, 2):
165 basis = np.array([[0, 0, 0], u, v])
166 if np.linalg.matrix_rank(basis) < 2:
167 continue
169 a = np.dot(u, atoms.get_cell())
170 b = np.dot(v, atoms.get_cell())
171 norm = np.linalg.norm(np.cross(a, b))
172 best = min(best, (norm, a, b), key=lambda x: x[0])
173 _, a, b = best
174 return a, b, orthogonal_basis(a, b)
177def isolate_monolayer(atoms, components, k, v):
179 a, b, basis = construct_inplane_basis(atoms, k, v)
181 # project atoms into new basis
182 c = np.cross(a, b)
183 c /= np.linalg.norm(c)
184 init_cell = np.dot(np.array([a, b, c]), basis.T)
186 positions, numbers = build_supercomponent(atoms, components, k, v)
187 scaled = np.linalg.solve(init_cell.T, np.dot(positions, basis.T).T).T
189 # move atoms into new cell
190 scaled[:, :2] %= 1.0
192 # subtract barycentre in z-direction
193 scaled[:, 2] -= np.mean(scaled, axis=0)[2]
195 # make cell large enough in z-direction
196 pos = np.dot(scaled, init_cell)
197 zmax = np.max(np.abs(pos[:, 2]))
198 cell = np.copy(init_cell)
199 cell[2] *= 4 * zmax
201 # construct a new atoms object containing the isolated chain
202 return Atoms(numbers=numbers, positions=pos, cell=cell, pbc=[1, 1, 0])
205def isolate_bulk(atoms, components, k, v):
206 positions, numbers = build_supercomponent(atoms, components, k, v,
207 anchor=False)
208 atoms = Atoms(numbers=numbers, positions=positions, cell=atoms.cell,
209 pbc=[1, 1, 1])
210 atoms.wrap(eps=0)
211 return atoms
214def isolate_cluster(atoms, components, k, v):
215 positions, numbers = build_supercomponent(atoms, components, k, v)
216 positions -= np.min(positions, axis=0)
217 cell = np.diag(np.max(positions, axis=0))
218 return Atoms(numbers=numbers, positions=positions, cell=cell, pbc=False)
221def isolate_components(atoms, kcutoff=None):
222 """Isolates components by dimensionality type.
224 Given a k-value cutoff the components (connected clusters) are
225 identified. For each component an Atoms object is created, which contains
226 that component only. The geometry of the resulting Atoms object depends
227 on the component dimensionality type:
229 0D: The cell is a tight box around the atoms. pbc=[0, 0, 0].
230 The cell has no physical meaning.
232 1D: The chain is aligned along the z-axis. pbc=[0, 0, 1].
233 The x and y cell directions have no physical meaning.
235 2D: The layer is aligned in the x-y plane. pbc=[1, 1, 0].
236 The z cell direction has no physical meaning.
238 3D: The original cell is used. pbc=[1, 1, 1].
240 Parameters:
242 atoms: ASE atoms object
243 The system to analyze.
244 kcutoff: float
245 The k-value cutoff to use. Default=None, in which case the
246 dimensionality scoring parameter is used to select the cutoff.
248 Returns:
250 components: dict
251 key: the component dimenionalities.
252 values: a list of Atoms objects for each dimensionality type.
253 """
254 data = {}
255 components, all_visited, ranks = traverse_graph(atoms, kcutoff)
257 for k, v in all_visited.items():
258 v = sorted(list(v))
260 # identify the components which constitute the component
261 key = tuple(np.unique([c for c, offset in v]))
262 dim = ranks[k]
264 if dim == 0:
265 data[('0D', key)] = isolate_cluster(atoms, components, k, v)
266 elif dim == 1:
267 data[('1D', key)] = isolate_chain(atoms, components, k, v)
268 elif dim == 2:
269 data[('2D', key)] = isolate_monolayer(atoms, components, k, v)
270 elif dim == 3:
271 data[('3D', key)] = isolate_bulk(atoms, components, k, v)
273 result = collections.defaultdict(list)
274 for (dim, _), atoms in data.items():
275 result[dim].append(atoms)
276 return result