Directory structure:
└── aimsgb/
    ├── COPYRIGHT.rst
    ├── README.rst
    ├── setup.cfg
    ├── setup.py
    ├── aimsgb/
    │   ├── __init__.py
    │   ├── agb.py
    │   ├── grain.py
    │   ├── grain_bound.py
    │   └── utils.py
    ├── html_notebooks/
    ├── tests/
    └── .github/
        └── workflows/
            └── testing.yml



================================================
FILE: README.rst
================================================
Introduction
============
aimsgb is an efficient and open-source Python library for generating atomic coordinates in periodic grain boundary models. It is designed to
construct various grain boundary structures from cubic and non-cubic initial
configurations. A convenient command line tool has also been provided to enable
easy and fast construction of tilt and twist boundaries by assigining the degree
of fit (Σ), rotation axis, grain boundary plane and initial crystal structure.

We also provide a web-based GUI to access aimsgb framework: `aimsgb.org
<https://aimsgb.org/>`_

Install aimsgb
==============
Method 1: Use Pip
-----------------
The easiest way to install aimsgb is to simply run a one-liner in pip::

   pip install aimsgb

Method 2: Use Git to install
----------------------------
1. Clone the latest version from github::

    git clone https://github.com/ksyang2013/aimsgb.git

2. Navigate to aimsgb folder::

    cd aimsgb

3. Type in the root of the repo::

    pip install .

4. or to install the package in development mode::

    pip install -e .


Usage
==================
Refer to the `documentation
<https://aimsgb-docs.readthedocs.io/en/latest/>`_ for more details.



================================================
FILE: aimsgb/__init__.py
================================================
from .grain import *
from .grain_bound import *


================================================
FILE: aimsgb/agb.py
================================================
#!/usr/bin/env python
from __future__ import division

import sys
import argparse
import numpy as np

from aimsgb import Grain, GBInformation, GrainBoundary

__author__ = "Jianli Cheng and Kesong YANG"
__copyright__ = "Copyright (C) 2018 The Regents of the University of California"
__maintainer__ = "Jianli Cheng"
__email__ = "jic198@ucsd.edu"
__status__ = "Production"
__date__ = "January 26, 2018"
aimsgb_descript = """**************************************************************************************************
       aimsgb - Ab-initio Interface Materials Simulation for Grain Boundaries (aimsgb.org)      
                        VERSION 0.1  2017 - 2018                                    
Aimsgb is one open-source python library to generate periodic grain boundary structures.    
A reference for the usage of aimsgb software is:
Jianli Cheng, Jian Luo, and Kesong Yang, Aimsgb: An Algorithm and Open-Source Python Library  
to Generate Periodic Grain Boundary Structures, Comput. Mater. Sci. 155, 92-103, (2018). 
**************************************************************************************************

Aimsgb Command Line Tools 
"""


def gb_list(args):
    axis = list(map(int, args.axis))
    sigma = args.sigma
    print(GBInformation(axis, sigma).__str__())


def gb(args):
    axis = list(map(int, args.axis))
    plane = args.plane
    initial_struct = None
    if args.filename:
        initial_struct = Grain.from_file(args.filename)
    elif args.mpid:
        initial_struct = Grain.from_mp_id(args.mpid)
    if initial_struct is None:
        raise ValueError("Please provide either filename or mpid.")
    
    gb = GrainBoundary(axis, args.sigma, plane, initial_struct, args.uc_a, args.uc_b)
    to_primitive = False if args.conventional else True
    structure = Grain.stack_grains(gb.grain_a, gb.grain_b, direction=gb.direction,
                                  to_primitive=to_primitive)
    structure.to(filename=args.out, fmt=args.fmt)
    print(f"CSL Matrix (det={np.linalg.det(gb.csl)}):\n{gb.csl}")
    print(f"{args.out} of sigma{gb.sigma}[{args.axis}]/({gb.plane_str}) is created")


