Coverage for /builds/ase/ase/ase/build/tools.py: 72.99%
211 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import numpy as np
5from ase.build.niggli import niggli_reduce_cell
8def cut(atoms, a=(1, 0, 0), b=(0, 1, 0), c=None, clength=None,
9 origo=(0, 0, 0), nlayers=None, extend=1.0, tolerance=0.01,
10 maxatoms=None):
11 """Cuts out a cell defined by *a*, *b*, *c* and *origo* from a
12 sufficiently repeated copy of *atoms*.
14 Typically, this function is used to create slabs of different
15 sizes and orientations. The vectors *a*, *b* and *c* are in scaled
16 coordinates and defines the returned cell and should normally be
17 integer-valued in order to end up with a periodic
18 structure. However, for systems with sub-translations, like fcc,
19 integer multiples of 1/2 or 1/3 might also make sense for some
20 directions (and will be treated correctly).
22 Parameters:
24 atoms: Atoms instance
25 This should correspond to a repeatable unit cell.
26 a: int | 3 floats
27 The a-vector in scaled coordinates of the cell to cut out. If
28 integer, the a-vector will be the scaled vector from *origo* to the
29 atom with index *a*.
30 b: int | 3 floats
31 The b-vector in scaled coordinates of the cell to cut out. If
32 integer, the b-vector will be the scaled vector from *origo* to the
33 atom with index *b*.
34 c: None | int | 3 floats
35 The c-vector in scaled coordinates of the cell to cut out.
36 if integer, the c-vector will be the scaled vector from *origo* to
37 the atom with index *c*.
38 If *None* it will be along cross(a, b) converted to real space
39 and normalised with the cube root of the volume. Note that this
40 in general is not perpendicular to a and b for non-cubic
41 systems. For cubic systems however, this is redused to
42 c = cross(a, b).
43 clength: None | float
44 If not None, the length of the c-vector will be fixed to
45 *clength* Angstroms. Should not be used together with
46 *nlayers*.
47 origo: int | 3 floats
48 Position of origo of the new cell in scaled coordinates. If
49 integer, the position of the atom with index *origo* is used.
50 nlayers: None | int
51 If *nlayers* is not *None*, the returned cell will have
52 *nlayers* atomic layers in the c-direction.
53 extend: 1 or 3 floats
54 The *extend* argument scales the effective cell in which atoms
55 will be included. It must either be three floats or a single
56 float scaling all 3 directions. By setting to a value just
57 above one, e.g. 1.05, it is possible to all the corner and
58 edge atoms in the returned cell. This will of cause make the
59 returned cell non-repeatable, but is very useful for
60 visualisation.
61 tolerance: float
62 Determines what is defined as a plane. All atoms within
63 *tolerance* Angstroms from a given plane will be considered to
64 belong to that plane.
65 maxatoms: None | int
66 This option is used to auto-tune *tolerance* when *nlayers* is
67 given for high zone axis systems. For high zone axis one
68 needs to reduce *tolerance* in order to distinguise the atomic
69 planes, resulting in the more atoms will be added and
70 eventually MemoryError. A too small *tolerance*, on the other
71 hand, might result in inproper splitting of atomic planes and
72 that too few layers are returned. If *maxatoms* is not None,
73 *tolerance* will automatically be gradually reduced until
74 *nlayers* atomic layers is obtained, when the number of atoms
75 exceeds *maxatoms*.
77 Example: Create an aluminium (111) slab with three layers.
79 >>> import ase
80 >>> from ase.spacegroup import crystal
81 >>> from ase.build.tools import cut
83 # First, a unit cell of Al
84 >>> a = 4.05
85 >>> aluminium = crystal('Al', [(0,0,0)], spacegroup=225,
86 ... cellpar=[a, a, a, 90, 90, 90])
88 # Then cut out the slab
89 >>> al111 = cut(aluminium, (1,-1,0), (0,1,-1), nlayers=3)
91 Example: Visualisation of the skutterudite unit cell
93 >>> from ase.spacegroup import crystal
94 >>> from ase.build.tools import cut
96 # Again, create a skutterudite unit cell
97 >>> a = 9.04
98 >>> skutterudite = crystal(
99 ... ('Co', 'Sb'),
100 ... basis=[(0.25,0.25,0.25), (0.0, 0.335, 0.158)],
101 ... spacegroup=204,
102 ... cellpar=[a, a, a, 90, 90, 90])
104 # Then use *origo* to put 'Co' at the corners and *extend* to
105 # include all corner and edge atoms.
106 >>> s = cut(skutterudite, origo=(0.25, 0.25, 0.25), extend=1.01)
107 >>> ase.view(s) # doctest:+SKIP
108 """
109 atoms = atoms.copy()
110 cell = atoms.cell
112 if isinstance(origo, int):
113 origo = atoms.get_scaled_positions()[origo]
114 origo = np.array(origo, dtype=float)
116 scaled = (atoms.get_scaled_positions() - origo) % 1.0
117 scaled %= 1.0 # needed to ensure that all numbers are *less* than one
118 atoms.set_scaled_positions(scaled)
120 if isinstance(a, int):
121 a = scaled[a] - origo
122 if isinstance(b, int):
123 b = scaled[b] - origo
124 if isinstance(c, int):
125 c = scaled[c] - origo
127 a = np.array(a, dtype=float)
128 b = np.array(b, dtype=float)
129 if c is None:
130 metric = np.dot(cell, cell.T)
131 vol = np.sqrt(np.linalg.det(metric))
132 h = np.cross(a, b)
133 H = np.linalg.solve(metric.T, h.T)
134 c = vol * H / vol**(1. / 3.)
135 c = np.array(c, dtype=float)
137 if nlayers:
138 # Recursive increase the length of c until we have at least
139 # *nlayers* atomic layers parallel to the a-b plane
140 while True:
141 at = cut(atoms, a, b, c, origo=origo, extend=extend,
142 tolerance=tolerance)
143 scaled = at.get_scaled_positions()
144 d = scaled[:, 2]
145 keys = np.argsort(d)
146 ikeys = np.argsort(keys)
147 tol = tolerance
148 while True:
149 mask = np.concatenate(([True], np.diff(d[keys]) > tol))
150 tags = np.cumsum(mask)[ikeys] - 1
151 levels = d[keys][mask]
152 if (maxatoms is None or len(at) < maxatoms or
153 len(levels) > nlayers):
154 break
155 tol *= 0.9
156 if len(levels) > nlayers:
157 break
158 c *= 2
160 at.cell[2] *= levels[nlayers]
161 return at[tags < nlayers]
163 newcell = np.dot(np.array([a, b, c]), cell)
164 if nlayers is None and clength is not None:
165 newcell[2, :] *= clength / np.linalg.norm(newcell[2])
167 # Create a new atoms object, repeated and translated such that
168 # it completely covers the new cell
169 scorners_newcell = np.array([[0., 0., 0.], [0., 0., 1.],
170 [0., 1., 0.], [0., 1., 1.],
171 [1., 0., 0.], [1., 0., 1.],
172 [1., 1., 0.], [1., 1., 1.]])
173 corners = np.dot(scorners_newcell, newcell * extend)
174 scorners = np.linalg.solve(cell.T, corners.T).T
175 rep = np.ceil(np.ptp(scorners, axis=0)).astype('int') + 1
176 trans = np.dot(np.floor(scorners.min(axis=0)), cell)
177 atoms = atoms.repeat(rep)
178 atoms.translate(trans)
179 atoms.set_cell(newcell)
181 # Mask out atoms outside new cell
182 stol = 0.1 * tolerance # scaled tolerance, XXX
183 maskcell = atoms.cell * extend
184 sp = np.linalg.solve(maskcell.T, (atoms.positions).T).T
185 mask = np.all(np.logical_and(-stol <= sp, sp < 1 - stol), axis=1)
186 atoms = atoms[mask]
187 return atoms
190class IncompatibleCellError(ValueError):
191 """Exception raised if stacking fails due to incompatible cells
192 between *atoms1* and *atoms2*."""
195def stack(atoms1, atoms2, axis=2, cell=None, fix=0.5,
196 maxstrain=0.5, distance=None, reorder=False,
197 output_strained=False):
198 """Return a new Atoms instance with *atoms2* stacked on top of
199 *atoms1* along the given axis. Periodicity in all directions is
200 ensured.
202 The size of the final cell is determined by *cell*, except
203 that the length alongh *axis* will be the sum of
204 *atoms1.cell[axis]* and *atoms2.cell[axis]*. If *cell* is None,
205 it will be interpolated between *atoms1* and *atoms2*, where
206 *fix* determines their relative weight. Hence, if *fix* equals
207 zero, the final cell will be determined purely from *atoms1* and
208 if *fix* equals one, it will be determined purely from
209 *atoms2*.
211 An ase.geometry.IncompatibleCellError exception is raised if the
212 cells of *atoms1* and *atoms2* are incompatible, e.g. if the far
213 corner of the unit cell of either *atoms1* or *atoms2* is
214 displaced more than *maxstrain*. Setting *maxstrain* to None
215 disables this check.
217 If *distance* is not None, the size of the final cell, along the
218 direction perpendicular to the interface, will be adjusted such
219 that the distance between the closest atoms in *atoms1* and
220 *atoms2* will be equal to *distance*. This option uses
221 scipy.optimize.fmin() and hence require scipy to be installed.
223 If *reorder* is True, then the atoms will be reordered such that
224 all atoms with the same symbol will follow sequencially after each
225 other, eg: 'Al2MnAl10Fe' -> 'Al12FeMn'.
227 If *output_strained* is True, then the strained versions of
228 *atoms1* and *atoms2* are returned in addition to the stacked
229 structure.
231 Example: Create an Ag(110)-Si(110) interface with three atomic layers
232 on each side.
234 >>> import ase
235 >>> from ase.spacegroup import crystal
236 >>> from ase.build.tools import cut, stack
237 >>>
238 >>> a_ag = 4.09
239 >>> ag = crystal(['Ag'], basis=[(0,0,0)], spacegroup=225,
240 ... cellpar=[a_ag, a_ag, a_ag, 90., 90., 90.])
241 >>> ag110 = cut(ag, (0, 0, 3), (-1.5, 1.5, 0), nlayers=3)
242 >>>
243 >>> a_si = 5.43
244 >>> si = crystal(['Si'], basis=[(0,0,0)], spacegroup=227,
245 ... cellpar=[a_si, a_si, a_si, 90., 90., 90.])
246 >>> si110 = cut(si, (0, 0, 2), (-1, 1, 0), nlayers=3)
247 >>>
248 >>> interface = stack(ag110, si110, maxstrain=1)
249 >>> ase.view(interface) # doctest: +SKIP
250 >>>
251 # Once more, this time adjusted such that the distance between
252 # the closest Ag and Si atoms will be 2.3 Angstrom (requires scipy).
253 >>> interface2 = stack(ag110, si110,
254 ... maxstrain=1, distance=2.3) # doctest:+ELLIPSIS
255 Optimization terminated successfully.
256 ...
257 >>> ase.view(interface2) # doctest: +SKIP
258 """
259 atoms1 = atoms1.copy()
260 atoms2 = atoms2.copy()
262 for atoms in [atoms1, atoms2]:
263 if not atoms.cell[axis].any():
264 atoms.center(vacuum=0.0, axis=axis)
266 if (np.sign(np.linalg.det(atoms1.cell)) !=
267 np.sign(np.linalg.det(atoms2.cell))):
268 raise IncompatibleCellError('Cells of *atoms1* and *atoms2* must have '
269 'same handedness.')
271 c1 = np.linalg.norm(atoms1.cell[axis])
272 c2 = np.linalg.norm(atoms2.cell[axis])
273 if cell is None:
274 cell1 = atoms1.cell.copy()
275 cell2 = atoms2.cell.copy()
276 cell1[axis] /= c1
277 cell2[axis] /= c2
278 cell = cell1 + fix * (cell2 - cell1)
279 cell[axis] /= np.linalg.norm(cell[axis])
280 cell1 = cell.copy()
281 cell2 = cell.copy()
282 cell1[axis] *= c1
283 cell2[axis] *= c2
285 if maxstrain:
286 strain1 = np.sqrt(((cell1 - atoms1.cell).sum(axis=0)**2).sum())
287 strain2 = np.sqrt(((cell2 - atoms2.cell).sum(axis=0)**2).sum())
288 if strain1 > maxstrain or strain2 > maxstrain:
289 raise IncompatibleCellError(
290 '*maxstrain* exceeded. *atoms1* strained %f and '
291 '*atoms2* strained %f.' % (strain1, strain2))
293 atoms1.set_cell(cell1, scale_atoms=True)
294 atoms2.set_cell(cell2, scale_atoms=True)
295 if output_strained:
296 atoms1_strained = atoms1.copy()
297 atoms2_strained = atoms2.copy()
299 if distance is not None:
300 from scipy.optimize import fmin
302 def mindist(pos1, pos2):
303 n1 = len(pos1)
304 n2 = len(pos2)
305 idx1 = np.arange(n1).repeat(n2)
306 idx2 = np.tile(np.arange(n2), n1)
307 return np.sqrt(((pos1[idx1] - pos2[idx2])**2).sum(axis=1).min())
309 def func(x):
310 t1, t2, h1, h2 = x[0:3], x[3:6], x[6], x[7]
311 pos1 = atoms1.positions + t1
312 pos2 = atoms2.positions + t2
313 d1 = mindist(pos1, pos2 + (h1 + 1.0) * atoms1.cell[axis])
314 d2 = mindist(pos2, pos1 + (h2 + 1.0) * atoms2.cell[axis])
315 return (d1 - distance)**2 + (d2 - distance)**2
317 atoms1.center()
318 atoms2.center()
319 x0 = np.zeros((8,))
320 x = fmin(func, x0)
321 t1, t2, h1, h2 = x[0:3], x[3:6], x[6], x[7]
322 atoms1.translate(t1)
323 atoms2.translate(t2)
324 atoms1.cell[axis] *= 1.0 + h1
325 atoms2.cell[axis] *= 1.0 + h2
327 atoms2.translate(atoms1.cell[axis])
328 atoms1.cell[axis] += atoms2.cell[axis]
329 atoms1.extend(atoms2)
331 if reorder:
332 atoms1 = sort(atoms1)
334 if output_strained:
335 return atoms1, atoms1_strained, atoms2_strained
336 else:
337 return atoms1
340def rotation_matrix(a1, a2, b1, b2):
341 """Returns a rotation matrix that rotates the vectors *a1* in the
342 direction of *a2* and *b1* in the direction of *b2*.
344 In the case that the angle between *a2* and *b2* is not the same
345 as between *a1* and *b1*, a proper rotation matrix will anyway be
346 constructed by first rotate *b2* in the *b1*, *b2* plane.
347 """
348 a1 = np.asarray(a1, dtype=float) / np.linalg.norm(a1)
349 b1 = np.asarray(b1, dtype=float) / np.linalg.norm(b1)
350 c1 = np.cross(a1, b1)
351 c1 /= np.linalg.norm(c1) # clean out rounding errors...
353 a2 = np.asarray(a2, dtype=float) / np.linalg.norm(a2)
354 b2 = np.asarray(b2, dtype=float) / np.linalg.norm(b2)
355 c2 = np.cross(a2, b2)
356 c2 /= np.linalg.norm(c2) # clean out rounding errors...
358 # Calculate rotated *b2*
359 theta = np.arccos(np.dot(a2, b2)) - np.arccos(np.dot(a1, b1))
360 b3 = np.sin(theta) * a2 + np.cos(theta) * b2
361 b3 /= np.linalg.norm(b3) # clean out rounding errors...
363 A1 = np.array([a1, b1, c1])
364 A2 = np.array([a2, b3, c2])
365 R = np.linalg.solve(A1, A2).T
366 return R
369def rotate(atoms, a1, a2, b1, b2, rotate_cell=True, center=(0, 0, 0)):
370 """Rotate *atoms*, such that *a1* will be rotated in the direction
371 of *a2* and *b1* in the direction of *b2*. The point at *center*
372 is fixed. Use *center='COM'* to fix the center of mass. If
373 *rotate_cell* is true, the cell will be rotated together with the
374 atoms.
376 Note that the 000-corner of the cell is by definition fixed at
377 origo. Hence, setting *center* to something other than (0, 0, 0)
378 will rotate the atoms out of the cell, even if *rotate_cell* is
379 True.
380 """
381 if isinstance(center, str) and center.lower() == 'com':
382 center = atoms.get_center_of_mass()
384 R = rotation_matrix(a1, a2, b1, b2)
385 atoms.positions[:] = np.dot(atoms.positions - center, R.T) + center
387 if rotate_cell:
388 atoms.cell[:] = np.dot(atoms.cell, R.T)
391def minimize_tilt_ij(atoms, modified=1, fixed=0, fold_atoms=True):
392 """Minimize the tilt angle for two given axes.
394 The problem is underdetermined. Therefore one can choose one axis
395 that is kept fixed.
396 """
398 orgcell_cc = atoms.get_cell()
399 pbc_c = atoms.get_pbc()
400 i = fixed
401 j = modified
402 if not (pbc_c[i] and pbc_c[j]):
403 raise RuntimeError('Axes have to be periodic')
405 prod_cc = np.dot(orgcell_cc, orgcell_cc.T)
406 cell_cc = 1. * orgcell_cc
407 nji = np.floor(- prod_cc[i, j] / prod_cc[i, i] + 0.5)
408 cell_cc[j] = orgcell_cc[j] + nji * cell_cc[i]
410 # sanity check
411 def volume(cell):
412 return np.abs(np.dot(cell[2], np.cross(cell[0], cell[1])))
413 V = volume(cell_cc)
414 assert abs(volume(orgcell_cc) - V) / V < 1.e-10
416 atoms.set_cell(cell_cc)
418 if fold_atoms:
419 atoms.wrap()
422def minimize_tilt(atoms, order=range(3), fold_atoms=True):
423 """Minimize the tilt angles of the unit cell."""
424 pbc_c = atoms.get_pbc()
426 for i1, c1 in enumerate(order):
427 for c2 in order[i1 + 1:]:
428 if pbc_c[c1] and pbc_c[c2]:
429 minimize_tilt_ij(atoms, c1, c2, fold_atoms)
432def update_cell_and_positions(atoms, new_cell, op):
433 """Helper method for transforming cell and positions of atoms object."""
434 scpos = np.linalg.solve(op, atoms.get_scaled_positions().T).T
436 # We do this twice because -1e-20 % 1 == 1:
437 scpos[:, atoms.pbc] %= 1.0
438 scpos[:, atoms.pbc] %= 1.0
440 atoms.set_cell(new_cell)
441 atoms.set_scaled_positions(scpos)
444def niggli_reduce(atoms):
445 """Convert the supplied atoms object's unit cell into its
446 maximally-reduced Niggli unit cell. Even if the unit cell is already
447 maximally reduced, it will be converted into its unique Niggli unit cell.
448 This will also wrap all atoms into the new unit cell.
450 References:
452 Niggli, P. "Krystallographische und strukturtheoretische Grundbegriffe.
453 Handbuch der Experimentalphysik", 1928, Vol. 7, Part 1, 108-176.
455 Krivy, I. and Gruber, B., "A Unified Algorithm for Determining the
456 Reduced (Niggli) Cell", Acta Cryst. 1976, A32, 297-298.
458 Grosse-Kunstleve, R.W.; Sauter, N. K.; and Adams, P. D. "Numerically
459 stable algorithms for the computation of reduced unit cells", Acta Cryst.
460 2004, A60, 1-6.
461 """
462 from ase.geometry.geometry import permute_axes
464 # Make sure non-periodic cell vectors are orthogonal
465 non_periodic_cv = atoms.cell[~atoms.pbc]
466 periodic_cv = atoms.cell[atoms.pbc]
467 if not np.isclose(np.dot(non_periodic_cv, periodic_cv.T), 0).all():
468 raise ValueError('Non-orthogonal cell along non-periodic dimensions')
470 input_atoms = atoms
472 # Permute axes, such that the non-periodic are along the last dimensions,
473 # since niggli_reduce_cell will change the order of axes.
474 permutation = np.argsort(~atoms.pbc)
475 ipermutation = np.empty_like(permutation)
476 ipermutation[permutation] = np.arange(len(permutation))
477 atoms = permute_axes(atoms, permutation)
479 # Perform the Niggli reduction on the cell
480 nonpbc = ~atoms.pbc
481 uncompleted_cell = atoms.cell.uncomplete(atoms.pbc)
482 new_cell, op = niggli_reduce_cell(uncompleted_cell)
483 new_cell[nonpbc] = atoms.cell[nonpbc]
484 update_cell_and_positions(atoms, new_cell, op)
486 # Undo the prior permutation.
487 atoms = permute_axes(atoms, ipermutation)
488 input_atoms.cell[:] = atoms.cell
489 input_atoms.positions[:] = atoms.positions
492def reduce_lattice(atoms, eps=2e-4):
493 """Reduce atoms object to canonical lattice.
495 This changes the cell and positions such that the atoms object has
496 the canonical form used for defining band paths but is otherwise
497 physically equivalent. The eps parameter is used as a tolerance
498 for determining the cell's Bravais lattice."""
499 from ase.lattice import identify_lattice
500 niggli_reduce(atoms)
501 lat, op = identify_lattice(atoms.cell, eps=eps)
502 update_cell_and_positions(atoms, lat.tocell(), np.linalg.inv(op))
505def sort(atoms, tags=None):
506 """Return a new Atoms object with sorted atomic order. The default
507 is to order according to chemical symbols, but if *tags* is not
508 None, it will be used instead. A stable sorting algorithm is used.
510 Example:
512 >>> from ase.build import bulk
513 >>> from ase.build.tools import sort
514 >>> # Two unit cells of NaCl:
515 >>> a = 5.64
516 >>> nacl = bulk('NaCl', 'rocksalt', a=a) * (2, 1, 1)
517 >>> nacl.get_chemical_symbols()
518 ['Na', 'Cl', 'Na', 'Cl']
519 >>> nacl_sorted = sort(nacl)
520 >>> nacl_sorted.get_chemical_symbols()
521 ['Cl', 'Cl', 'Na', 'Na']
522 >>> np.all(nacl_sorted.cell == nacl.cell)
523 True
524 """
525 if tags is None:
526 tags = atoms.get_chemical_symbols()
527 else:
528 tags = list(tags)
529 deco = sorted([(tag, i) for i, tag in enumerate(tags)])
530 indices = [i for tag, i in deco]
531 return atoms[indices]