Quasi-crystal
=============

In this tutorial, we show how to construct the quasi-crystal, in which we also need to shift
twist, reshape and merge the cells. Taking bi-layer graphene quasi-crystal as an example, a
quasi-crystal with 12-fold symmetry is formed by twisting one layer by :math:`30^\circ` with respect
to the center of :math:`\mathbf{c} = 2/3\ \mathbf{a}_1 + 2/3\ \mathbf{a}_2`, where
:math:`\mathbf{a}_1` and :math:`\mathbf{a}_2` are the lattice vectors of the primitive cell of
fixed layer. The script can be found at ``examples/advanced/quasi_crystal.py``. We begin with
importing the packages:

.. code-block:: python

    import math
    from typing import Union, Tuple

    import numpy as np
    from numpy.linalg import norm

    import tbplas as tb


The ``main`` function for constructing the quasi-crystal is defined as:

.. code-block:: python

    def main():
        prim_cell = tb.make_graphene_diamond()
        dim = (33, 33, 1)
        angle = 30 / 180 * math.pi
        center = np.array((2./3, 2./3, 0))
        radius = 3.0
        shift = 0.3
        cell = make_quasi_crystal_pc(prim_cell, dim=dim, angle=angle, center=center,
                                     radius=radius, shift=shift)
        cell.plot(with_cells=False, with_orbitals=False, hop_as_arrows=False,
                  hop_eng_cutoff=0.3)

Firstly, we define the geometric parameters. We need a large cell to hold the quasi-crystal, whose
dimension is given in ``dim``. ``angle`` is the twisting angle and ``center`` is the fractional
coordinate of twisting center. The radius of the quasi-crystal is controlled by ``radius``, while
``shift`` specifies the interlayer distance. The function ``make_quasi_crystal_pc`` is defined as:

.. code-block:: python

    def make_quasi_crystal_pc(prim_cell: tb.PrimitiveCell,
                              dim: Tuple[int, int, int],
                              angle: float,
                              center: np.ndarray,
                              radius: float = 3.0,
                              shift: float = 0.3) -> tb.PrimitiveCell:
        """
        Create quasi-crystal at primitive level.

        :param prim_cell: primitive cell from which the quasi-crystal is built
        :param dim: dimension of the extended primitive cell
        :param angle: twisting angle in RADIAN
        :param center: fractional coordinate of the twisting center in the primitive cell
        :param radius: radius of quasi-crystal in nm
        :param shift: distance of shift along z-axis in NANOMETER
        :return: the quasi-crystal
        """
        # Get the Cartesian coordinate of rotation center
        center = np.array([dim[0]//2, dim[1]//2, 0]) + center
        center = np.matmul(center, prim_cell.lat_vec)

        # Build fixed and twisted layers
        layer_fixed = tb.extend_prim_cell(prim_cell, dim=dim)
        layer_twisted = tb.extend_prim_cell(prim_cell, dim=dim)

        # Rotate and shift twisted layer
        tb.spiral_prim_cell(layer_twisted, angle=angle, center=center,
                            shift=shift)

        # Remove unnecessary orbitals
        cutoff_pc(layer_fixed, center=center, radius=radius)
        cutoff_pc(layer_twisted, center=center, radius=radius)

        # Reset the lattice of twisted layer
        layer_twisted.reset_lattice(layer_fixed.lat_vec, layer_fixed.origin,
                                    unit=tb.NM, fix_orb=True)

        # Build inter-cell hopping terms
        inter_hop = tb.PCInterHopping(layer_fixed, layer_twisted)
        neighbors = tb.find_neighbors(layer_fixed, layer_twisted,
                                      a_max=0, b_max=0, max_distance=0.75)
        for term in neighbors:
            i, j = term.pair
            inter_hop.add_hopping(term.rn, i, j, calc_hop(term.rij))

        # Extend hopping terms
        extend_hop(layer_fixed)
        extend_hop(layer_twisted)
        final_cell = tb.merge_prim_cell(layer_fixed, layer_twisted, inter_hop)
        return final_cell

where we firstly get the Cartesian coordinate of the twisting center, then build the fixed and twisted
layers and twist the twisted layer with respect to the center, as is done in :ref:`hetero_model`.
The twisting operation is done by the :func:`.spiral_prim_cell` function, where the Cartesian
coordinate of the center is given in the ``center`` argument. Then we remove unnecessary orbitals to
produce a round quasi-crystal with finite radius. This is done by calling the ``cutoff_pc`` function,
which is defined as:

.. code-block:: python

    def cutoff_pc(prim_cell: tb.PrimitiveCell,
                  center: np.ndarray,
                  radius: float = 3.0) -> None:
        """
        Cutoff primitive cell up to given radius with respect to given center.

        :param prim_cell: supercell to cut
        :param center: Cartesian coordinate of center in nm
        :param radius: cutoff radius in nm
        :return: None. The incoming supercell is modified.
        """
        idx_remove = []
        orb_pos = prim_cell.orb_pos_nm
        for i, pos in enumerate(orb_pos):
            if norm(pos[:2] - center[:2]) > radius:
                idx_remove.append(i)
        prim_cell.remove_orbitals(idx_remove)
        prim_cell.trim()

where we loop over orbital positions to collect the indices of unnecessary orbitals, then call the
``remove_orbitals`` and ``trim`` functions. Before merging the fixed and twisted layers, we need to
reset the lattice vectors and origin of twisted layer to that of fixed layer by calling the
``reset_lattice`` function. Then we can merge them safely by calling the :func:`.merge_prim_cell`
function. Finally, we extend the hoppings and visualize the quasi-crystal, where the ``extend_hop``
function is defined in :ref:`hetero_model`. The output is shown in following figure:

.. figure:: images/quasi_crystal/quasi_crystal.png
    :align: center
    :scale: 50%

    Plot of the quasi-crystal formed from the incommensurate :math:`30^\circ` twisted bi-layer
    graphene with a radius of 3 nm.