def main():
    parser = argparse.ArgumentParser(description=aimsgb_descript,
                                     formatter_class=argparse.RawTextHelpFormatter)
    subparsers = parser.add_subparsers(help="command", dest="command")

    parser_gb_list = subparsers.add_parser(
        "list", help="Show the values of sigma, theta, GB plane and CSL matrix "
        "from a given rotation axis\nEXAMPLE: aimsgb list 001",
        formatter_class=argparse.RawTextHelpFormatter)
    parser_gb_list.add_argument("axis", metavar="rotation axis", type=str,
                                help="The rotation axis of GB\n"
                                "EXAMPLE: aimsgb list 001")
    parser_gb_list.add_argument("sigma", default=30, type=int, nargs="?",
                                help="Set the sigma limit for on screen output "
                                     "(default: %(default)s)\n"
                                     "EXAMPLE: aimsgb list 001 100")
    parser_gb_list.set_defaults(func=gb_list)

    parser_gb = subparsers.add_parser(
        "gb", help="Build grain boundary based on rotation axis, sigma, GB plane, "
                   "and input structure file.\nThe user can also specify many "
                   "optional arguments, such grain size and interface terminations."
                   "\nEXAMPLE1: aimsgb gb 001 5 1 2 0 -f POSCAR"
                   "\nEXAMPLE2: aimsgb gb 001 5 1 2 0 -mp mp-13 -ua 2 -ub 3"
                   "\nEXAMPLE3: aimsgb gb 001 5 1 2 0 -f POSCAR -ua 3 -ub 2 -dl 0b1t1b0t",
        formatter_class=argparse.RawTextHelpFormatter)
    parser_gb.add_argument("axis", metavar="rotation axis", type=str,
                           help="The rotation axis of GB, EXAMPLE: 110")
    parser_gb.add_argument("sigma", type=int,
                           help="The sigma value for grain boundary, EXAMPLE: 3")
    parser_gb.add_argument("plane", type=int, nargs=3,
                           help="The GB plane for grain boundary, EXAMPLE: 1 1 0")
    parser_gb.add_argument("-f", "--filename", type=str,
                           help="The initial structure file for grain boundary.")
    parser_gb.add_argument("-mp", "--mpid", type=str,
                           help="The mpid to get initial structure from Materials Project.")
    parser_gb.add_argument("-o", "--out", type=str, default="POSCAR", 
                           help="The output filename. (default: %(default)s)")
    parser_gb.add_argument("-ua", "--uc_a", default=1, type=int,
                           help="The size (uc) for grain A. (default: %(default)s)"
                                "\nEXAMPLE: aimsgb gb 001 5 1 2 0 -f POSCAR -ua 2")
    parser_gb.add_argument("-ub", "--uc_b", default=1, type=int,
                           help="The size (uc) for grain B. (default: %(default)s)"
                                "\nEXAMPLE: aimsgb gb 001 5 1 2 0 -f POSCAR -ub 2")
    parser_gb.add_argument("-dl", "--delete_layer", default="0b0t0b0t", type=str,
                           help="Delete bottom or top layers for each grain. "
                                "(default: %(default)s)"
                                "\nEXAMPLE: aimsgb gb 001 5 1 2 0 -f POSCAR -dl 0b1t1b0t")
    parser_gb.add_argument("-v", "--vacuum", default=0.0, type=float,
                           help="Set vacuum thickness for grain boundary. "
                                "(default: %(default)s angstrom)"
                                "\nEXAMPLE: aimsgb gb 001 5 1 2 0 -f POSCAR -v 20")
    parser_gb.add_argument("-t", "--tol", default=0.25, type=float,
                           help="Tolerance factor to determine if two "
                                "atoms are at the same plane. "
                                "(default: %(default)s angstrom)"
                                "\nEXAMPLE: aimsgb gb 001 5 1 2 0 -f POSCAR -dl 0b1t1b0t -t 0.5")
    parser_gb.add_argument("-ad", "--add_if_dist", default=0.0, type=float,
                           help="Add extra distance between two grains. "
                                "(default: %(default)s angstrom)"
                                "\nEXAMPLE: aimsgb gb 001 5 1 2 0 -f POSCAR -ad 0.5")
    parser_gb.add_argument("-c", "--conventional", action="store_true",
                           help="Get conventional GB, not primitive"
                                "\nEXAMPLE: aimsgb gb 001 5 1 2 0 -f POSCAR -c")
    parser_gb.add_argument("-fmt", "--fmt", default="poscar", const="poscar",
                           nargs="?", choices=["poscar", "cif", "cssr", "json"],
                           help="Choose the output format. (default: %(default)s)"
                                "\nEXAMPLE: aimsgb gb 001 5 1 2 0 -f POSCAR -fmt cif")
    parser_gb.set_defaults(func=gb)

    args = parser.parse_args()

    try:
        getattr(args, "func")
    except AttributeError:
        parser.print_help()
        sys.exit(0)
    args.func(args)


if __name__ == "__main__":
    main()



================================================
FILE: aimsgb/grain.py
================================================
import copy
import re
import warnings
import numpy as np
from numpy import sin, radians
from functools import reduce, wraps
from itertools import groupby
from aimsgb.utils import reduce_vector
from pymatgen.core.structure import Structure, Lattice, PeriodicSite
from pymatgen.transformations.advanced_transformations import CubicSupercellTransformation
from pymatgen.analysis.structure_matcher import StructureMatcher

__author__ = "Jianli CHENG and Kesong YANG"
__copyright__ = "Copyright 2018 University of California San Diego"
__maintainer__ = "Jianli CHENG"
__email__ = "jic198@ucsd.edu"
__status__ = "Production"
__date__ = "January 26, 2018"


