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