Coverage for /builds/kinetik161/ase/ase/geometry/dimensionality/isolation.py: 99.34%
152 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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 (analyze_dimensionality,
20 rank_determination)
21from ase.geometry.dimensionality.bond_generator import next_bond
22from ase.geometry.dimensionality.interval_analysis import merge_intervals
25def orthogonal_basis(X, Y=None):
27 is_1d = Y is None
28 b = np.zeros((3, 3))
29 b[0] = X
30 if not is_1d:
31 b[1] = Y
32 b = complete_cell(b)
34 Q = np.linalg.qr(b.T)[0].T
35 if np.dot(b[0], Q[0]) < 0:
36 Q[0] = -Q[0]
38 if np.dot(b[2], Q[1]) < 0:
39 Q[1] = -Q[1]
41 if np.linalg.det(Q) < 0:
42 Q[2] = -Q[2]
44 if is_1d:
45 Q = Q[[1, 2, 0]]
47 return Q
50def select_cutoff(atoms):
51 intervals = analyze_dimensionality(atoms, method='RDA', merge=False)
52 dimtype = max(merge_intervals(intervals), key=lambda x: x.score).dimtype
53 m = next(e for e in intervals if e.dimtype == dimtype)
54 if m.b == float("inf"):
55 return m.a + 0.1
56 else:
57 return (m.a + m.b) / 2
60def traverse_graph(atoms, kcutoff):
61 if kcutoff is None:
62 kcutoff = select_cutoff(atoms)
64 rda = rank_determination.RDA(len(atoms))
65 for (k, i, j, offset) in next_bond(atoms):
66 if k > kcutoff:
67 break
68 rda.insert_bond(i, j, offset)
70 rda.check()
71 return rda.graph.find_all(), rda.all_visited, rda.ranks
74def build_supercomponent(atoms, components, k, v, anchor=True):
75 # build supercomponent by mapping components into visited cells
76 positions = []
77 numbers = []
79 for c, offset in dict(v[::-1]).items():
80 indices = np.where(components == c)[0]
81 ps = atoms.positions + np.dot(offset, atoms.get_cell())
82 positions.extend(ps[indices])
83 numbers.extend(atoms.numbers[indices])
84 positions = np.array(positions)
85 numbers = np.array(numbers)
87 # select an 'anchor' atom, which will lie at the origin
88 anchor_index = next(i for i in range(len(atoms)) if components[i] == k)
89 if anchor:
90 positions -= atoms.positions[anchor_index]
91 return positions, numbers
94def select_chain_rotation(scaled):
96 best = (-1, [1, 0, 0])
97 for s in scaled:
98 vhat = np.array([s[0], s[1], 0])
99 norm = np.linalg.norm(vhat)
100 if norm < 1E-6:
101 continue
102 vhat /= norm
103 obj = np.sum(np.dot(scaled, vhat)**2)
104 best = max(best, (obj, vhat), key=lambda x: x[0])
105 _, vhat = best
106 cost, sint, _ = vhat
107 rot = np.array([[cost, -sint, 0], [sint, cost, 0], [0, 0, 1]])
108 return np.dot(scaled, rot)
111def isolate_chain(atoms, components, k, v):
113 # identify the vector along the chain; this is the new cell vector
114 basis_points = np.array([offset for c, offset in v if c == k])
115 assert len(basis_points) >= 2
116 assert (0, 0, 0) in [tuple(e) for e in basis_points]
118 sizes = np.linalg.norm(basis_points, axis=1)
119 index = np.argsort(sizes)[1]
120 basis = basis_points[index]
121 vector = np.dot(basis, atoms.get_cell())
122 norm = np.linalg.norm(vector)
123 vhat = vector / norm
125 # project atoms into new basis
126 positions, numbers = build_supercomponent(atoms, components, k, v)
127 scaled = np.dot(positions, orthogonal_basis(vhat).T / norm)
129 # move atoms into new cell
130 scaled[:, 2] %= 1.0
132 # subtract barycentre in x and y directions
133 scaled[:, :2] -= np.mean(scaled, axis=0)[:2]
135 # pick a good chain rotation (i.e. non-random)
136 scaled = select_chain_rotation(scaled)
138 # make cell large enough in x and y directions
139 init_cell = norm * np.eye(3)
140 pos = np.dot(scaled, init_cell)
141 rmax = np.max(np.linalg.norm(pos[:, :2], axis=1))
142 rmax = max(1, rmax)
143 cell = np.diag([4 * rmax, 4 * rmax, norm])
145 # construct a new atoms object containing the isolated chain
146 return Atoms(numbers=numbers, positions=pos, cell=cell, pbc=[0, 0, 1])
149def construct_inplane_basis(atoms, k, v):
151 basis_points = np.array([offset for c, offset in v if c == k])
152 assert len(basis_points) >= 3
153 assert (0, 0, 0) in [tuple(e) for e in basis_points]
155 sizes = np.linalg.norm(basis_points, axis=1)
156 indices = np.argsort(sizes)
157 basis_points = basis_points[indices]
159 # identify primitive basis
160 best = (float("inf"), None)
161 for u, v in itertools.combinations(basis_points, 2):
163 basis = np.array([[0, 0, 0], u, v])
164 if np.linalg.matrix_rank(basis) < 2:
165 continue
167 a = np.dot(u, atoms.get_cell())
168 b = np.dot(v, atoms.get_cell())
169 norm = np.linalg.norm(np.cross(a, b))
170 best = min(best, (norm, a, b), key=lambda x: x[0])
171 _, a, b = best
172 return a, b, orthogonal_basis(a, b)
175def isolate_monolayer(atoms, components, k, v):
177 a, b, basis = construct_inplane_basis(atoms, k, v)
179 # project atoms into new basis
180 c = np.cross(a, b)
181 c /= np.linalg.norm(c)
182 init_cell = np.dot(np.array([a, b, c]), basis.T)
184 positions, numbers = build_supercomponent(atoms, components, k, v)
185 scaled = np.linalg.solve(init_cell.T, np.dot(positions, basis.T).T).T
187 # move atoms into new cell
188 scaled[:, :2] %= 1.0
190 # subtract barycentre in z-direction
191 scaled[:, 2] -= np.mean(scaled, axis=0)[2]
193 # make cell large enough in z-direction
194 pos = np.dot(scaled, init_cell)
195 zmax = np.max(np.abs(pos[:, 2]))
196 cell = np.copy(init_cell)
197 cell[2] *= 4 * zmax
199 # construct a new atoms object containing the isolated chain
200 return Atoms(numbers=numbers, positions=pos, cell=cell, pbc=[1, 1, 0])
203def isolate_bulk(atoms, components, k, v):
204 positions, numbers = build_supercomponent(atoms, components, k, v,
205 anchor=False)
206 atoms = Atoms(numbers=numbers, positions=positions, cell=atoms.cell,
207 pbc=[1, 1, 1])
208 atoms.wrap(eps=0)
209 return atoms
212def isolate_cluster(atoms, components, k, v):
213 positions, numbers = build_supercomponent(atoms, components, k, v)
214 positions -= np.min(positions, axis=0)
215 cell = np.diag(np.max(positions, axis=0))
216 return Atoms(numbers=numbers, positions=positions, cell=cell, pbc=False)
219def isolate_components(atoms, kcutoff=None):
220 """Isolates components by dimensionality type.
222 Given a k-value cutoff the components (connected clusters) are
223 identified. For each component an Atoms object is created, which contains
224 that component only. The geometry of the resulting Atoms object depends
225 on the component dimensionality type:
227 0D: The cell is a tight box around the atoms. pbc=[0, 0, 0].
228 The cell has no physical meaning.
230 1D: The chain is aligned along the z-axis. pbc=[0, 0, 1].
231 The x and y cell directions have no physical meaning.
233 2D: The layer is aligned in the x-y plane. pbc=[1, 1, 0].
234 The z cell direction has no physical meaning.
236 3D: The original cell is used. pbc=[1, 1, 1].
238 Parameters:
240 atoms: ASE atoms object
241 The system to analyze.
242 kcutoff: float
243 The k-value cutoff to use. Default=None, in which case the
244 dimensionality scoring parameter is used to select the cutoff.
246 Returns:
248 components: dict
249 key: the component dimenionalities.
250 values: a list of Atoms objects for each dimensionality type.
251 """
252 data = {}
253 components, all_visited, ranks = traverse_graph(atoms, kcutoff)
255 for k, v in all_visited.items():
256 v = sorted(list(v))
258 # identify the components which constitute the component
259 key = tuple(np.unique([c for c, offset in v]))
260 dim = ranks[k]
262 if dim == 0:
263 data[('0D', key)] = isolate_cluster(atoms, components, k, v)
264 elif dim == 1:
265 data[('1D', key)] = isolate_chain(atoms, components, k, v)
266 elif dim == 2:
267 data[('2D', key)] = isolate_monolayer(atoms, components, k, v)
268 elif dim == 3:
269 data[('3D', key)] = isolate_bulk(atoms, components, k, v)
271 result = collections.defaultdict(list)
272 for (dim, _), atoms in data.items():
273 result[dim].append(atoms)
274 return result