class Grain(Structure):
    """
    We use the Structure class from pymatgen and add several new functions.
    """
    @wraps(Structure.__init__)
    def __init__(self, *args, **kwargs):
        super(Structure, self).__init__(*args, **kwargs)
        self._sites = list(self._sites)

    @staticmethod
    def get_b_from_a(grain_a):
        """
        Generate grain_b structure from a grain_a structure by rotating the sites in 
        grain_a by 180 degree.
        """
        for i in range(3):
            grain_b = grain_a.copy()
            anchor = grain_b.lattice.get_cartesian_coords(np.array([.0, .0, .0]))
            axis = [0, 0]
            axis.insert(i, 1)
            grain_b.rotate_sites(theta=np.radians(180), axis=axis, anchor=anchor)
            if grain_a != grain_b:
                break
        # grain_b = grain_a.copy()
        # csl_t = csl.transpose()
        # if sum(abs(csl_t[0]) - abs(csl_t[1])) > 0:
        #     axis = (0, 1, 0)
        # else:
        #     axis = (1, 0, 0)
        # anchor = grain_b.lattice.get_cartesian_coords(np.array([.0, .0, .0]))
        # # print(axis, anchor)
        # # exit()
        # grain_b.rotate_sites(theta=np.radians(180), axis=axis, anchor=anchor)
        return grain_b

    @classmethod
    def from_mp_id(cls, mp_id):
        """
        Get a structure from Materials Project database.

        Args:
            mp_id (str): Materials Project ID.

        Returns:
            A structure object.
        """
        from mp_api.client import MPRester

        mpr = MPRester()
        s = mpr.get_structure_by_material_id(mp_id, conventional_unit_cell=True)
        return cls.from_dict(s.as_dict())
    
    def get_orthogonal_matrix(self):
        warnings.warn("The lattice system of the grain is not orthogonal. "
                          "aimsgb will find a supercell of the grain structure "
                          "that is orthogonalized. This may take a while. ")
        cst = CubicSupercellTransformation(force_90_degrees=True,
                                           min_length=min(self.lattice.abc))
        _s = self.copy()
        _s = cst.apply_transformation(_s)
        matrix = [reduce_vector(i) for i in cst.transformation_matrix]
        _s = self.copy()
        _s.make_supercell(matrix)
        sm = StructureMatcher(attempt_supercell=True, primitive_cell=False)
        matrix = sm.get_supercell_matrix(_s, self)
        return np.array([reduce_vector(i) for i in matrix])
    
    def fix_sites_in_layers(self, layer_indices, tol=0.25, direction=2):
        """
        Fix sites in certain layers. The layer by layer direction is given by direction.
        This function is useful for selective dynamics calculations in VASP.

        Args:
            layer_indices (list): A list of layer indices.
            tol (float): Tolerance factor in Angstrom to determnine if sites are 
                in the same layer. Default to 0.25.
            direction (int): Direction to sort the sites by layers. 0: a, 1: b, 2: c
        """
        layers = self.sort_sites_in_layers(tol=tol, direction=direction)
        sd_sites = []
        for i, l in enumerate(layers):
            if i in layer_indices:
                sd_sites.extend(zip([[False, False, False]] * len(l), [_i[1] for _i in l]))
            else:
                sd_sites.extend(zip([[True, True, True]] * len(l), [_i[1] for _i in l]))
        values = [i[0] for i in sorted(sd_sites, key=lambda x: x[1])]
        self.add_site_property("selective_dynamics", values)

    def make_supercell(self, scaling_matrix):
        """
        Create a supercell. Very similar to pymatgen's Structure.make_supercell
        However, we need to make sure that all fractional coordinates that equal
        to 1 will become 0 and the lattice are redefined so that x_c = [0, 0, c]
    
        Args:
            scaling_matrix (3x3 matrix): The scaling matrix to make supercell.
        """
        s = self * scaling_matrix
        for i, site in enumerate(s):
            f_coords = np.mod(site.frac_coords, 1)
            # The following for loop is probably not necessary. But I will leave
            # it here for now.
            for j, v in enumerate(f_coords):
                if abs(v - 1) < 1e-6:
                    f_coords[j] = 0
            s[i] = PeriodicSite(site.specie, f_coords, site.lattice,
                                properties=site.properties)
        self._sites = s.sites
        self._lattice = s.lattice
        new_lat = Lattice.from_parameters(*s.lattice.parameters)
        self.lattice = new_lat

    def delete_bt_layer(self, bt, tol=0.25, direction=2):
        """
        Delete bottom or top layer of the structure.

        Args:
            bt (str): Specify whether it's a top or bottom layer delete. "b"
                means bottom layer and "t" means top layer.
            tol (float): Tolerance factor in Angstrom to determnine if sites are 
                in the same layer. Default to 0.25.
            direction (int): Direction to sort the sites by layers. 0: a, 1: b, 2: c
        """
        if bt == "t":
            l1, l2 = (-1, -2)
        else:
            l1, l2 = (0, 1)

        l = self.lattice.abc[direction]
        layers = self.sort_sites_in_layers(tol=tol, direction=direction)
        l_dist = abs(layers[l1][0][0].coords[direction] - layers[l2][0][0].coords[direction])
        l_vector = [1, 1]
        l_vector.insert(direction, (l - l_dist) / l)
        new_lat = Lattice(self.lattice.matrix * np.array(l_vector)[:, None])

        layers.pop(l1)
        sites = reduce(lambda x, y: np.concatenate((x, y), axis=0), layers)
        new_sites = []
        l_dist = 0 if bt == "t" else l_dist
        l_vector = [0, 0]
        l_vector.insert(direction, l_dist)
        for site, _ in sites:
            new_sites.append(PeriodicSite(site.specie, site.coords - l_vector,
                                          new_lat, coords_are_cartesian=True))
        self._sites = new_sites
        self._lattice = new_lat

    def sort_sites_in_layers(self, tol=0.25, direction=2):
        """
        Sort the sites in a structure by layers.

        Args:
            tol (float): Tolerance factor in Angstrom to determnine if sites are 
                in the same layer. Default to 0.25.
            direction (int): Direction to sort the sites by layers. 0: a, 1: b, 2: c

        Returns:
            Lists with a list of (site, index) in the same plane as one list.
        """
        sites_indices = sorted(zip(self.sites, range(len(self))), 
                               key=lambda x: x[0].frac_coords[direction])
        layers = []
        for k, g in groupby(sites_indices, key=lambda x: x[0].frac_coords[direction]):
            layers.append(list(g))
        new_layers = []
        k = -1
        for i in range(len(layers)):
            if i > k:
                tmp = layers[i]
                for j in range(i + 1, len(layers)):
                    if self.lattice.abc[direction] * abs(
                                    layers[j][0][0].frac_coords[direction] -
                                    layers[i][0][0].frac_coords[direction]) < tol:
                        tmp.extend(layers[j])
                        k = j
                    else:
                        break
                new_layers.append(sorted(tmp))
        # check if the 1st layer and last layer are actually the same layer
        # use the fractional as cartesian doesn't work for unorthonormal
        if self.lattice.abc[direction] * abs(
                                new_layers[0][0][0].frac_coords[direction] + 1 -
                                new_layers[-1][0][0].frac_coords[direction]) < tol:
            tmp = new_layers[0] + new_layers[-1]
            new_layers = new_layers[1:-1]
            new_layers.append(sorted(tmp))
        return new_layers

    # def set_orthogonal_grain(self):
    #     a, b, c = self.lattice.abc
    #     self.lattice = Lattice.orthorhombic(a, b, c)

    def set_orthogonal_grain(self, direction=2):
        """This method was adopted from pymatgen.surface.Slab.get_orthogonal_c_slab.
        **Note that this breaks inherent symmetries in the grain structure.**

        Args:
            direction (int): The lattice vector that is forced to be orthogonal. 0: a, 1: b, 2: c

        Returns:
            (Grain) An orthogonalized grain structure.
        """
        abc = list(self.lattice.matrix)
        _abc = copy.deepcopy(abc)
        _abc.pop(direction)
        new_c = np.cross(*_abc)
        new_c /= np.linalg.norm(new_c)
        new_c = np.dot(abc[direction], new_c) * new_c
        _abc.insert(direction, new_c)
        new_latt = Lattice(_abc)
        return Grain(
            lattice=new_latt, 
            species=self.species_and_occu, 
            coords=self.cart_coords,
            coords_are_cartesian=True,
            site_properties=self.site_properties
        )


    def build_grains(self, csl, direction, uc_a=1, uc_b=1):
        """
        Build structures for grain A and B from the coincidnet site lattice (CSL) matrix, 
        number of unit cell of grain A and number of unit cell of grain B. Each grain
        is essentially a supercell of the initial structure.

        Args:
            csl (3x3 matrix): CSL matrix (scaling matrix)
            direction (int): Stacking direction of GB. 0: a, 1: b, 2: c
            uc_a (int): Number of unit cell of grain A. Default to 1.
            uc_b (int): Number of unit cell of grain B. Default to 1.

        Returns:
            Grain objects for grain A and B.
        """
        csl_t = csl.transpose()
        # rotate along a longer axis between a and b
        grain_a = self.copy()
        grain_a.make_supercell(csl_t)
        # grain_a.to(filename='POSCAR_a')
        # exit()

        if not grain_a.lattice.is_orthogonal:
            warnings.warn("The lattice system of the grain is not orthogonal. The lattice "
                          "will be forced to be orthogonal. **Note that this breaks inherent symmetries of the grain.**")
            grain_a = grain_a.set_orthogonal_grain(direction)

        # grain_a.to(filename='POSCAR_a')
        # exit()
        temp_a = grain_a.copy()
        scale_vector = [1, 1]
        scale_vector.insert(direction, uc_b)
        temp_a.make_supercell(scale_vector)
        grain_b = self.get_b_from_a(temp_a)
        # make sure that all fractional coordinates that equal to 1 will become 0
        grain_b.make_supercell([1, 1, 1])

        scale_vector = [1, 1]
        scale_vector.insert(direction, uc_a)
        grain_a.make_supercell(scale_vector)

        # grain_b.to(filename='POSCAR_b')
        # exit()
        return grain_a, grain_b

    @classmethod
    def stack_grains(cls, grain_a, grain_b, vacuum=0.0, gap=0.0, direction=2,
                     delete_layer="0b0t0b0t", tol=0.25, to_primitive=True):
        """
        Build an interface structure by stacking two grains along a given direction.
        The grain_b a- and b-vectors will be forced to be the grain_a's
        a- and b-vectors.

        Args:
            grain_a (Grain): Substrate for the interface structure
            grain_b (Grain): Film for the interface structure
            vacuum (float): Vacuum space above the film in Angstroms. Default to 0.0
            gap (float): Gap between substrate and film in Angstroms. Default to 0.0
            direction (int): Stacking direction of the interface structure. 0: a, 1: b, 2: c.
            delete_layer (str): Delete top and bottom layers of the substrate and film.
                8 characters in total. The first 4 characters is for the substrate and
                the other 4 is for the film. "b" means bottom layer and "t" means
                top layer. Integer represents the number of layers to be deleted.
                Default to "0b0t0b0t", which means no deletion of layers. The
                direction of top and bottom layers is based on the given direction.
            tol (float): Tolerance factor in Angstrom to determnine if sites are 
                in the same layer. Default to 0.25.
            to_primitive (bool): Whether to get primitive structure of GB. Default to true.

        Returns:
             GB structure (Grain)
        """
        delete_layer = delete_layer.lower()
        delete = re.findall('(\d+)(\w)', delete_layer)
        if len(delete) != 4:
            raise ValueError(f"'{delete_layer}' is not supported. Please make sure the format "
                             "is 0b0t0b0t.")
        for i, v in enumerate(delete):
            for _ in range(int(v[0])):
                if i <= 1:
                    grain_a.delete_bt_layer(v[1], tol, direction)
                else:
                    grain_b.delete_bt_layer(v[1], tol, direction)
        abc_a = list(grain_a.lattice.abc)
        abc_b, angles = np.reshape(grain_b.lattice.parameters, (2, 3))
        if direction == 1:
            l = (abc_a[direction] + gap) * sin(radians(angles[2]))
        else:
            l = abc_a[direction] + gap
        abc_a[direction] += abc_b[direction] + 2 * gap + vacuum
        new_lat = Lattice.from_parameters(*abc_a, *angles)
        a_fcoords = new_lat.get_fractional_coords(grain_a.cart_coords)

        grain_a = Grain(new_lat, grain_a.species, a_fcoords, site_properties=grain_a.site_properties)
        l_vector = [0, 0]
        l_vector.insert(direction, l)
        b_fcoords = new_lat.get_fractional_coords(
            grain_b.cart_coords + l_vector)
        grain_b = Grain(new_lat, grain_b.species, b_fcoords, site_properties=grain_b.site_properties)

        structure = Grain.from_sites(grain_a[:] + grain_b[:])
        structure = structure.get_sorted_structure()
        if to_primitive:
            structure = structure.get_primitive_structure()

        return cls.from_dict(structure.as_dict())



