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]