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