================================================
FILE: aimsgb/grain_bound.py
================================================
from __future__ import division

import re
import warnings
import collections
from tabulate import tabulate
from fractions import Fraction
from math import sqrt, atan, degrees, pi
import numpy as np
from numpy import sin, cos, ceil, radians, inner, identity
from numpy.linalg import inv, det, norm, solve
from pymatgen.core.lattice import Lattice
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from aimsgb import Grain
from aimsgb.utils import reduce_vector, co_prime, plus_minus_gen, \
    is_integer, get_smallest_multiplier, reduce_integer, transpose_matrix

__author__ = "Jianli Cheng, Kesong Yang"
__copyright__ = "Copyright 2018, Yanggroup"
__maintainer__ = "Jianli Cheng"
__email__ = "jic198@ucsd.edu"
__status__ = "Production"
__date__ = "January 26, 2018"

# SIGMA_SYMBOL = u'\u03A3'
UNIMODULAR_MATRIX = np.array([identity(3),
                              [[1, 0, 1],
                               [0, 1, 0],
                               [0, 1, 1]],
                              [[1, 0, 1],
                               [0, 1, 0],
                               [0, 1, -1]],
                              [[1, 0, 1],
                               [0, 1, 0],
                               [-1, 1, 0]],
                              [[1, 0, 1],
                               [1, 1, 0],
                               [1, 1, 1]]])
STRUCTURE_MATRIX = np.array([identity(3),
                             [[0.5, -0.5, 0],
                              [0.5, 0.5, 0],
                              [0.5, 0.5, 1]],
                             [[0.5, 0.5, 0],
                              [0, 0.5, 0.5],
                              [0.5, 0, 0.5]]])
HEXAGONAL_001_CSL = {"5": {"012": [[1, -5, 0], [-2, -4, 0], [0, 0, -1]],
                           "013": [[1, -7, 0], [-3, -5, 0], [0, 0, -1]]},
                     "13": {"023": [[ 2, -8, 0], [-3, -7, 0], [0, 0, -1]], 
                            "015": [[1, -11, 0], [-5, -7, 0], [0, 0, -1]]},
                     "17": {"014": [[1, 7, 0], [-3, 5, 0], [0, 0, -1]], 
                            "035": [[2, -8, 0], [-3, -7, 0], [0, 0, -1]]},
                     }


