You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
656 lines
23 KiB
656 lines
23 KiB
import numpy as np |
|
|
|
|
|
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 |
|
|
|
|
|
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 |
|
|
|
|
|
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() |
|
|
|
|
|
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) |
|
|
|
|
|
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)) |
|
|
|
|
|
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]
|
|
|