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