@transpose_matrix
def reduce_csl(csl):
    """
    Reduce CSL matrix
    Args:
        csl: 3x3 matrix

    Returns:
        3x3 CSL matrix
    """
    csl = csl.round().astype(int)
    return np.array([reduce_vector(i) for i in csl])


@transpose_matrix
def o_lattice_to_csl(o_lattice, n):
    """
    The algorithm was borrowed from gosam project with slight changes.
    Link to the project: https://github.com/wojdyr/gosam

    There are two major steps: (1) Manipulate the columns of O-lattice to get
    an integral basis matrix for CSL: make two columns become integers and
    the remaining column can be multiplied by n whereby the determinant
    becomes sigma. (2) Simplify CSL so its vectors acquire the shortest length:
    decrease the integers in the matrix while keeping the determinant the same
    by adding other column vectors (row vectors in the following example) to a
    column vector. If after the addition or subtraction, the maximum value or
    absolute summation of added or subtracted vector is smaller than the
    original, then proceed the addition or subtraction.
    0 0 -1      0 0 -1      0 0 -1      0 0 -1      0 0 -1
    1 2 -1  ->  1 2 0   ->  1 2 0   ->  1 2 0   ->  1 2 0
    1 -3 2      1 -3 2      1 -3 1      1 -3 0      2 -1 0

    Args:
     o_lattice (3x3 array): O-lattice in crystal coordinates
     n (int): Number of O-lattice units per CSL unit

    Returns:
     CSL matrix (3x3 array) in crystal coordinates
    """
    csl = o_lattice.copy()
    if n < 0:
        csl[0] *= -1
        n *= -1
    while True:
        m = [get_smallest_multiplier(i) for i in csl]
        m_prod = np.prod(m)
        if m_prod <= n:
            for i in range(3):
                csl[i] *= m[i]
            if m_prod < n:
                assert n % m_prod == 0
                csl[0] *= n / m_prod
            break
        else:
            changed = False
            for i in range(3):
                for j in range(3):
                    if changed or i == j or m[i] == 1 or m[j] == 1:
                        continue
                    a, b = (i, j) if m[i] <= m[j] else (j, i)
                    for k in plus_minus_gen(1, m[b]):
                        handle = csl[a] + k * csl[b]
                        if get_smallest_multiplier(handle) < m[a]:
                            csl[a] += k * csl[b]
                            changed = True
                            break
            if not changed:
                # This situation rarely happens. Not sure if this solution is
                # legit, as det not equals to sigma. E.g. Sigma 115[113]
                for i in range(3):
                    csl[i] *= m[i]
                break
    csl = csl.round().astype(int)

    # Reshape CSL
    def simplify(l1, l2):
        x = abs(l1 + l2)
        y = abs(l1)
        changed = False
        while max(x) < max(y) or (max(x) == max(y) and sum(x) < sum(y)):
            l1 += l2
            changed = True
            x = abs(l1 + l2)
            y = abs(l1)
        return changed

    while True:
        changed = False
        for i in range(3):
            for j in range(3):
                if i != j and not changed:
                    changed = simplify(csl[i], csl[j]) or simplify(csl[i], -csl[j])
                    if changed:
                        break
        if not changed:
            break
    return csl


@transpose_matrix
def orthogonalize_csl(csl, axis):
    """
    (1) Set the 3rd column of csl same as the rotation axis. The algorithm was
    borrowed from gosam project with slight changes. Link to the project:
    https://github.com/wojdyr/gosam
    (2) Orthogonalize CSL, which is essentially a Gram-Schmidt process. At the
    same time, we want to make sure the column vectors of orthogonalized csl
    has the smallest value possible. That's why we compared two different ways.
    csl = [v1, v2, v3], vi is the column vector
    u1 = v3/||v3||, y2 = v1 - (v1 . u1)u1
    u2 = y2/||y2||, y3 = v3 - [(v3 . u1)u1 + (v3 . u2)u2]
    u3 = y3/||y3||
    """
    axis = np.array(axis)
    c = solve(csl.transpose(), axis)
    if not is_integer(c):
        mult = get_smallest_multiplier(c)
        c *= mult
    c = c.round().astype(int)
    ind = min([(i, v) for i, v in enumerate(c) if not np.allclose(v, 0)],
              key=lambda x: abs(x[1]))[0]
    if ind != 2:
        csl[ind], csl[2] = csl[2].copy(), -csl[ind]
        c[ind], c[2] = c[2], -c[ind]

    csl[2] = np.dot(c, csl)
    if c[2] < 0:
        csl[1] *= -1

    def get_integer(vec):
        # Used vec = np.array(vec, dtype=float) before, but does not work for
        # [5.00000000e-01, -5.00000000e-01,  2.22044605e-16] in Sigma3[112]
        vec = np.round(vec, 12)
        vec_sign = np.array([1 if abs(i) == i else -1 for i in vec])
        vec = list(abs(vec))
        new_vec = []
        if 0.0 in vec:
            zero_ind = vec.index(0)
            vec.pop(zero_ind)
            if 0.0 in vec:
                new_vec = [get_smallest_multiplier(vec) * i for i in vec]
            else:
                frac = Fraction(vec[0] / vec[1]).limit_denominator()
                new_vec = [frac.numerator, frac.denominator]
            new_vec.insert(zero_ind, 0)
        else:
            for i in range(len(vec) - 1):
                frac = Fraction(vec[i] / vec[i + 1]).limit_denominator()
                new_vec.extend([frac.numerator, frac.denominator])
            if new_vec[1] == new_vec[2]:
                new_vec = [new_vec[0], new_vec[1], new_vec[3]]
            else:
                new_vec = reduce_vector([new_vec[0] * new_vec[2],
                                         new_vec[1] * new_vec[2],
                                         new_vec[3] * new_vec[1]])
        assert is_integer(new_vec)
        return new_vec * vec_sign

    u1 = csl[2] / norm(csl[2])
    y2_1 = csl[1] - np.dot(csl[1], u1) * u1
    c0_1 = get_integer(y2_1)
    y2_2 = csl[0] - np.dot(csl[0], u1) * u1
    c0_2 = get_integer(y2_2)
    if sum(abs(c0_1)) > sum(abs(c0_2)):
        u2 = y2_2 / norm(y2_2)
        y3 = csl[1] - np.dot(csl[1], u1) * u1 - np.dot(csl[1], u2) * u2
        csl[1] = get_integer(y3)
        csl[0] = c0_2
    else:
        u2 = y2_1 / norm(y2_1)
        y3 = csl[0] - np.dot(csl[0], u1) * u1 - np.dot(csl[0], u2) * u2
        csl[1] = c0_1
        csl[0] = get_integer(y3)
    for i in range(3):
        for j in range(i + 1, 3):
            if not np.allclose(np.dot(csl[i], csl[j]), 0):
                raise ValueError(f"Non-orthogonal basis: {csl}")
    return csl.round().astype(int)


