import numpy as np
[docs]def cut(atoms, a=(1, 0, 0), b=(0, 1, 0), c=None, clength=None,
origo=(0, 0, 0), nlayers=None, extend=1.0, tolerance=0.01,
maxatoms=None):
"""Cuts out a cell defined by *a*, *b*, *c* and *origo* from a
sufficiently repeated copy of *atoms*.
Typically, this function is used to create slabs of different
sizes and orientations. The vectors *a*, *b* and *c* are in scaled
coordinates and defines the returned cell and should normally be
integer-valued in order to end up with a periodic
structure. However, for systems with sub-translations, like fcc,
integer multiples of 1/2 or 1/3 might also make sense for some
directions (and will be treated correctly).
Parameters:
atoms: Atoms instance
This should correspond to a repeatable unit cell.
a: int | 3 floats
The a-vector in scaled coordinates of the cell to cut out. If
integer, the a-vector will be the scaled vector from *origo* to the
atom with index *a*.
b: int | 3 floats
The b-vector in scaled coordinates of the cell to cut out. If
integer, the b-vector will be the scaled vector from *origo* to the
atom with index *b*.
c: None | int | 3 floats
The c-vector in scaled coordinates of the cell to cut out.
if integer, the c-vector will be the scaled vector from *origo* to
the atom with index *c*.
If *None* it will be along cross(a, b) converted to real space
and normalised with the cube root of the volume. Note that this
in general is not perpendicular to a and b for non-cubic
systems. For cubic systems however, this is redused to
c = cross(a, b).
clength: None | float
If not None, the length of the c-vector will be fixed to
*clength* Angstroms. Should not be used together with
*nlayers*.
origo: int | 3 floats
Position of origo of the new cell in scaled coordinates. If
integer, the position of the atom with index *origo* is used.
nlayers: None | int
If *nlayers* is not *None*, the returned cell will have
*nlayers* atomic layers in the c-direction.
extend: 1 or 3 floats
The *extend* argument scales the effective cell in which atoms
will be included. It must either be three floats or a single
float scaling all 3 directions. By setting to a value just
above one, e.g. 1.05, it is possible to all the corner and
edge atoms in the returned cell. This will of cause make the
returned cell non-repeatable, but is very useful for
visualisation.
tolerance: float
Determines what is defined as a plane. All atoms within
*tolerance* Angstroms from a given plane will be considered to
belong to that plane.
maxatoms: None | int
This option is used to auto-tune *tolerance* when *nlayers* is
given for high zone axis systems. For high zone axis one
needs to reduce *tolerance* in order to distinguise the atomic
planes, resulting in the more atoms will be added and
eventually MemoryError. A too small *tolerance*, on the other
hand, might result in inproper splitting of atomic planes and
that too few layers are returned. If *maxatoms* is not None,
*tolerance* will automatically be gradually reduced until
*nlayers* atomic layers is obtained, when the number of atoms
exceeds *maxatoms*.
Example:
>>> import ase
>>> from ase.spacegroup import crystal
>>>
# Create an aluminium (111) slab with three layers
#
# First an unit cell of Al
>>> a = 4.05
>>> aluminium = crystal('Al', [(0,0,0)], spacegroup=225,
... cellpar=[a, a, a, 90, 90, 90])
>>>
# Then cut out the slab
>>> al111 = cut(aluminium, (1,-1,0), (0,1,-1), nlayers=3)
>>>
# Visualisation of the skutterudite unit cell
#
# Again, create a skutterudite unit cell
>>> a = 9.04
>>> skutterudite = crystal(
... ('Co', 'Sb'),
... basis=[(0.25,0.25,0.25), (0.0, 0.335, 0.158)],
... spacegroup=204,
... cellpar=[a, a, a, 90, 90, 90])
>>>
# Then use *origo* to put 'Co' at the corners and *extend* to
# include all corner and edge atoms.
>>> s = cut(skutterudite, origo=(0.25, 0.25, 0.25), extend=1.01)
>>> ase.view(s) # doctest: +SKIP
"""
atoms = atoms.copy()
cell = atoms.cell
if isinstance(origo, int):
origo = atoms.get_scaled_positions()[origo]
origo = np.array(origo, dtype=float)
scaled = (atoms.get_scaled_positions() - origo) % 1.0
scaled %= 1.0 # needed to ensure that all numbers are *less* than one
atoms.set_scaled_positions(scaled)
if isinstance(a, int):
a = scaled[a] - origo
if isinstance(b, int):
b = scaled[b] - origo
if isinstance(c, int):
c = scaled[c] - origo
a = np.array(a, dtype=float)
b = np.array(b, dtype=float)
if c is None:
metric = np.dot(cell, cell.T)
vol = np.sqrt(np.linalg.det(metric))
h = np.cross(a, b)
H = np.linalg.solve(metric.T, h.T)
c = vol * H / vol**(1. / 3.)
c = np.array(c, dtype=float)
if nlayers:
# Recursive increase the length of c until we have at least
# *nlayers* atomic layers parallel to the a-b plane
while True:
at = cut(atoms, a, b, c, origo=origo, extend=extend,
tolerance=tolerance)
scaled = at.get_scaled_positions()
d = scaled[:, 2]
keys = np.argsort(d)
ikeys = np.argsort(keys)
tol = tolerance
while True:
mask = np.concatenate(([True], np.diff(d[keys]) > tol))
tags = np.cumsum(mask)[ikeys] - 1
levels = d[keys][mask]
if (maxatoms is None or len(at) < maxatoms or
len(levels) > nlayers):
break
tol *= 0.9
if len(levels) > nlayers:
break
c *= 2
at.cell[2] *= levels[nlayers]
return at[tags < nlayers]
newcell = np.dot(np.array([a, b, c]), cell)
if nlayers is None and clength is not None:
newcell[2, :] *= clength / np.linalg.norm(newcell[2])
# Create a new atoms object, repeated and translated such that
# it completely covers the new cell
scorners_newcell = np.array([[0., 0., 0.], [0., 0., 1.],
[0., 1., 0.], [0., 1., 1.],
[1., 0., 0.], [1., 0., 1.],
[1., 1., 0.], [1., 1., 1.]])
corners = np.dot(scorners_newcell, newcell * extend)
scorners = np.linalg.solve(cell.T, corners.T).T
rep = np.ceil(scorners.ptp(axis=0)).astype('int') + 1
trans = np.dot(np.floor(scorners.min(axis=0)), cell)
atoms = atoms.repeat(rep)
atoms.translate(trans)
atoms.set_cell(newcell)
# Mask out atoms outside new cell
stol = 0.1 * tolerance # scaled tolerance, XXX
maskcell = atoms.cell * extend
sp = np.linalg.solve(maskcell.T, (atoms.positions).T).T
mask = np.all(np.logical_and(-stol <= sp, sp < 1 - stol), axis=1)
atoms = atoms[mask]
return atoms
class IncompatibleCellError(ValueError):
"""Exception raised if stacking fails due to incompatible cells
between *atoms1* and *atoms2*."""
pass
[docs]def stack(atoms1, atoms2, axis=2, cell=None, fix=0.5,
maxstrain=0.5, distance=None, reorder=False,
output_strained=False):
"""Return a new Atoms instance with *atoms2* stacked on top of
*atoms1* along the given axis. Periodicity in all directions is
ensured.
The size of the final cell is determined by *cell*, except
that the length alongh *axis* will be the sum of
*atoms1.cell[axis]* and *atoms2.cell[axis]*. If *cell* is None,
it will be interpolated between *atoms1* and *atoms2*, where
*fix* determines their relative weight. Hence, if *fix* equals
zero, the final cell will be determined purely from *atoms1* and
if *fix* equals one, it will be determined purely from
*atoms2*.
An ase.geometry.IncompatibleCellError exception is raised if the
cells of *atoms1* and *atoms2* are incompatible, e.g. if the far
corner of the unit cell of either *atoms1* or *atoms2* is
displaced more than *maxstrain*. Setting *maxstrain* to None
disables this check.
If *distance* is not None, the size of the final cell, along the
direction perpendicular to the interface, will be adjusted such
that the distance between the closest atoms in *atoms1* and
*atoms2* will be equal to *distance*. This option uses
scipy.optimize.fmin() and hence require scipy to be installed.
If *reorder* is True, then the atoms will be reordered such that
all atoms with the same symbol will follow sequencially after each
other, eg: 'Al2MnAl10Fe' -> 'Al12FeMn'.
If *output_strained* is True, then the strained versions of
*atoms1* and *atoms2* are returned in addition to the stacked
structure.
Example:
>>> import ase
>>> from ase.spacegroup import crystal
>>>
# Create an Ag(110)-Si(110) interface with three atomic layers
# on each side.
>>> a_ag = 4.09
>>> ag = crystal(['Ag'], basis=[(0,0,0)], spacegroup=225,
... cellpar=[a_ag, a_ag, a_ag, 90., 90., 90.])
>>> ag110 = cut(ag, (0, 0, 3), (-1.5, 1.5, 0), nlayers=3)
>>>
>>> a_si = 5.43
>>> si = crystal(['Si'], basis=[(0,0,0)], spacegroup=227,
... cellpar=[a_si, a_si, a_si, 90., 90., 90.])
>>> si110 = cut(si, (0, 0, 2), (-1, 1, 0), nlayers=3)
>>>
>>> interface = stack(ag110, si110, maxstrain=1)
>>> ase.view(interface) # doctest: +SKIP
>>>
# Once more, this time adjusted such that the distance between
# the closest Ag and Si atoms will be 2.3 Angstrom (requires scipy).
>>> interface2 = stack(ag110, si110,
... maxstrain=1, distance=2.3) # doctest:+ELLIPSIS
Optimization terminated successfully.
...
>>> ase.view(interface2) # doctest: +SKIP
"""
atoms1 = atoms1.copy()
atoms2 = atoms2.copy()
for atoms in [atoms1, atoms2]:
if not atoms.cell[axis].any():
atoms.center(vacuum=0.0, axis=axis)
if (np.sign(np.linalg.det(atoms1.cell)) !=
np.sign(np.linalg.det(atoms2.cell))):
raise IncompatibleCellError('Cells of *atoms1* and *atoms2* must have '
'same handedness.')
c1 = np.linalg.norm(atoms1.cell[axis])
c2 = np.linalg.norm(atoms2.cell[axis])
if cell is None:
cell1 = atoms1.cell.copy()
cell2 = atoms2.cell.copy()
cell1[axis] /= c1
cell2[axis] /= c2
cell = cell1 + fix * (cell2 - cell1)
cell[axis] /= np.linalg.norm(cell[axis])
cell1 = cell.copy()
cell2 = cell.copy()
cell1[axis] *= c1
cell2[axis] *= c2
if maxstrain:
strain1 = np.sqrt(((cell1 - atoms1.cell).sum(axis=0)**2).sum())
strain2 = np.sqrt(((cell2 - atoms2.cell).sum(axis=0)**2).sum())
if strain1 > maxstrain or strain2 > maxstrain:
raise IncompatibleCellError(
'*maxstrain* exceeded. *atoms1* strained %f and '
'*atoms2* strained %f.' % (strain1, strain2))
atoms1.set_cell(cell1, scale_atoms=True)
atoms2.set_cell(cell2, scale_atoms=True)
if output_strained:
atoms1_strained = atoms1.copy()
atoms2_strained = atoms2.copy()
if distance is not None:
from scipy.optimize import fmin
def mindist(pos1, pos2):
n1 = len(pos1)
n2 = len(pos2)
idx1 = np.arange(n1).repeat(n2)
idx2 = np.tile(np.arange(n2), n1)
return np.sqrt(((pos1[idx1] - pos2[idx2])**2).sum(axis=1).min())
def func(x):
t1, t2, h1, h2 = x[0:3], x[3:6], x[6], x[7]
pos1 = atoms1.positions + t1
pos2 = atoms2.positions + t2
d1 = mindist(pos1, pos2 + (h1 + 1.0) * atoms1.cell[axis])
d2 = mindist(pos2, pos1 + (h2 + 1.0) * atoms2.cell[axis])
return (d1 - distance)**2 + (d2 - distance)**2
atoms1.center()
atoms2.center()
x0 = np.zeros((8,))
x = fmin(func, x0)
t1, t2, h1, h2 = x[0:3], x[3:6], x[6], x[7]
atoms1.translate(t1)
atoms2.translate(t2)
atoms1.cell[axis] *= 1.0 + h1
atoms2.cell[axis] *= 1.0 + h2
atoms2.translate(atoms1.cell[axis])
atoms1.cell[axis] += atoms2.cell[axis]
atoms1.extend(atoms2)
if reorder:
atoms1 = sort(atoms1)
if output_strained:
return atoms1, atoms1_strained, atoms2_strained
else:
return atoms1
def rotation_matrix(a1, a2, b1, b2):
"""Returns a rotation matrix that rotates the vectors *a1* in the
direction of *a2* and *b1* in the direction of *b2*.
In the case that the angle between *a2* and *b2* is not the same
as between *a1* and *b1*, a proper rotation matrix will anyway be
constructed by first rotate *b2* in the *b1*, *b2* plane.
"""
a1 = np.asarray(a1, dtype=float) / np.linalg.norm(a1)
b1 = np.asarray(b1, dtype=float) / np.linalg.norm(b1)
c1 = np.cross(a1, b1)
c1 /= np.linalg.norm(c1) # clean out rounding errors...
a2 = np.asarray(a2, dtype=float) / np.linalg.norm(a2)
b2 = np.asarray(b2, dtype=float) / np.linalg.norm(b2)
c2 = np.cross(a2, b2)
c2 /= np.linalg.norm(c2) # clean out rounding errors...
# Calculate rotated *b2*
theta = np.arccos(np.dot(a2, b2)) - np.arccos(np.dot(a1, b1))
b3 = np.sin(theta) * a2 + np.cos(theta) * b2
b3 /= np.linalg.norm(b3) # clean out rounding errors...
A1 = np.array([a1, b1, c1])
A2 = np.array([a2, b3, c2])
R = np.linalg.solve(A1, A2).T
return R
[docs]def rotate(atoms, a1, a2, b1, b2, rotate_cell=True, center=(0, 0, 0)):
"""Rotate *atoms*, such that *a1* will be rotated in the direction
of *a2* and *b1* in the direction of *b2*. The point at *center*
is fixed. Use *center='COM'* to fix the center of mass. If
*rotate_cell* is true, the cell will be rotated together with the
atoms.
Note that the 000-corner of the cell is by definition fixed at
origo. Hence, setting *center* to something other than (0, 0, 0)
will rotate the atoms out of the cell, even if *rotate_cell* is
True.
"""
if isinstance(center, str) and center.lower() == 'com':
center = atoms.get_center_of_mass()
R = rotation_matrix(a1, a2, b1, b2)
atoms.positions[:] = np.dot(atoms.positions - center, R.T) + center
if rotate_cell:
atoms.cell[:] = np.dot(atoms.cell, R.T)
def minimize_tilt_ij(atoms, modified=1, fixed=0, fold_atoms=True):
"""Minimize the tilt angle for two given axes.
The problem is underdetermined. Therefore one can choose one axis
that is kept fixed.
"""
orgcell_cc = atoms.get_cell()
pbc_c = atoms.get_pbc()
i = fixed
j = modified
if not (pbc_c[i] and pbc_c[j]):
raise RuntimeError('Axes have to be periodic')
prod_cc = np.dot(orgcell_cc, orgcell_cc.T)
cell_cc = 1. * orgcell_cc
nji = np.floor(- prod_cc[i, j] / prod_cc[i, i] + 0.5)
cell_cc[j] = orgcell_cc[j] + nji * cell_cc[i]
# sanity check
def volume(cell):
return np.abs(np.dot(cell[2], np.cross(cell[0], cell[1])))
V = volume(cell_cc)
assert(abs(volume(orgcell_cc) - V) / V < 1.e-10)
atoms.set_cell(cell_cc)
if fold_atoms:
atoms.wrap()
[docs]def minimize_tilt(atoms, order=range(3), fold_atoms=True):
"""Minimize the tilt angles of the unit cell."""
pbc_c = atoms.get_pbc()
for i1, c1 in enumerate(order):
for c2 in order[i1 + 1:]:
if pbc_c[c1] and pbc_c[c2]:
minimize_tilt_ij(atoms, c1, c2, fold_atoms)
def niggli_reduce_cell(cell, epsfactor=None):
from ase.geometry import cellpar_to_cell
if epsfactor is None:
epsfactor = 1e-5
eps = epsfactor * abs(np.linalg.det(cell))**(1./3.)
cell = np.asarray(cell)
I3 = np.eye(3, dtype=int)
I6 = np.eye(6, dtype=int)
C = I3.copy()
D = I6.copy()
g0 = np.zeros(6, dtype=float)
g0[0] = np.dot(cell[0], cell[0])
g0[1] = np.dot(cell[1], cell[1])
g0[2] = np.dot(cell[2], cell[2])
g0[3] = 2 * np.dot(cell[1], cell[2])
g0[4] = 2 * np.dot(cell[0], cell[2])
g0[5] = 2 * np.dot(cell[0], cell[1])
g = np.dot(D, g0)
def lt(x, y, eps=eps):
return x < y - eps
def gt(x, y, eps=eps):
return lt(y, x, eps)
def eq(x, y, eps=eps):
return not (lt(x, y, eps) or gt(x, y, eps))
for _ in range(10000):
if (gt(g[0], g[1])
or (eq(g[0], g[1]) and gt(abs(g[3]), abs(g[4])))):
C = np.dot(C, -I3[[1, 0, 2]])
D = np.dot(I6[[1, 0, 2, 4, 3, 5]], D)
g = np.dot(D, g0)
continue
elif (gt(g[1], g[2])
or (eq(g[1], g[2]) and gt(abs(g[4]), abs(g[5])))):
C = np.dot(C, -I3[[0, 2, 1]])
D = np.dot(I6[[0, 2, 1, 3, 5, 4]], D)
g = np.dot(D, g0)
continue
lmn = np.array(gt(g[3:], 0, eps=eps/2), dtype=int)
lmn -= np.array(lt(g[3:], 0, eps=eps/2), dtype=int)
if lmn.prod() == 1:
ijk = lmn.copy()
for idx in range(3):
if ijk[idx] == 0:
ijk[idx] = 1
else:
ijk = np.ones(3, dtype=int)
if np.any(lmn != -1):
r = None
for idx in range(3):
if lmn[idx] == 1:
ijk[idx] = -1
elif lmn[idx] == 0:
r = idx
if ijk.prod() == -1:
ijk[r] = -1
C *= ijk[np.newaxis]
D[3] *= ijk[1] * ijk[2]
D[4] *= ijk[0] * ijk[2]
D[5] *= ijk[0] * ijk[1]
g = np.dot(D, g0)
if (gt(abs(g[3]), g[1])
or (eq(g[3], g[1]) and lt(2 * g[4], g[5]))
or (eq(g[3], -g[1]) and lt(g[5], 0))):
s = int(np.sign(g[3]))
A = I3.copy()
A[1, 2] = -s
C = np.dot(C, A)
B = I6.copy()
B[2, 1] = 1
B[2, 3] = -s
B[3, 1] = -2 * s
B[4, 5] = -s
D = np.dot(B, D)
g = np.dot(D, g0)
elif (gt(abs(g[4]), g[0])
or (eq(g[4], g[0]) and lt(2 * g[3], g[5]))
or (eq(g[4], -g[0]) and lt(g[5], 0))):
s = int(np.sign(g[4]))
A = I3.copy()
A[0, 2] = -s
C = np.dot(C, A)
B = I6.copy()
B[2, 0] = 1
B[2, 4] = -s
B[3, 5] = -s
B[4, 0] = -2 * s
D = np.dot(B, D)
g = np.dot(D, g0)
elif (gt(abs(g[5]), g[0])
or (eq(g[5], g[0]) and lt(2 * g[3], g[4]))
or (eq(g[5], -g[0]) and lt(g[4], 0))):
s = int(np.sign(g[5]))
A = I3.copy()
A[0, 1] = -s
C = np.dot(C, A)
B = I6.copy()
B[1, 0] = 1
B[1, 5] = -s
B[3, 4] = -s
B[5, 0] = -2 * s
D = np.dot(B, D)
g = np.dot(D, g0)
elif (lt(g[[0, 1, 3, 4, 5]].sum(), 0)
or (eq(g[[0, 1, 3, 4, 5]].sum(), 0)
and gt(2 * (g[0] + g[4]) + g[5], 0))):
A = I3.copy()
A[:, 2] = 1
C = np.dot(C, A)
B = I6.copy()
B[2, :] = 1
B[3, 1] = 2
B[3, 5] = 1
B[4, 0] = 2
B[4, 5] = 1
D = np.dot(B, D)
g = np.dot(D, g0)
else:
break
else:
raise RuntimeError('Niggli reduction not done in 10000 steps!\n'
'cell={}\n'
'operation={}'
.format(cell.tolist(), C.tolist()))
abc = np.sqrt(g[:3])
# Prevent division by zero e.g. for cell==zeros((3, 3)):
abcprod = max(abc.prod(), 1e-100)
cosangles = abc * g[3:] / (2 * abcprod)
angles = 180 * np.arccos(cosangles) / np.pi
newcell = np.array(cellpar_to_cell(np.concatenate([abc, angles])),
dtype=float)
return newcell, C
def update_cell_and_positions(atoms, new_cell, op):
"""Helper method for transforming cell and positions of atoms object."""
scpos = np.linalg.solve(op, atoms.get_scaled_positions().T).T
scpos %= 1.0
scpos %= 1.0
atoms.set_cell(new_cell)
atoms.set_scaled_positions(scpos)
[docs]def niggli_reduce(atoms):
"""Convert the supplied atoms object's unit cell into its
maximally-reduced Niggli unit cell. Even if the unit cell is already
maximally reduced, it will be converted into its unique Niggli unit cell.
This will also wrap all atoms into the new unit cell.
References:
Niggli, P. "Krystallographische und strukturtheoretische Grundbegriffe.
Handbuch der Experimentalphysik", 1928, Vol. 7, Part 1, 108-176.
Krivy, I. and Gruber, B., "A Unified Algorithm for Determining the
Reduced (Niggli) Cell", Acta Cryst. 1976, A32, 297-298.
Grosse-Kunstleve, R.W.; Sauter, N. K.; and Adams, P. D. "Numerically
stable algorithms for the computation of reduced unit cells", Acta Cryst.
2004, A60, 1-6.
"""
assert all(atoms.pbc), 'Can only reduce 3d periodic unit cells!'
new_cell, op = niggli_reduce_cell(atoms.cell)
update_cell_and_positions(atoms, new_cell, op)
def reduce_lattice(atoms, eps=2e-4):
"""Reduce atoms object to canonical lattice.
This changes the cell and positions such that the atoms object has
the canonical form used for defining band paths but is otherwise
physically equivalent. The eps parameter is used as a tolerance
for determining the cell's Bravais lattice."""
from ase.geometry.bravais_type_engine import identify_lattice
niggli_reduce(atoms)
lat, op = identify_lattice(atoms.cell, eps=eps)
update_cell_and_positions(atoms, lat.tocell(), np.linalg.inv(op))
[docs]def sort(atoms, tags=None):
"""Return a new Atoms object with sorted atomic order. The default
is to order according to chemical symbols, but if *tags* is not
None, it will be used instead. A stable sorting algorithm is used.
Example:
>>> from ase.build import bulk
>>> # Two unit cells of NaCl:
>>> a = 5.64
>>> nacl = bulk('NaCl', 'rocksalt', a=a) * (2, 1, 1)
>>> nacl.get_chemical_symbols()
['Na', 'Cl', 'Na', 'Cl']
>>> nacl_sorted = sort(nacl)
>>> nacl_sorted.get_chemical_symbols()
['Cl', 'Cl', 'Na', 'Na']
>>> np.all(nacl_sorted.cell == nacl.cell)
True
"""
if tags is None:
tags = atoms.get_chemical_symbols()
else:
tags = list(tags)
deco = sorted([(tag, i) for i, tag in enumerate(tags)])
indices = [i for tag, i in deco]
return atoms[indices]