class GBInformation(dict):
    """
    GBInformation object essentially consists of a dictionary with information
    including sigma, CSL matrix, GB plane, rotation angle and rotation matrix
    """
    def __init__(self, axis, max_sigma, specific=False):
        """
        Creates a GBInformation object.
        Args:
            axis ([u, v, w]): Rotation axis.
            max_sigma (int): The largest sigma value
            specific (bool): Whether collecting information for a specific sigma
        """
        super(GBInformation, self).__init__()
        if max_sigma < 3:
            raise ValueError("Sigma should be larger than 2. '1' or '2' "
                             "means a layer by layer epitaxial film.")

        axis = np.array(reduce_vector(axis))
        self.axis = axis
        self.max_sigma = max_sigma
        self.specific = specific
        self.update(self.get_gb_info())

    def __str__(self):
        axis_str = "".join(map(str, self.axis))
        outs = [f"Grain boundary information for rotation axis: {axis_str}",
                f"Show the sigma values up to {self.max_sigma} (Note: * means twist GB, Theta is the rotation angle)"
                ]
        data = []
        to_s = lambda x: f"{x:.2f}"
        for key, item in sorted(self.items()):
            for i, value in enumerate(item["GB plane"]):
                count = -1
                for v in value:
                    count += 1
                    plane_str = f"({' '.join(map(str, v))})"
                    if v == list(self.axis):
                        plane_str += "*"
                    if count == 0:
                        row = [key, to_s(self[key]["Theta"][i])]
                    else:
                        row = [None, None]
                    csl = [" ".join('%2s' % k for k in j)
                           for j in self[key]["CSL matrix"][i]]
                    row.extend([plane_str, csl[count]])
                    data.append(row)
        outs.append(tabulate(data, numalign="center", tablefmt='orgtbl',
                             headers=["Sigma", "Theta", "GB plane", "CSL matrix"]))
        return "\n".join(outs)

    def get_gb_info(self):
        """
        Calculate sigma, rotation angle, GB plane, rotation matrix and CSL matrix
        The algorithm for getting sigma, m, n, theta is from H. Grimmer:
        https://doi.org/10.1107/S0108767384000246

        Returns:
            gb_info(dict)
        """
        max_m = int(ceil(sqrt(4 * self.max_sigma)))
        gb_info = collections.defaultdict(dict)
        sigma_theta = collections.defaultdict(list)

        for m in range(max_m):
            for n in range(max_m):
                if not co_prime(m, n):
                    continue
                sigma = self.get_sigma(m, n)
                if self.specific and sigma != self.max_sigma:
                    continue
                if not sigma or sigma > self.max_sigma:
                    continue
                theta = self.get_theta(m, n)
                sigma_theta[sigma].append([theta, m, n])

        if not sigma_theta:
            raise ValueError("Cannot find any matching GB. Most likely there "
                             f"is no sigma {self.max_sigma}{self.axis} GB.")
        for sigma in sigma_theta:
            sigma_theta[sigma] = sorted(sigma_theta[sigma], key=lambda t: t[0])
            min_theta = sigma_theta[sigma][0][0]
            rot_matrix = self.get_rotate_matrix(min_theta)
            csl_matrix = self.get_csl_matrix(sigma, rot_matrix)
            csl = orthogonalize_csl(csl_matrix, self.axis)
            # Sometime when getting CSL from O-lattice, det not equals to sigma.
            # That's why it needs to be reduced. E.g. Sigma 115[113]
            csl = reduce_csl(csl)
            all_csl = [csl]
            if sorted(self.axis) == [0, 0, 1]:
                ind = 0
                while True:
                    m, n = sigma_theta[sigma][ind][1:]
                    ext_csl = csl.copy()
                    for i, v in enumerate(ext_csl):
                        if sorted(v) != [0, 0, 1]:
                            if abs(v[0]) > abs(v[1]):
                                ext_csl[i] = [v[0] / abs(v[0]) * m, v[1] / abs(v[1]) * n, 0]
                            else:
                                ext_csl[i] = [v[0] / abs(v[0]) * n, v[1] / abs(v[1]) * m, 0]
                    if (csl == ext_csl).all():
                        ind += 1
                    else:
                        break
                all_csl = [csl, ext_csl]
                if ind:
                    gb_info[sigma] = {"Theta": [min_theta, 90 - min_theta]}
                else:
                    gb_info[sigma] = {"Theta": [90 - min_theta, min_theta]}
            else:
                gb_info[sigma] = {"Theta": [min_theta]}

            gb_info[sigma].update({"GB plane": [[list(j) for j in i.transpose()]
                                             for i in all_csl],
                                   "Rotation matrix": rot_matrix, "CSL matrix": all_csl})
        return gb_info

    def get_theta(self, m, n):
        """
        Calculate rotation angle from m, n

        Args:
            m (int)
            n (int)

        Returns:
            theta
        """
        try:
            return degrees(2 * atan(n * sqrt(inner(self.axis, self.axis)) / m))
        except ZeroDivisionError:
            return degrees(pi)

    def get_sigma(self, m, n):
        """
        Calculate sigma from m, n

        Args:
            m (int)
            n (int)

        Returns:
            sigma
        """
        sigma = m ** 2 + n ** 2 * inner(self.axis, self.axis)
        alph = 1
        while sigma != 0 and sigma % 2 == 0:
            alph *= 2
            sigma //= 2
        return sigma if sigma > 1 else None

    def get_rotate_matrix(self, angle):
        """
        use Rodrigues' rotation formula to get rotation matrix
        https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
        """
        rotate_axis = np.array(self.axis / sqrt(inner(self.axis, self.axis)))
        angle = radians(angle)
        k_matrix = np.array([[0, -rotate_axis[2], rotate_axis[1]],
                            [rotate_axis[2], 0, -rotate_axis[0]],
                            [-rotate_axis[1], rotate_axis[0], 0]])
        return identity(3) + k_matrix * sin(angle) + \
               np.dot(k_matrix, k_matrix) * (1 - cos(angle))

    def get_csl_matrix(self, sigma, rotate_matrix):
        """
        Calculate CSL matrix from sigma and rotation matrix. The algorithm is
        from H. Grimmer et al. https://doi.org/10.1107/S056773947400043X
        For the structure matrix, we will use the identity matrix. Since for
        the initial structure that is not primitive cubic, we will transform it
        to conventional standard cell.
        Args:
            sigma (int): Degree of fit
            rotate_matrix (3x3 matrix): Rotation matrix

        Returns:
            3x3 CSL matrix
        """
        s = STRUCTURE_MATRIX[0]
        for u in UNIMODULAR_MATRIX:
            t = np.eye(3) - np.dot(np.dot(np.dot(u, inv(s)), inv(rotate_matrix)), s)
            if abs(det(t)) > 1e-6:
                break
        o_lattice = np.round(inv(t), 12)
        n = np.round(sigma / det(o_lattice), 6)
        csl_matrix = o_lattice_to_csl(o_lattice, n)
        return csl_matrix


class GrainBoundary(object):
    """
    A grain boundary (GB) object. The initial structure can be cubic or non-cubic
    crystal. If non-cubic, the crystal will be transferred to conventional cell.
    The generated GB could be either tilted or twisted based on the given GB
    plane. If the GB plane is parallel to the rotation axis, the generated GB
    will be a twisted one. Otherwise, tilted.
    """
    def __init__(self, axis, sigma, plane, initial_struct, uc_a=1, uc_b=1):
        """
        Build grain boundary based on rotation axis, sigma, GB plane, grain size,
        if_model and vacuum thickness.

        Args:
            axis ([u, v, w]): Rotation axis.
            sigma (int): The area ratio between the unit cell of CSL and the
                given crystal lattice.
            plane ([h, k, l]): Miller index of GB plane. If the GB plane is parallel
                to the rotation axis, the generated GB will be a twist GB. If they
                are perpendicular, the generated GB will be a tilt GB.
            initial_struct (Grain): Initial input structure. Must be an
                object of Grain
            uc_a (int): Number of unit cell of grain A. Default to 1.
            uc_b (int): Number of unit cell of grain B. Default to 1.
        """
        if not isinstance(initial_struct, Grain):
            raise ValueError("The input 'initial_struct' must be an object "
                             "of Grain.")
        self.axis = axis
        self.plane = list(plane)
        self.plane_str = " ".join(map(str, self.plane))

        if sigma % 2 == 0:
            reduce_sigma = reduce_integer(sigma)
            warnings.warn(
                f"{sigma} is an even number. However sigma must be an odd number. "
                f"We will choose sigma={reduce_sigma}.",
                RuntimeWarning)
            sigma = reduce_sigma
        self.sigma = sigma
        self.gb_info = GBInformation(self.axis, self.sigma, specific=True)
        self.csl = None
        for i, v in enumerate(self.gb_info[self.sigma]["GB plane"]):
            if self.plane in v:
                self.csl = self.gb_info[self.sigma]["CSL matrix"][i]
        if self.csl is None:
            avail_plane = "', '".join(["', '".join([" ".join(map(str, j)) for j in i])
                                    for i in self.gb_info[self.sigma]["GB plane"]])
            raise ValueError(f"The given GB plane '{self.plane_str}' cannot be realized. Choose "
                             f"the plane in ['{avail_plane}']")
        self.direction = None
        for i, v in enumerate(self.csl.transpose()):
            if self.plane == list(v):
                self.direction = i
        new_s = SpacegroupAnalyzer(initial_struct).get_conventional_standard_structure()
        initial_struct = Grain.from_sites(new_s[:])

        if SpacegroupAnalyzer(initial_struct).get_lattice_type() in ['hexagonal', 'rhombohedral']:
            axis_str = "".join(map(str, sorted(self.axis)))
            if axis_str == '001':
                _plane_str = "".join(map(str, sorted(map(abs, self.plane))))
                if HEXAGONAL_001_CSL.get(str(self.sigma), {}).get(_plane_str):
                    self.csl = np.array(HEXAGONAL_001_CSL.get(str(self.sigma), {}).get(_plane_str))

        self._grain_a, self._grain_b = initial_struct.build_grains(
            self.csl, self.direction, uc_a, uc_b)

    @property
    def rot_matrix(self):
        """
        Rotation matrix for calculating CSL matrix
        """
        return self.gb_info[self.sigma]["Rotation matrix"]

    @property
    def theta(self):
        """
        Rotation angle for calculating rotation matrix and to bring two grains
        into a perfect match
        """
        return self.gb_info[self.sigma]["Theta"]

    # @property
    # def csl(self):
    #     """
    #     CSL matrix
    #     """
    #     for i, v in enumerate(self.gb_info[self.sigma]["GB plane"]):
    #         if self.plane in v:
    #             return self.gb_info[self.sigma]["CSL matrix"][i]
    #     avail_plane = ", ".join([", ".join([" ".join(map(str, j)) for j in i])
    #                              for i in self.gb_info[self.sigma]["GB plane"]])
    #     raise ValueError(f"The given GB plane '{self.plane_str}' cannot be realized. Choose "
    #                      f"the plane in [{avail_plane}]")

    @property
    def grain_a(self):
        """
        Grain class instance for grain A
        """
        return self._grain_a

    @property
    def grain_b(self):
        """
        Grain class instance for grain B
        """
        return self._grain_b

    def build_gb(self, vacuum=0.0, add_if_dist=0.0, to_primitive=True,
                 delete_layer="0b0t0b0t", tol=0.25):
        """
        Build the GB based on the given crystal, uc of grain A and B, if_model,
        vacuum thickness, distance between two grains and tolerance factor.
        Args:
            vacuum (float), Angstrom: Vacuum thickness for GB.
                Default to 0.0
            add_if_dist (float), Angstrom: Add extra distance at the interface
                between two grains.
                Default to 0.0
            to_primitive (bool): Whether to get primitive structure of GB.
                Default to true.
            delete_layer (str): Delete interface layers on both sides of each grain.
                8 characters in total. The first 4 characters is for grain A and
                the other 4 is for grain B. "b" means bottom layer and "t" means
                top layer. Integer represents the number of layers to be deleted.
                Default to "0b0t0b0t", which means no deletion of layers. The
                direction of top and bottom layers is based on direction.
            tol (float), Angstrom: Tolerance factor to determine whether two
                atoms are at the same plane.
                Default to 0.25
        Returns:
             GB structure (Grain)
        """
        warnings.warn("The build_gb method is deprecated. Please use Grain.stack_grains methode instead.")
        
        ind = self.direction
        delete_layer = delete_layer.lower()
        delete = re.findall('(\d+)(\w)', delete_layer)
        if len(delete) != 4:
            raise ValueError(f"'{delete_layer}' is not supported. Please make sure the format "
                             "is 0b0t0b0t.")
        for i, v in enumerate(delete):
            for j in range(int(v[0])):
                if i <= 1:
                    self.grain_a.delete_bt_layer(v[1], tol, ind)
                else:
                    self.grain_b.delete_bt_layer(v[1], tol, ind)
        abc_a = list(self.grain_a.lattice.abc)
        abc_b, angles = np.reshape(self.grain_b.lattice.parameters, (2, 3))
        if ind == 1:
            l = (abc_a[ind] + add_if_dist) * sin(radians(angles[2]))
        else:
            l = abc_a[ind] + add_if_dist
        abc_a[ind] += abc_b[ind] + 2 * add_if_dist + vacuum
        new_lat = Lattice.from_parameters(*abc_a, *angles)
        a_fcoords = new_lat.get_fractional_coords(self.grain_a.cart_coords)

        grain_a = Grain(new_lat, self.grain_a.species, a_fcoords)
        l_vector = [0, 0]
        l_vector.insert(ind, l)
        b_fcoords = new_lat.get_fractional_coords(
            self.grain_b.cart_coords + l_vector)
        grain_b = Grain(new_lat, self.grain_b.species, b_fcoords)

        gb = Grain.from_sites(grain_a[:] + grain_b[:])
        gb = gb.get_sorted_structure()
        if to_primitive:
            gb = gb.get_primitive_structure()
        return gb



================================================
FILE: aimsgb/utils.py
================================================
from __future__ import division

import numpy as np
from numpy.linalg import norm
from functools import wraps, reduce
try:
    from math import gcd as pygcd
except ImportError:
    from fractions import gcd as pygcd

__author__ = "Jianli Cheng, Kesong Yang"
__copyright__ = "Copyright 2018, Yanggroup"
__maintainer__ = "Jianli Cheng"
__email__ = "jic198@ucsd.edu"
__status__ = "Production"
__date__ = "January 26, 2018"


def reduce_vector(vector):
    """
    Reduce a vector
    """
    d = abs(reduce(gcd, vector))
    vector = tuple([int(i / d) for i in vector])

    return vector


def co_prime(a, b):
    """
    Check if two integers are co-prime
    """
    return gcd(a, b) in (0, 1)


def plus_minus_gen(start, end):
    """
    Generate a list of plus and minus alternating integers
    """
    for i in range(start, end):
        yield i
        yield -i


def is_integer(a, tol=1e-5):
    """
    Check whether the number is integer
    """
    return norm(abs(a - np.round(a))) < tol


def get_smallest_multiplier(a, max_n=10000):
    """
    Get the smallest multiplier to make the list with all integers
    Args:
        a (list): A list of numbers
        max_n (int): The up limit to search multiplier

    Returns:
        The smallest integer multiplier
    """
    a = np.array(a)
    for i in range(1, max_n):
        if is_integer(i * a):
            return i
    raise ValueError("Cannot find an integer matrix with multiplier "
                     "searched already up to %s" % max_n)


def reduce_integer(integer):
    """
    Get the odd number for an integer
    """
    while gcd(integer, 2) != 1:
        integer //= 2
    return integer


def gcd(*numbers):
    """
    Get a greatest common divisor for a list of numbers
    """
    n = numbers[0]
    for i in numbers:
        n = pygcd(n, i)
    return n


def transpose_matrix(func):
    """
    Transpose the first argument and the return value.
    Args:
        func: The function that uses transpose_matrix as a decorator.
    """

    @wraps(func)
    def transpose(*args, **kwargs):
        args_list = list(args)
        args_list[0] = args_list[0].transpose()
        matrix = func(*args_list, **kwargs)
        return matrix.transpose()

    return transpose





================================================
FILE: .github/workflows/testing.yml
================================================
name: Testing

on: [push, pull_request]
  
jobs:
  build:
    strategy:
      matrix:
        os: [ubuntu-latest]
        python-version: ["3.9", "3.10"]

    runs-on: ${{ matrix.os }}

    steps:
    - uses: actions/checkout@v3
    - name: Set up Python ${{ matrix.python-version }}
      uses: actions/setup-python@v4
      with:
        python-version: ${{ matrix.python-version }}
        cache: 'pip'
        cache-dependency-path: '**/setup.py'
    - name: Install dependencies
      run: |
        pip install -e .
    - name: pytest
      run: |
        pytest

  release:
    runs-on: ubuntu-latest
    needs: build
    if: github.event_name == 'release' && needs.build.result == 'success'
    permissions:
      # For pypi trusted publishing
      id-token: write
    steps:
      - name: Check out repo
        uses: actions/checkout@v3

      - name: Build and upload dist
        run: |
          pip install build
          python -m build

      - name: Publish to PyPi
        uses: pypa/gh-action-pypi-publish@release/v1
        with:
          skip-existing: true
          verbose: true        
