xref: /dports/science/py-pymatgen/pymatgen-2022.0.15/pymatgen/symmetry/kpath.py

# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.

"""
Provides classes for generating high-symmetry k-paths using different conventions.
"""

import abc
import itertools
import operator
from math import ceil, cos, e, pi, sin, tan
from warnings import warn

import networkx as nx
import numpy as np
import spglib
from monty.dev import requires
from scipy.linalg import sqrtm

from pymatgen.core.lattice import Lattice
from pymatgen.core.operations import MagSymmOp, SymmOp
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer

try:
    from seekpath import get_path  # type: ignore
except ImportError:
    get_path = None

__author__ = "Geoffroy Hautier, Katherine Latimer, Jason Munro"
__copyright__ = "Copyright 2020, The Materials Project"
__version__ = "0.1"
__maintainer__ = "Jason Munro"
__email__ = "jmunro@lbl.gov"
__status__ = "Development"
__date__ = "March 2020"


class KPathBase(metaclass=abc.ABCMeta):
    """
    This is the base class for classes used to generate high-symmetry
    paths in reciprocal space (k-paths) for band structure calculations.
    """

    @abc.abstractmethod
    def __init__(self, structure, symprec=0.01, angle_tolerance=5, atol=1e-5, *args, **kwargs):
        """
        Args:
        structure (Structure): Structure object
        symprec (float): Tolerance for symmetry finding
        angle_tolerance (float): Angle tolerance for symmetry finding.
        atol (float): Absolute tolerance used to compare structures
            and determine symmetric equivalence of points and lines
            in the BZ.
        """

        self._structure = structure
        self._latt = self._structure.lattice
        self._rec_lattice = self._structure.lattice.reciprocal_lattice
        self._kpath = None

        self._symprec = symprec
        self._atol = atol
        self._angle_tolerance = angle_tolerance

    @property
    def structure(self):
        """
        Returns:
        The input structure
        """
        return self._structure

    @property
    def lattice(self):
        """
        Returns:
        The real space lattice
        """
        return self._latt

    @property
    def rec_lattice(self):
        """
        Returns:
        The reciprocal space lattice
        """
        return self._rec_lattice

    @property
    def kpath(self):
        """
        Returns:
        The symmetry line path in reciprocal space
        """
        return self._kpath

    def get_kpoints(self, line_density=20, coords_are_cartesian=True):
        """
        Returns:
        the kpoints along the paths in cartesian coordinates
        together with the labels for symmetry points -Wei.
        """
        list_k_points = []
        sym_point_labels = []
        for b in self.kpath["path"]:
            for i in range(1, len(b)):
                start = np.array(self.kpath["kpoints"][b[i - 1]])
                end = np.array(self.kpath["kpoints"][b[i]])
                distance = np.linalg.norm(
                    self._rec_lattice.get_cartesian_coords(start) - self._rec_lattice.get_cartesian_coords(end)
                )
                nb = int(ceil(distance * line_density))
                if nb == 0:
                    continue
                sym_point_labels.extend([b[i - 1]] + [""] * (nb - 1) + [b[i]])
                list_k_points.extend(
                    [
                        self._rec_lattice.get_cartesian_coords(start)
                        + float(i)
                        / float(nb)
                        * (self._rec_lattice.get_cartesian_coords(end) - self._rec_lattice.get_cartesian_coords(start))
                        for i in range(0, nb + 1)
                    ]
                )
        if coords_are_cartesian:
            return list_k_points, sym_point_labels

        frac_k_points = [self._rec_lattice.get_fractional_coords(k) for k in list_k_points]
        return frac_k_points, sym_point_labels


class KPathSetyawanCurtarolo(KPathBase):
    """
    This class looks for path along high symmetry lines in
    the Brillouin Zone.
    It is based on Setyawan, W., & Curtarolo, S. (2010).
    High-throughput electronic band structure calculations:
    Challenges and tools. Computational Materials Science,
    49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010
    It should be used with primitive structures that
    comply with the definition from the paper.
    The symmetry is determined by spglib through the
    SpacegroupAnalyzer class. The analyzer can be used to
    produce the correct primitive structure (method
    get_primitive_standard_structure(international_monoclinic=False)).
    A warning will signal possible compatibility problems
    with the given structure. KPoints from get_kpoints() method
    are returned in the reciprocal cell basis defined in the paper.
    """

    def __init__(self, structure, symprec=0.01, angle_tolerance=5, atol=1e-5):
        """
        Args:
        structure (Structure): Structure object
        symprec (float): Tolerance for symmetry finding
        angle_tolerance (float): Angle tolerance for symmetry finding.
        atol (float): Absolute tolerance used to compare the input
            structure with the one expected as primitive standard.
            A warning will be issued if the lattices don't match.
        """
        if "magmom" in structure.site_properties.keys():
            warn(
                "'magmom' entry found in site properties but will be ignored \
                  for the Setyawan and Curtarolo convention."
            )

        super().__init__(structure, symprec=symprec, angle_tolerance=angle_tolerance, atol=atol)

        self._sym = SpacegroupAnalyzer(structure, symprec=symprec, angle_tolerance=angle_tolerance)
        self._prim = self._sym.get_primitive_standard_structure(international_monoclinic=False)
        self._conv = self._sym.get_conventional_standard_structure(international_monoclinic=False)
        self._rec_lattice = self._prim.lattice.reciprocal_lattice

        # Note: this warning will be issued for space groups 38-41, since the primitive cell must be
        # reformatted to match Setyawan/Curtarolo convention in order to work with the current k-path
        # generation scheme.
        if not np.allclose(self._structure.lattice.matrix, self._prim.lattice.matrix, atol=atol):
            warn(
                "The input structure does not match the expected standard primitive! "
                "The path can be incorrect. Use at your own risk."
            )

        lattice_type = self._sym.get_lattice_type()
        spg_symbol = self._sym.get_space_group_symbol()

        if lattice_type == "cubic":
            if "P" in spg_symbol:
                self._kpath = self.cubic()
            elif "F" in spg_symbol:
                self._kpath = self.fcc()
            elif "I" in spg_symbol:
                self._kpath = self.bcc()
            else:
                warn("Unexpected value for spg_symbol: %s" % spg_symbol)

        elif lattice_type == "tetragonal":
            if "P" in spg_symbol:
                self._kpath = self.tet()
            elif "I" in spg_symbol:
                a = self._conv.lattice.abc[0]
                c = self._conv.lattice.abc[2]
                if c < a:
                    self._kpath = self.bctet1(c, a)
                else:
                    self._kpath = self.bctet2(c, a)
            else:
                warn("Unexpected value for spg_symbol: %s" % spg_symbol)

        elif lattice_type == "orthorhombic":
            a = self._conv.lattice.abc[0]
            b = self._conv.lattice.abc[1]
            c = self._conv.lattice.abc[2]

            if "P" in spg_symbol:
                self._kpath = self.orc()

            elif "F" in spg_symbol:
                if 1 / a ** 2 > 1 / b ** 2 + 1 / c ** 2:
                    self._kpath = self.orcf1(a, b, c)
                elif 1 / a ** 2 < 1 / b ** 2 + 1 / c ** 2:
                    self._kpath = self.orcf2(a, b, c)
                else:
                    self._kpath = self.orcf3(a, b, c)

            elif "I" in spg_symbol:
                self._kpath = self.orci(a, b, c)

            elif "C" in spg_symbol or "A" in spg_symbol:
                self._kpath = self.orcc(a, b, c)
            else:
                warn("Unexpected value for spg_symbol: %s" % spg_symbol)

        elif lattice_type == "hexagonal":
            self._kpath = self.hex()

        elif lattice_type == "rhombohedral":
            alpha = self._prim.lattice.parameters[3]
            if alpha < 90:
                self._kpath = self.rhl1(alpha * pi / 180)
            else:
                self._kpath = self.rhl2(alpha * pi / 180)

        elif lattice_type == "monoclinic":
            a, b, c = self._conv.lattice.abc
            alpha = self._conv.lattice.parameters[3]
            # beta = self._conv.lattice.parameters[4]

            if "P" in spg_symbol:
                self._kpath = self.mcl(b, c, alpha * pi / 180)

            elif "C" in spg_symbol:
                kgamma = self._rec_lattice.parameters[5]
                if kgamma > 90:
                    self._kpath = self.mclc1(a, b, c, alpha * pi / 180)
                if kgamma == 90:
                    self._kpath = self.mclc2(a, b, c, alpha * pi / 180)
                if kgamma < 90:
                    if b * cos(alpha * pi / 180) / c + b ** 2 * sin(alpha * pi / 180) ** 2 / a ** 2 < 1:
                        self._kpath = self.mclc3(a, b, c, alpha * pi / 180)
                    if b * cos(alpha * pi / 180) / c + b ** 2 * sin(alpha * pi / 180) ** 2 / a ** 2 == 1:
                        self._kpath = self.mclc4(a, b, c, alpha * pi / 180)
                    if b * cos(alpha * pi / 180) / c + b ** 2 * sin(alpha * pi / 180) ** 2 / a ** 2 > 1:
                        self._kpath = self.mclc5(a, b, c, alpha * pi / 180)
            else:
                warn("Unexpected value for spg_symbol: %s" % spg_symbol)

        elif lattice_type == "triclinic":
            kalpha = self._rec_lattice.parameters[3]
            kbeta = self._rec_lattice.parameters[4]
            kgamma = self._rec_lattice.parameters[5]
            if kalpha > 90 and kbeta > 90 and kgamma > 90:
                self._kpath = self.tria()
            if kalpha < 90 and kbeta < 90 and kgamma < 90:
                self._kpath = self.trib()
            if kalpha > 90 and kbeta > 90 and kgamma == 90:
                self._kpath = self.tria()
            if kalpha < 90 and kbeta < 90 and kgamma == 90:
                self._kpath = self.trib()

        else:
            warn("Unknown lattice type %s" % lattice_type)

    @property
    def conventional(self):
        """
        Returns:
            The conventional cell structure
        """
        return self._conv

    @property
    def prim(self):
        """
        Returns:
            The primitive cell structure
        """
        return self._prim

    @property
    def prim_rec(self):
        """
        Returns:
            The primitive reciprocal cell structure
        """
        return self._rec_lattice

    def cubic(self):
        """
        CUB Path
        """
        self.name = "CUB"
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "X": np.array([0.0, 0.5, 0.0]),
            "R": np.array([0.5, 0.5, 0.5]),
            "M": np.array([0.5, 0.5, 0.0]),
        }
        path = [["\\Gamma", "X", "M", "\\Gamma", "R", "X"], ["M", "R"]]
        return {"kpoints": kpoints, "path": path}

    def fcc(self):
        """
        FCC Path
        """
        self.name = "FCC"
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "K": np.array([3.0 / 8.0, 3.0 / 8.0, 3.0 / 4.0]),
            "L": np.array([0.5, 0.5, 0.5]),
            "U": np.array([5.0 / 8.0, 1.0 / 4.0, 5.0 / 8.0]),
            "W": np.array([0.5, 1.0 / 4.0, 3.0 / 4.0]),
            "X": np.array([0.5, 0.0, 0.5]),
        }
        path = [
            ["\\Gamma", "X", "W", "K", "\\Gamma", "L", "U", "W", "L", "K"],
            ["U", "X"],
        ]
        return {"kpoints": kpoints, "path": path}

    def bcc(self):
        """
        BCC Path
        """
        self.name = "BCC"
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "H": np.array([0.5, -0.5, 0.5]),
            "P": np.array([0.25, 0.25, 0.25]),
            "N": np.array([0.0, 0.0, 0.5]),
        }
        path = [["\\Gamma", "H", "N", "\\Gamma", "P", "H"], ["P", "N"]]
        return {"kpoints": kpoints, "path": path}

    def tet(self):
        """
        TET Path
        """
        self.name = "TET"
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "A": np.array([0.5, 0.5, 0.5]),
            "M": np.array([0.5, 0.5, 0.0]),
            "R": np.array([0.0, 0.5, 0.5]),
            "X": np.array([0.0, 0.5, 0.0]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            ["\\Gamma", "X", "M", "\\Gamma", "Z", "R", "A", "Z"],
            ["X", "R"],
            ["M", "A"],
        ]
        return {"kpoints": kpoints, "path": path}

    def bctet1(self, c, a):
        """
        BCT1 Path
        """
        self.name = "BCT1"
        eta = (1 + c ** 2 / a ** 2) / 4.0
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "M": np.array([-0.5, 0.5, 0.5]),
            "N": np.array([0.0, 0.5, 0.0]),
            "P": np.array([0.25, 0.25, 0.25]),
            "X": np.array([0.0, 0.0, 0.5]),
            "Z": np.array([eta, eta, -eta]),
            "Z_1": np.array([-eta, 1 - eta, eta]),
        }
        path = [["\\Gamma", "X", "M", "\\Gamma", "Z", "P", "N", "Z_1", "M"], ["X", "P"]]
        return {"kpoints": kpoints, "path": path}

    def bctet2(self, c, a):
        """
        BCT2 Path
        """
        self.name = "BCT2"
        eta = (1 + a ** 2 / c ** 2) / 4.0
        zeta = a ** 2 / (2 * c ** 2)
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "N": np.array([0.0, 0.5, 0.0]),
            "P": np.array([0.25, 0.25, 0.25]),
            "\\Sigma": np.array([-eta, eta, eta]),
            "\\Sigma_1": np.array([eta, 1 - eta, -eta]),
            "X": np.array([0.0, 0.0, 0.5]),
            "Y": np.array([-zeta, zeta, 0.5]),
            "Y_1": np.array([0.5, 0.5, -zeta]),
            "Z": np.array([0.5, 0.5, -0.5]),
        }
        path = [
            [
                "\\Gamma",
                "X",
                "Y",
                "\\Sigma",
                "\\Gamma",
                "Z",
                "\\Sigma_1",
                "N",
                "P",
                "Y_1",
                "Z",
            ],
            ["X", "P"],
        ]
        return {"kpoints": kpoints, "path": path}

    def orc(self):
        """
        ORC Path
        """
        self.name = "ORC"
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "R": np.array([0.5, 0.5, 0.5]),
            "S": np.array([0.5, 0.5, 0.0]),
            "T": np.array([0.0, 0.5, 0.5]),
            "U": np.array([0.5, 0.0, 0.5]),
            "X": np.array([0.5, 0.0, 0.0]),
            "Y": np.array([0.0, 0.5, 0.0]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            ["\\Gamma", "X", "S", "Y", "\\Gamma", "Z", "U", "R", "T", "Z"],
            ["Y", "T"],
            ["U", "X"],
            ["S", "R"],
        ]
        return {"kpoints": kpoints, "path": path}

    def orcf1(self, a, b, c):
        """
        ORFC1 Path
        """
        self.name = "ORCF1"
        zeta = (1 + a ** 2 / b ** 2 - a ** 2 / c ** 2) / 4
        eta = (1 + a ** 2 / b ** 2 + a ** 2 / c ** 2) / 4

        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "A": np.array([0.5, 0.5 + zeta, zeta]),
            "A_1": np.array([0.5, 0.5 - zeta, 1 - zeta]),
            "L": np.array([0.5, 0.5, 0.5]),
            "T": np.array([1, 0.5, 0.5]),
            "X": np.array([0.0, eta, eta]),
            "X_1": np.array([1, 1 - eta, 1 - eta]),
            "Y": np.array([0.5, 0.0, 0.5]),
            "Z": np.array([0.5, 0.5, 0.0]),
        }
        path = [
            ["\\Gamma", "Y", "T", "Z", "\\Gamma", "X", "A_1", "Y"],
            ["T", "X_1"],
            ["X", "A", "Z"],
            ["L", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}

    def orcf2(self, a, b, c):
        """
        ORFC2 Path
        """
        self.name = "ORCF2"
        phi = (1 + c ** 2 / b ** 2 - c ** 2 / a ** 2) / 4
        eta = (1 + a ** 2 / b ** 2 - a ** 2 / c ** 2) / 4
        delta = (1 + b ** 2 / a ** 2 - b ** 2 / c ** 2) / 4
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "C": np.array([0.5, 0.5 - eta, 1 - eta]),
            "C_1": np.array([0.5, 0.5 + eta, eta]),
            "D": np.array([0.5 - delta, 0.5, 1 - delta]),
            "D_1": np.array([0.5 + delta, 0.5, delta]),
            "L": np.array([0.5, 0.5, 0.5]),
            "H": np.array([1 - phi, 0.5 - phi, 0.5]),
            "H_1": np.array([phi, 0.5 + phi, 0.5]),
            "X": np.array([0.0, 0.5, 0.5]),
            "Y": np.array([0.5, 0.0, 0.5]),
            "Z": np.array([0.5, 0.5, 0.0]),
        }
        path = [
            ["\\Gamma", "Y", "C", "D", "X", "\\Gamma", "Z", "D_1", "H", "C"],
            ["C_1", "Z"],
            ["X", "H_1"],
            ["H", "Y"],
            ["L", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}

    def orcf3(self, a, b, c):
        """
        ORFC3 Path
        """
        self.name = "ORCF3"
        zeta = (1 + a ** 2 / b ** 2 - a ** 2 / c ** 2) / 4
        eta = (1 + a ** 2 / b ** 2 + a ** 2 / c ** 2) / 4
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "A": np.array([0.5, 0.5 + zeta, zeta]),
            "A_1": np.array([0.5, 0.5 - zeta, 1 - zeta]),
            "L": np.array([0.5, 0.5, 0.5]),
            "T": np.array([1, 0.5, 0.5]),
            "X": np.array([0.0, eta, eta]),
            "X_1": np.array([1, 1 - eta, 1 - eta]),
            "Y": np.array([0.5, 0.0, 0.5]),
            "Z": np.array([0.5, 0.5, 0.0]),
        }
        path = [
            ["\\Gamma", "Y", "T", "Z", "\\Gamma", "X", "A_1", "Y"],
            ["X", "A", "Z"],
            ["L", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}

    def orci(self, a, b, c):
        """
        ORCI Path
        """
        self.name = "ORCI"
        zeta = (1 + a ** 2 / c ** 2) / 4
        eta = (1 + b ** 2 / c ** 2) / 4
        delta = (b ** 2 - a ** 2) / (4 * c ** 2)
        mu = (a ** 2 + b ** 2) / (4 * c ** 2)
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "L": np.array([-mu, mu, 0.5 - delta]),
            "L_1": np.array([mu, -mu, 0.5 + delta]),
            "L_2": np.array([0.5 - delta, 0.5 + delta, -mu]),
            "R": np.array([0.0, 0.5, 0.0]),
            "S": np.array([0.5, 0.0, 0.0]),
            "T": np.array([0.0, 0.0, 0.5]),
            "W": np.array([0.25, 0.25, 0.25]),
            "X": np.array([-zeta, zeta, zeta]),
            "X_1": np.array([zeta, 1 - zeta, -zeta]),
            "Y": np.array([eta, -eta, eta]),
            "Y_1": np.array([1 - eta, eta, -eta]),
            "Z": np.array([0.5, 0.5, -0.5]),
        }
        path = [
            ["\\Gamma", "X", "L", "T", "W", "R", "X_1", "Z", "\\Gamma", "Y", "S", "W"],
            ["L_1", "Y"],
            ["Y_1", "Z"],
        ]
        return {"kpoints": kpoints, "path": path}

    def orcc(self, a, b, c):
        """
        ORCC Path
        """
        self.name = "ORCC"
        zeta = (1 + a ** 2 / b ** 2) / 4
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "A": np.array([zeta, zeta, 0.5]),
            "A_1": np.array([-zeta, 1 - zeta, 0.5]),
            "R": np.array([0.0, 0.5, 0.5]),
            "S": np.array([0.0, 0.5, 0.0]),
            "T": np.array([-0.5, 0.5, 0.5]),
            "X": np.array([zeta, zeta, 0.0]),
            "X_1": np.array([-zeta, 1 - zeta, 0.0]),
            "Y": np.array([-0.5, 0.5, 0]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            [
                "\\Gamma",
                "X",
                "S",
                "R",
                "A",
                "Z",
                "\\Gamma",
                "Y",
                "X_1",
                "A_1",
                "T",
                "Y",
            ],
            ["Z", "T"],
        ]
        return {"kpoints": kpoints, "path": path}

    def hex(self):
        """
        HEX Path
        """
        self.name = "HEX"
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "A": np.array([0.0, 0.0, 0.5]),
            "H": np.array([1.0 / 3.0, 1.0 / 3.0, 0.5]),
            "K": np.array([1.0 / 3.0, 1.0 / 3.0, 0.0]),
            "L": np.array([0.5, 0.0, 0.5]),
            "M": np.array([0.5, 0.0, 0.0]),
        }
        path = [
            ["\\Gamma", "M", "K", "\\Gamma", "A", "L", "H", "A"],
            ["L", "M"],
            ["K", "H"],
        ]
        return {"kpoints": kpoints, "path": path}

    def rhl1(self, alpha):
        """
        RHL1 Path
        """
        self.name = "RHL1"
        eta = (1 + 4 * cos(alpha)) / (2 + 4 * cos(alpha))
        nu = 3.0 / 4.0 - eta / 2.0
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "B": np.array([eta, 0.5, 1.0 - eta]),
            "B_1": np.array([1.0 / 2.0, 1.0 - eta, eta - 1.0]),
            "F": np.array([0.5, 0.5, 0.0]),
            "L": np.array([0.5, 0.0, 0.0]),
            "L_1": np.array([0.0, 0.0, -0.5]),
            "P": np.array([eta, nu, nu]),
            "P_1": np.array([1.0 - nu, 1.0 - nu, 1.0 - eta]),
            "P_2": np.array([nu, nu, eta - 1.0]),
            "Q": np.array([1.0 - nu, nu, 0.0]),
            "X": np.array([nu, 0.0, -nu]),
            "Z": np.array([0.5, 0.5, 0.5]),
        }
        path = [
            ["\\Gamma", "L", "B_1"],
            ["B", "Z", "\\Gamma", "X"],
            ["Q", "F", "P_1", "Z"],
            ["L", "P"],
        ]
        return {"kpoints": kpoints, "path": path}

    def rhl2(self, alpha):
        """
        RHL2 Path
        """
        self.name = "RHL2"
        eta = 1 / (2 * tan(alpha / 2.0) ** 2)
        nu = 3.0 / 4.0 - eta / 2.0
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "F": np.array([0.5, -0.5, 0.0]),
            "L": np.array([0.5, 0.0, 0.0]),
            "P": np.array([1 - nu, -nu, 1 - nu]),
            "P_1": np.array([nu, nu - 1.0, nu - 1.0]),
            "Q": np.array([eta, eta, eta]),
            "Q_1": np.array([1.0 - eta, -eta, -eta]),
            "Z": np.array([0.5, -0.5, 0.5]),
        }
        path = [["\\Gamma", "P", "Z", "Q", "\\Gamma", "F", "P_1", "Q_1", "L", "Z"]]
        return {"kpoints": kpoints, "path": path}

    def mcl(self, b, c, beta):
        """
        MCL Path
        """
        self.name = "MCL"
        eta = (1 - b * cos(beta) / c) / (2 * sin(beta) ** 2)
        nu = 0.5 - eta * c * cos(beta) / b
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "A": np.array([0.5, 0.5, 0.0]),
            "C": np.array([0.0, 0.5, 0.5]),
            "D": np.array([0.5, 0.0, 0.5]),
            "D_1": np.array([0.5, 0.5, -0.5]),
            "E": np.array([0.5, 0.5, 0.5]),
            "H": np.array([0.0, eta, 1.0 - nu]),
            "H_1": np.array([0.0, 1.0 - eta, nu]),
            "H_2": np.array([0.0, eta, -nu]),
            "M": np.array([0.5, eta, 1.0 - nu]),
            "M_1": np.array([0.5, 1 - eta, nu]),
            "M_2": np.array([0.5, 1 - eta, nu]),
            "X": np.array([0.0, 0.5, 0.0]),
            "Y": np.array([0.0, 0.0, 0.5]),
            "Y_1": np.array([0.0, 0.0, -0.5]),
            "Z": np.array([0.5, 0.0, 0.0]),
        }
        path = [
            ["\\Gamma", "Y", "H", "C", "E", "M_1", "A", "X", "H_1"],
            ["M", "D", "Z"],
            ["Y", "D"],
        ]
        return {"kpoints": kpoints, "path": path}

    def mclc1(self, a, b, c, alpha):
        """
        MCLC1 Path
        """
        self.name = "MCLC1"
        zeta = (2 - b * cos(alpha) / c) / (4 * sin(alpha) ** 2)
        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
        psi = 0.75 - a ** 2 / (4 * b ** 2 * sin(alpha) ** 2)
        phi = psi + (0.75 - psi) * b * cos(alpha) / c
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "N": np.array([0.5, 0.0, 0.0]),
            "N_1": np.array([0.0, -0.5, 0.0]),
            "F": np.array([1 - zeta, 1 - zeta, 1 - eta]),
            "F_1": np.array([zeta, zeta, eta]),
            "F_2": np.array([-zeta, -zeta, 1 - eta]),
            "I": np.array([phi, 1 - phi, 0.5]),
            "I_1": np.array([1 - phi, phi - 1, 0.5]),
            "L": np.array([0.5, 0.5, 0.5]),
            "M": np.array([0.5, 0.0, 0.5]),
            "X": np.array([1 - psi, psi - 1, 0.0]),
            "X_1": np.array([psi, 1 - psi, 0.0]),
            "X_2": np.array([psi - 1, -psi, 0.0]),
            "Y": np.array([0.5, 0.5, 0.0]),
            "Y_1": np.array([-0.5, -0.5, 0.0]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            ["\\Gamma", "Y", "F", "L", "I"],
            ["I_1", "Z", "F_1"],
            ["Y", "X_1"],
            ["X", "\\Gamma", "N"],
            ["M", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}

    def mclc2(self, a, b, c, alpha):
        """
        MCLC2 Path
        """
        self.name = "MCLC2"
        zeta = (2 - b * cos(alpha) / c) / (4 * sin(alpha) ** 2)
        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
        psi = 0.75 - a ** 2 / (4 * b ** 2 * sin(alpha) ** 2)
        phi = psi + (0.75 - psi) * b * cos(alpha) / c
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "N": np.array([0.5, 0.0, 0.0]),
            "N_1": np.array([0.0, -0.5, 0.0]),
            "F": np.array([1 - zeta, 1 - zeta, 1 - eta]),
            "F_1": np.array([zeta, zeta, eta]),
            "F_2": np.array([-zeta, -zeta, 1 - eta]),
            "F_3": np.array([1 - zeta, -zeta, 1 - eta]),
            "I": np.array([phi, 1 - phi, 0.5]),
            "I_1": np.array([1 - phi, phi - 1, 0.5]),
            "L": np.array([0.5, 0.5, 0.5]),
            "M": np.array([0.5, 0.0, 0.5]),
            "X": np.array([1 - psi, psi - 1, 0.0]),
            "X_1": np.array([psi, 1 - psi, 0.0]),
            "X_2": np.array([psi - 1, -psi, 0.0]),
            "Y": np.array([0.5, 0.5, 0.0]),
            "Y_1": np.array([-0.5, -0.5, 0.0]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            ["\\Gamma", "Y", "F", "L", "I"],
            ["I_1", "Z", "F_1"],
            ["N", "\\Gamma", "M"],
        ]
        return {"kpoints": kpoints, "path": path}

    def mclc3(self, a, b, c, alpha):
        """
        MCLC3 Path
        """
        self.name = "MCLC3"
        mu = (1 + b ** 2 / a ** 2) / 4.0
        delta = b * c * cos(alpha) / (2 * a ** 2)
        zeta = mu - 0.25 + (1 - b * cos(alpha) / c) / (4 * sin(alpha) ** 2)
        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
        phi = 1 + zeta - 2 * mu
        psi = eta - 2 * delta
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "F": np.array([1 - phi, 1 - phi, 1 - psi]),
            "F_1": np.array([phi, phi - 1, psi]),
            "F_2": np.array([1 - phi, -phi, 1 - psi]),
            "H": np.array([zeta, zeta, eta]),
            "H_1": np.array([1 - zeta, -zeta, 1 - eta]),
            "H_2": np.array([-zeta, -zeta, 1 - eta]),
            "I": np.array([0.5, -0.5, 0.5]),
            "M": np.array([0.5, 0.0, 0.5]),
            "N": np.array([0.5, 0.0, 0.0]),
            "N_1": np.array([0.0, -0.5, 0.0]),
            "X": np.array([0.5, -0.5, 0.0]),
            "Y": np.array([mu, mu, delta]),
            "Y_1": np.array([1 - mu, -mu, -delta]),
            "Y_2": np.array([-mu, -mu, -delta]),
            "Y_3": np.array([mu, mu - 1, delta]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            ["\\Gamma", "Y", "F", "H", "Z", "I", "F_1"],
            ["H_1", "Y_1", "X", "\\Gamma", "N"],
            ["M", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}

    def mclc4(self, a, b, c, alpha):
        """
        MCLC4 Path
        """
        self.name = "MCLC4"
        mu = (1 + b ** 2 / a ** 2) / 4.0
        delta = b * c * cos(alpha) / (2 * a ** 2)
        zeta = mu - 0.25 + (1 - b * cos(alpha) / c) / (4 * sin(alpha) ** 2)
        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
        phi = 1 + zeta - 2 * mu
        psi = eta - 2 * delta
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "F": np.array([1 - phi, 1 - phi, 1 - psi]),
            "F_1": np.array([phi, phi - 1, psi]),
            "F_2": np.array([1 - phi, -phi, 1 - psi]),
            "H": np.array([zeta, zeta, eta]),
            "H_1": np.array([1 - zeta, -zeta, 1 - eta]),
            "H_2": np.array([-zeta, -zeta, 1 - eta]),
            "I": np.array([0.5, -0.5, 0.5]),
            "M": np.array([0.5, 0.0, 0.5]),
            "N": np.array([0.5, 0.0, 0.0]),
            "N_1": np.array([0.0, -0.5, 0.0]),
            "X": np.array([0.5, -0.5, 0.0]),
            "Y": np.array([mu, mu, delta]),
            "Y_1": np.array([1 - mu, -mu, -delta]),
            "Y_2": np.array([-mu, -mu, -delta]),
            "Y_3": np.array([mu, mu - 1, delta]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            ["\\Gamma", "Y", "F", "H", "Z", "I"],
            ["H_1", "Y_1", "X", "\\Gamma", "N"],
            ["M", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}

    def mclc5(self, a, b, c, alpha):
        """
        MCLC5 Path
        """
        self.name = "MCLC5"
        zeta = (b ** 2 / a ** 2 + (1 - b * cos(alpha) / c) / sin(alpha) ** 2) / 4
        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
        mu = eta / 2 + b ** 2 / (4 * a ** 2) - b * c * cos(alpha) / (2 * a ** 2)
        nu = 2 * mu - zeta
        rho = 1 - zeta * a ** 2 / b ** 2
        omega = (4 * nu - 1 - b ** 2 * sin(alpha) ** 2 / a ** 2) * c / (2 * b * cos(alpha))
        delta = zeta * c * cos(alpha) / b + omega / 2 - 0.25
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "F": np.array([nu, nu, omega]),
            "F_1": np.array([1 - nu, 1 - nu, 1 - omega]),
            "F_2": np.array([nu, nu - 1, omega]),
            "H": np.array([zeta, zeta, eta]),
            "H_1": np.array([1 - zeta, -zeta, 1 - eta]),
            "H_2": np.array([-zeta, -zeta, 1 - eta]),
            "I": np.array([rho, 1 - rho, 0.5]),
            "I_1": np.array([1 - rho, rho - 1, 0.5]),
            "L": np.array([0.5, 0.5, 0.5]),
            "M": np.array([0.5, 0.0, 0.5]),
            "N": np.array([0.5, 0.0, 0.0]),
            "N_1": np.array([0.0, -0.5, 0.0]),
            "X": np.array([0.5, -0.5, 0.0]),
            "Y": np.array([mu, mu, delta]),
            "Y_1": np.array([1 - mu, -mu, -delta]),
            "Y_2": np.array([-mu, -mu, -delta]),
            "Y_3": np.array([mu, mu - 1, delta]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            ["\\Gamma", "Y", "F", "L", "I"],
            ["I_1", "Z", "H", "F_1"],
            ["H_1", "Y_1", "X", "\\Gamma", "N"],
            ["M", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}

    def tria(self):
        """
        TRI1a Path
        """
        self.name = "TRI1a"
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "L": np.array([0.5, 0.5, 0.0]),
            "M": np.array([0.0, 0.5, 0.5]),
            "N": np.array([0.5, 0.0, 0.5]),
            "R": np.array([0.5, 0.5, 0.5]),
            "X": np.array([0.5, 0.0, 0.0]),
            "Y": np.array([0.0, 0.5, 0.0]),
            "Z": np.array([0.0, 0.0, 0.5]),
        }
        path = [
            ["X", "\\Gamma", "Y"],
            ["L", "\\Gamma", "Z"],
            ["N", "\\Gamma", "M"],
            ["R", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}

    def trib(self):
        """
        TRI1b Path
        """
        self.name = "TRI1b"
        kpoints = {
            "\\Gamma": np.array([0.0, 0.0, 0.0]),
            "L": np.array([0.5, -0.5, 0.0]),
            "M": np.array([0.0, 0.0, 0.5]),
            "N": np.array([-0.5, -0.5, 0.5]),
            "R": np.array([0.0, -0.5, 0.5]),
            "X": np.array([0.0, -0.5, 0.0]),
            "Y": np.array([0.5, 0.0, 0.0]),
            "Z": np.array([-0.5, 0.0, 0.5]),
        }
        path = [
            ["X", "\\Gamma", "Y"],
            ["L", "\\Gamma", "Z"],
            ["N", "\\Gamma", "M"],
            ["R", "\\Gamma"],
        ]
        return {"kpoints": kpoints, "path": path}


class KPathSeek(KPathBase):
    """
    This class looks for path along high symmetry lines in
    the Brillouin Zone.
    It is based on Hinuma, Y., Pizzi, G., Kumagai, Y., Oba, F.,
    & Tanaka, I. (2017). Band structure diagram paths based on
    crystallography. Computational Materials Science, 128, 140–184.
    https://doi.org/10.1016/j.commatsci.2016.10.015
    It should be used with primitive structures that
    comply with the definition from the paper.
    The symmetry is determined by spglib through the
    SpacegroupAnalyzer class. KPoints from get_kpoints() method
    are returned in the reciprocal cell basis defined in the paper.
    """

    @requires(
        get_path is not None,
        "SeeK-path is required to use the convention by Hinuma et al.",
    )
    def __init__(self, structure, symprec=0.01, angle_tolerance=5, atol=1e-5, system_is_tri=True):
        """
        Args:
            structure (Structure): Structure object
            symprec (float): Tolerance for symmetry finding
            angle_tolerance (float): Angle tolerance for symmetry finding.
            atol (float): Absolute tolerance used to determine edge cases
                for settings of structures.
            system_is_tri (boolean): Indicates if the system is time-reversal
                invariant.
        """

        super().__init__(structure, symprec=symprec, angle_tolerance=angle_tolerance, atol=atol)

        positions = structure.frac_coords

        sp = structure.site_properties
        species = [site.species for site in structure]
        site_data = species

        if not system_is_tri:
            warn("Non-zero 'magmom' data will be used to define unique atoms in the cell.")
            site_data = zip(species, [tuple(vec) for vec in sp["magmom"]])

        unique_species = []
        numbers = []

        for species, g in itertools.groupby(site_data):
            if species in unique_species:
                ind = unique_species.index(species)
                numbers.extend([ind + 1] * len(tuple(g)))
            else:
                unique_species.append(species)
                numbers.extend([len(unique_species)] * len(tuple(g)))

        cell = (self._latt.matrix, positions, numbers)

        lattice, scale_pos, atom_num = spglib.standardize_cell(
            cell, to_primitive=False, no_idealize=True, symprec=symprec
        )

        spg_struct = (lattice, scale_pos, atom_num)
        spath_dat = get_path(spg_struct, system_is_tri, "hpkot", atol, symprec, angle_tolerance)

        self._tmat = self._trans_sc_to_Hin(spath_dat["bravais_lattice_extended"])
        self._rec_lattice = Lattice(spath_dat["reciprocal_primitive_lattice"])

        spath_data_formatted = [[spath_dat["path"][0][0]]]
        count = 0
        for pnum in range(len(spath_dat["path"]) - 1):
            if spath_dat["path"][pnum][1] == spath_dat["path"][pnum + 1][0]:
                spath_data_formatted[count].append(spath_dat["path"][pnum][1])
            else:
                spath_data_formatted[count].append(spath_dat["path"][pnum][1])
                spath_data_formatted.append([])
                count += 1
                spath_data_formatted[count].append(spath_dat["path"][pnum + 1][0])

        spath_data_formatted[-1].append(spath_dat["path"][-1][1])

        self._kpath = {
            "kpoints": spath_dat["point_coords"],
            "path": spath_data_formatted,
        }

    @staticmethod
    def _trans_sc_to_Hin(sub_class):

        if sub_class in [
            "cP1",
            "cP2",
            "cF1",
            "cF2",
            "cI1",
            "tP1",
            "oP1",
            "hP1",
            "hP2",
            "tI1",
            "tI2",
            "oF1",
            "oF3",
            "oI1",
            "oI3",
            "oC1",
            "hR1",
            "hR2",
            "aP1",
            "aP2",
            "aP3",
            "oA1",
        ]:
            return np.eye(3)
        if sub_class == "oF2":
            return np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]])
        if sub_class == "oI2":
            return np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]])
        if sub_class == "oI3":
            return np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]])
        if sub_class == "oA2":
            return np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]])
        if sub_class == "oC2":
            return np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]])
        if sub_class in ["mP1", "mC1", "mC2", "mC3"]:
            return np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]])
        raise RuntimeError("Sub-classification of crystal not found!")


class KPathLatimerMunro(KPathBase):
    """
    This class looks for a path along high symmetry lines in the
    Brillouin zone. It is based on the method outlined in:

    npj Comput Mater 6, 112 (2020). 10.1038/s41524-020-00383-7

    The user should ensure that the lattice of the input structure
    is as reduced as possible, i.e. that there is no linear
    combination of lattice vectors which can produce a vector of
    lesser magnitude than the given set (this is required to
    obtain the correct Brillouin zone within the current
    implementaiton). This is checked during initialization and a
    warning is issued if the condition is not fulfilled.
    In the case of magnetic structures, care must also be taken to
    provide the magnetic primitive cell (i.e. that which reproduces
    the entire crystal, including the correct magnetic ordering,
    upon application of lattice translations). There is no way to
    programatically check for this, so if the input structure is
    incorrect, the class will output the incorrect kpath without
    any warning being issued.
    """

    def __init__(
        self,
        structure,
        has_magmoms=False,
        magmom_axis=None,
        symprec=0.01,
        angle_tolerance=5,
        atol=1e-5,
    ):
        """
        Args:
            structure (Structure): Structure object
            has_magmoms (boolean): Whether the input structure contains
                magnetic moments as site properties with the key 'magmom.'
                Values may be in the form of 3-component vectors given in
                the basis of the input lattice vectors, or as scalars, in
                which case the spin axis will default to a_3, the third
                real-space lattice vector (this triggers a warning).
            magmom_axis (list or numpy array): 3-component vector specifying
                direction along which magnetic moments given as scalars
                should point. If all magnetic moments are provided as
                vectors then this argument is not used.
            symprec (float): Tolerance for symmetry finding
            angle_tolerance (float): Angle tolerance for symmetry finding.
            atol (float): Absolute tolerance used to determine symmetric
                equivalence of points and lines on the BZ.
        """
        super().__init__(structure, symprec=symprec, angle_tolerance=angle_tolerance, atol=atol)

        # Check to see if input lattice is reducible. Ref: B Gruber in Acta. Cryst. Vol. A29,
        # pp. 433-440 ('The Relationship between Reduced Cells in a General Bravais lattice').
        # The correct BZ will still be obtained if the lattice vectors are reducible by any
        # linear combination of themselves with coefficients of absolute value less than 2,
        # hence a missing factor of 2 as compared to the reference.
        reducible = []
        for i in range(3):
            for j in range(3):
                if i != j:
                    if (
                        np.absolute(np.dot(self._latt.matrix[i], self._latt.matrix[j]))
                        > np.dot(self._latt.matrix[i], self._latt.matrix[i])
                        and np.absolute(
                            np.dot(self._latt.matrix[i], self._latt.matrix[j])
                            - np.dot(self._latt.matrix[i], self._latt.matrix[i])
                        )
                        > atol
                    ):
                        reducible.append(True)
                    else:
                        reducible.append(False)
        if np.any(reducible):
            print("reducible")
            warn(
                "The lattice of the input structure is not fully reduced!"
                "The path can be incorrect. Use at your own risk."
            )

        if magmom_axis is None:
            magmom_axis = np.array([0, 0, 1])
            axis_specified = False
        else:
            axis_specified = True

        self._kpath = self._get_ksymm_kpath(has_magmoms, magmom_axis, axis_specified, symprec, angle_tolerance, atol)

    @property
    def mag_type(self):
        """
        Returns:
        The type of magnetic space group as a string.
        Current implementation does not distinguish
        between types 3 and 4, so return value is '3/4'.
        If has_magmoms is False, returns '0'.
        """
        return self._mag_type

    def _get_ksymm_kpath(self, has_magmoms, magmom_axis, axis_specified, symprec, angle_tolerance, atol):
        ID = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
        # parity, aka the inversion operation (not calling it
        PAR = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, -1]])
        # INV to avoid confusion with np.linalg.inv() function)

        # 1: Get lattices of real and reciprocal structures, and reciprocal
        # point group, and Brillouin zone (BZ)

        V = self._latt.matrix.T  # fractional real space to cartesian real space
        # fractional reciprocal space to cartesian reciprocal space
        W = self._rec_lattice.matrix.T
        # fractional real space to fractional reciprocal space
        A = np.dot(np.linalg.inv(W), V)

        if has_magmoms:
            grey_struct = self._structure.copy()
            grey_struct.remove_site_property("magmom")
            sga = SpacegroupAnalyzer(grey_struct, symprec=symprec, angle_tolerance=angle_tolerance)
            grey_ops = sga.get_symmetry_operations()
            self._structure = self._convert_all_magmoms_to_vectors(magmom_axis, axis_specified)
            mag_ops = self._get_magnetic_symmetry_operations(self._structure, grey_ops, atol)

            D = [
                SymmOp.from_rotation_and_translation(
                    rotation_matrix=op.rotation_matrix,
                    translation_vec=op.translation_vector,
                )
                for op in mag_ops
                if op.time_reversal == 1
            ]

            fD = [
                SymmOp.from_rotation_and_translation(
                    rotation_matrix=op.rotation_matrix,
                    translation_vec=op.translation_vector,
                )
                for op in mag_ops
                if op.time_reversal == -1
            ]

            if np.array([m == np.array([0, 0, 0]) for m in self._structure.site_properties["magmom"]]).all():
                fD = D
                D = []

            if len(fD) == 0:  # no operations contain time reversal; type 1
                self._mag_type = "1"
                isomorphic_point_group = [d.rotation_matrix for d in D]
                recip_point_group = self._get_reciprocal_point_group(isomorphic_point_group, ID, A)
            elif len(D) == 0:  # all operations contain time reversal / all magmoms zero; type 2
                self._mag_type = "2"
                isomorphic_point_group = [d.rotation_matrix for d in fD]
                recip_point_group = self._get_reciprocal_point_group(isomorphic_point_group, PAR, A)
            else:  # half and half; type 3 or 4
                self._mag_type = "3/4"
                f = self._get_coset_factor(D + fD, D)
                isomorphic_point_group = [d.rotation_matrix for d in D]
                recip_point_group = self._get_reciprocal_point_group(
                    isomorphic_point_group, np.dot(PAR, f.rotation_matrix), A
                )
        else:
            self._mag_type = "0"
            if "magmom" in self._structure.site_properties:
                warn(
                    "The parameter has_magmoms is False, but site_properties contains the key magmom."
                    "This property will be removed and could result in different symmetry operations."
                )
                self._structure.remove_site_property("magmom")
            sga = SpacegroupAnalyzer(self._structure)
            ops = sga.get_symmetry_operations()
            isomorphic_point_group = [op.rotation_matrix for op in ops]

            recip_point_group = self._get_reciprocal_point_group(isomorphic_point_group, PAR, A)

        self._rpg = recip_point_group

        # 2: Get all vertices, edge- and face- center points of BZ ("key points")

        key_points, bz_as_key_point_inds, face_center_inds = self._get_key_points()

        # 3: Find symmetry-equivalent points, which can be mapped to each other by a combination of point group
        # operations and integer translations by lattice vectors. The integers will only be -1, 0, or 1, since
        # we are restricting to the BZ.

        key_points_inds_orbits = self._get_key_point_orbits(key_points=key_points)

        # 4: Get all lines on BZ between adjacent key points and between gamma
        # and key points ("key lines")

        key_lines = self._get_key_lines(key_points=key_points, bz_as_key_point_inds=bz_as_key_point_inds)

        # 5: Find symmetry-equivalent key lines, defined as endpoints of first line being equivalent
        # to end points of second line, and a random point in between being equivalent to the mapped
        # random point.

        key_lines_inds_orbits = self._get_key_line_orbits(
            key_points=key_points,
            key_lines=key_lines,
            key_points_inds_orbits=key_points_inds_orbits,
        )

        # 6 & 7: Get little groups for key points (group of symmetry elements present at that point).
        # Get little groups for key lines (group of symmetry elements present at every point
        # along the line). This is implemented by testing the symmetry at a point e/pi of the
        # way between the two endpoints.

        little_groups_points, little_groups_lines = self._get_little_groups(
            key_points=key_points,
            key_points_inds_orbits=key_points_inds_orbits,
            key_lines_inds_orbits=key_lines_inds_orbits,
        )

        # 8: Choose key lines for k-path. Loose criteria set: choose any points / segments
        # with spatial symmetry greater than the general position (AKA more symmetry operations
        # than just the identity or identity * TR in the little group).
        # This function can be edited to alter high-symmetry criteria for choosing points and lines

        point_orbits_in_path, line_orbits_in_path = self._choose_path(
            key_points=key_points,
            key_points_inds_orbits=key_points_inds_orbits,
            key_lines_inds_orbits=key_lines_inds_orbits,
            little_groups_points=little_groups_points,
            little_groups_lines=little_groups_lines,
        )

        # 10: Consolidate selected segments into a single irreducible section of the Brilouin zone (as determined
        # by the reciprocal point and lattice symmetries). This is accomplished by identifying the boundary
        # planes of the IRBZ. Also, get labels for points according to distance away from axes.

        IRBZ_points_inds = self._get_IRBZ(recip_point_group, W, key_points, face_center_inds, atol)
        lines_in_path_inds = []
        for ind in line_orbits_in_path:
            for tup in key_lines_inds_orbits[ind]:
                if tup[0] in IRBZ_points_inds and tup[1] in IRBZ_points_inds:
                    lines_in_path_inds.append(tup)
                    break
        G = nx.Graph(lines_in_path_inds)
        lines_in_path_inds = list(nx.edge_dfs(G))
        points_in_path_inds = [ind for tup in lines_in_path_inds for ind in tup]
        points_in_path_inds_unique = list(set(points_in_path_inds))

        orbit_cosines = []
        for i, orbit in enumerate(key_points_inds_orbits[:-1]):

            orbit_cosines.append(
                sorted(
                    sorted(
                        [
                            (
                                j,
                                np.round(
                                    np.dot(key_points[k], self.LabelPoints(j))
                                    / (np.linalg.norm(key_points[k]) * np.linalg.norm(self.LabelPoints(j))),
                                    decimals=3,
                                ),
                            )
                            for k in orbit
                            for j in range(26)
                        ],
                        key=operator.itemgetter(0),
                    ),
                    key=operator.itemgetter(1),
                    reverse=True,
                )
            )

        orbit_labels = self._get_orbit_labels(orbit_cosines, key_points_inds_orbits, atol)
        key_points_labels = ["" for i in range(len(key_points))]
        for i, orbit in enumerate(key_points_inds_orbits):
            for point_ind in orbit:
                key_points_labels[point_ind] = self.LabelSymbol(int(orbit_labels[i]))

        kpoints = {}
        reverse_kpoints = {}
        for point_ind in points_in_path_inds_unique:
            point_label = key_points_labels[point_ind]
            if point_label not in kpoints.keys():
                kpoints[point_label] = key_points[point_ind]
                reverse_kpoints[point_ind] = point_label
            else:
                existing_labels = [key for key in kpoints.keys() if point_label in key]
                if "'" not in point_label:
                    existing_labels[:] = [label for label in existing_labels if "'" not in label]
                if len(existing_labels) == 1:
                    max_occurence = 0
                else:
                    if "'" not in point_label:
                        max_occurence = max([int(label[3:-1]) for label in existing_labels[1:]])
                    else:
                        max_occurence = max([int(label[4:-1]) for label in existing_labels[1:]])
                kpoints[point_label + "_{" + str(max_occurence + 1) + "}"] = key_points[point_ind]
                reverse_kpoints[point_ind] = point_label + "_{" + str(max_occurence + 1) + "}"

        path = []
        i = 0
        start_of_subpath = True
        while i < len(points_in_path_inds):
            if start_of_subpath:
                path.append([reverse_kpoints[points_in_path_inds[i]]])
                i += 1
                start_of_subpath = False
            elif points_in_path_inds[i] == points_in_path_inds[i + 1]:
                path[-1].append(reverse_kpoints[points_in_path_inds[i]])
                i += 2
            else:
                path[-1].append(reverse_kpoints[points_in_path_inds[i]])
                i += 1
                start_of_subpath = True
            if i == len(points_in_path_inds) - 1:
                path[-1].append(reverse_kpoints[points_in_path_inds[i]])
                i += 1

        return {"kpoints": kpoints, "path": path}

    def _choose_path(
        self,
        key_points,
        key_points_inds_orbits,
        key_lines_inds_orbits,
        little_groups_points,
        little_groups_lines,
    ):
        #
        # This function can be edited to alter high-symmetry criteria for choosing points and lines
        #

        ID = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
        PAR = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, -1]])

        gamma_ind = len(key_points) - 1

        line_orbits_in_path = []
        point_orbits_in_path = []
        for (i, little_group) in enumerate(little_groups_lines):
            add_rep = False
            nC2 = 0
            nC3 = 0
            nsig = 0
            for j, opind in enumerate(little_group):
                op = self._rpg[opind]
                if not (op == ID).all():
                    if (np.dot(op, op) == ID).all():
                        if np.linalg.det(op) == 1:
                            nC2 += 1
                            break
                        if not (op == PAR).all():
                            nsig += 1
                            break
                    elif (np.dot(op, np.dot(op, op)) == ID).all():
                        nC3 += 1
                        break
            if nC2 > 0 or nC3 > 0 or nsig > 0:
                add_rep = True

            if add_rep:
                line_orbits_in_path.append(i)
                l = key_lines_inds_orbits[i][0]
                ind0 = l[0]
                ind1 = l[1]
                found0 = False
                found1 = False
                for (j, orbit) in enumerate(key_points_inds_orbits):
                    if ind0 in orbit:
                        point_orbits_in_path.append(j)
                        found0 = True
                    if ind1 in orbit:
                        point_orbits_in_path.append(j)
                        found1 = True
                    if found0 and found1:
                        break
        point_orbits_in_path = list(set(point_orbits_in_path))

        # Choose remaining unconnected key points for k-path. The ones that remain are
        # those with inversion symmetry. Connect them to gamma.

        unconnected = []

        for i in range(len(key_points_inds_orbits)):
            if i not in point_orbits_in_path:
                unconnected.append(i)

        for ind in unconnected:
            connect = False
            for op_ind in little_groups_points[ind]:
                op = self._rpg[op_ind]
                if (op == ID).all():
                    pass
                elif (op == PAR).all():
                    connect = True
                    break
                elif np.linalg.det(op) == 1:
                    if (np.dot(op, np.dot(op, op)) == ID).all():
                        pass
                    else:
                        connect = True
                        break
                else:
                    pass
            if connect:
                l = (key_points_inds_orbits[ind][0], gamma_ind)
                for (j, orbit) in enumerate(key_lines_inds_orbits):
                    if l in orbit:
                        line_orbits_in_path.append(j)
                        break
                if gamma_ind not in point_orbits_in_path:
                    point_orbits_in_path.append(gamma_ind)
                point_orbits_in_path.append(ind)

        return point_orbits_in_path, line_orbits_in_path

    def _get_key_points(self):
        decimals = ceil(-1 * np.log10(self._atol)) - 1
        bz = self._rec_lattice.get_wigner_seitz_cell()

        key_points = []
        face_center_inds = []
        bz_as_key_point_inds = []

        # pymatgen gives BZ in cartesian coordinates; convert to fractional in
        # the primitive basis for reciprocal space
        for (i, facet) in enumerate(bz):
            for (j, vert) in enumerate(facet):
                vert = self._rec_lattice.get_fractional_coords(vert)
                bz[i][j] = vert
        pop = []
        for i, facet in enumerate(bz):
            rounded_facet = np.around(facet, decimals=decimals)
            u, indices = np.unique(rounded_facet, axis=0, return_index=True)
            if len(u) in [1, 2]:
                pop.append(i)
            else:
                bz[i] = [facet[j] for j in np.sort(indices)]
        bz = [bz[i] for i in range(len(bz)) if i not in pop]

        # use vertex points to calculate edge- and face- centers
        for (i, facet) in enumerate(bz):
            bz_as_key_point_inds.append([])
            for (j, vert) in enumerate(facet):
                edge_center = (vert + facet[j + 1]) / 2.0 if j != len(facet) - 1 else (vert + facet[0]) / 2.0
                duplicatevert = False
                duplicateedge = False
                for (k, point) in enumerate(key_points):
                    if np.allclose(vert, point, atol=self._atol):
                        bz_as_key_point_inds[i].append(k)
                        duplicatevert = True
                        break
                for (k, point) in enumerate(key_points):
                    if np.allclose(edge_center, point, atol=self._atol):
                        bz_as_key_point_inds[i].append(k)
                        duplicateedge = True
                        break
                if not duplicatevert:
                    key_points.append(vert)
                    bz_as_key_point_inds[i].append(len(key_points) - 1)
                if not duplicateedge:
                    key_points.append(edge_center)
                    bz_as_key_point_inds[i].append(len(key_points) - 1)
            if len(facet) == 4:  # parallelogram facet
                face_center = (facet[0] + facet[1] + facet[2] + facet[3]) / 4.0
                key_points.append(face_center)
                face_center_inds.append(len(key_points) - 1)
                bz_as_key_point_inds[i].append(len(key_points) - 1)
            else:  # hexagonal facet
                face_center = (facet[0] + facet[1] + facet[2] + facet[3] + facet[4] + facet[5]) / 6.0
                key_points.append(face_center)
                face_center_inds.append(len(key_points) - 1)
                bz_as_key_point_inds[i].append(len(key_points) - 1)

        # add gamma point
        key_points.append(np.array([0, 0, 0]))
        return key_points, bz_as_key_point_inds, face_center_inds

    def _get_key_point_orbits(self, key_points):
        key_points_copy = dict(zip(range(len(key_points) - 1), key_points[0 : len(key_points) - 1]))
        # gamma not equivalent to any on BZ and is last point added to
        # key_points
        key_points_inds_orbits = []

        i = 0
        while len(key_points_copy) > 0:
            key_points_inds_orbits.append([])
            k0ind = list(key_points_copy.keys())[0]
            k0 = key_points_copy[k0ind]
            key_points_inds_orbits[i].append(k0ind)
            key_points_copy.pop(k0ind)

            for op in self._rpg:
                to_pop = []
                k1 = np.dot(op, k0)
                for ind_key in key_points_copy:
                    diff = k1 - key_points_copy[ind_key]
                    if self._all_ints(diff, atol=self._atol):
                        key_points_inds_orbits[i].append(ind_key)
                        to_pop.append(ind_key)

                for key in to_pop:
                    key_points_copy.pop(key)
            i += 1

        key_points_inds_orbits.append([len(key_points) - 1])

        return key_points_inds_orbits

    @staticmethod
    def _get_key_lines(key_points, bz_as_key_point_inds):
        key_lines = []
        gamma_ind = len(key_points) - 1

        for (i, facet_as_key_point_inds) in enumerate(bz_as_key_point_inds):
            facet_as_key_point_inds_bndy = facet_as_key_point_inds[: len(facet_as_key_point_inds) - 1]
            # not the face center point (don't need to check it since it's not
            # shared with other facets)
            face_center_ind = facet_as_key_point_inds[-1]
            for (j, ind) in enumerate(facet_as_key_point_inds_bndy):
                if (
                    min(ind, facet_as_key_point_inds_bndy[j - 1]),
                    max(ind, facet_as_key_point_inds_bndy[j - 1]),
                ) not in key_lines:
                    key_lines.append(
                        (
                            min(ind, facet_as_key_point_inds_bndy[j - 1]),
                            max(ind, facet_as_key_point_inds_bndy[j - 1]),
                        )
                    )
                k = j + 1 if j != len(facet_as_key_point_inds_bndy) - 1 else 0
                if (
                    min(ind, facet_as_key_point_inds_bndy[k]),
                    max(ind, facet_as_key_point_inds_bndy[k]),
                ) not in key_lines:
                    key_lines.append(
                        (
                            min(ind, facet_as_key_point_inds_bndy[k]),
                            max(ind, facet_as_key_point_inds_bndy[k]),
                        )
                    )
                if (ind, gamma_ind) not in key_lines:
                    key_lines.append((ind, gamma_ind))
                key_lines.append((min(ind, face_center_ind), max(ind, face_center_ind)))
            key_lines.append((face_center_ind, gamma_ind))

        return key_lines

    def _get_key_line_orbits(self, key_points, key_lines, key_points_inds_orbits):
        key_lines_copy = dict(zip(range(len(key_lines)), key_lines))
        key_lines_inds_orbits = []

        i = 0
        while len(key_lines_copy) > 0:
            key_lines_inds_orbits.append([])
            l0ind = list(key_lines_copy.keys())[0]
            l0 = key_lines_copy[l0ind]
            key_lines_inds_orbits[i].append(l0)
            key_lines_copy.pop(l0ind)
            to_pop = []
            p00 = key_points[l0[0]]
            p01 = key_points[l0[1]]
            pmid0 = p00 + e / pi * (p01 - p00)
            for ind_key in key_lines_copy:

                l1 = key_lines_copy[ind_key]
                p10 = key_points[l1[0]]
                p11 = key_points[l1[1]]
                equivptspar = False
                equivptsperp = False
                equivline = False

                if (
                    np.array([l0[0] in orbit and l1[0] in orbit for orbit in key_points_inds_orbits]).any()
                    and np.array([l0[1] in orbit and l1[1] in orbit for orbit in key_points_inds_orbits]).any()
                ):
                    equivptspar = True
                elif (
                    np.array([l0[1] in orbit and l1[0] in orbit for orbit in key_points_inds_orbits]).any()
                    and np.array([l0[0] in orbit and l1[1] in orbit for orbit in key_points_inds_orbits]).any()
                ):
                    equivptsperp = True

                if equivptspar:
                    pmid1 = p10 + e / pi * (p11 - p10)
                    for op in self._rpg:
                        if not equivline:
                            p00pr = np.dot(op, p00)
                            diff0 = p10 - p00pr
                            if self._all_ints(diff0, atol=self._atol):
                                pmid0pr = np.dot(op, pmid0) + diff0
                                p01pr = np.dot(op, p01) + diff0
                                if np.allclose(p11, p01pr, atol=self._atol) and np.allclose(
                                    pmid1, pmid0pr, atol=self._atol
                                ):
                                    equivline = True

                elif equivptsperp:
                    pmid1 = p11 + e / pi * (p10 - p11)
                    for op in self._rpg:
                        if not equivline:
                            p00pr = np.dot(op, p00)
                            diff0 = p11 - p00pr
                            if self._all_ints(diff0, atol=self._atol):
                                pmid0pr = np.dot(op, pmid0) + diff0
                                p01pr = np.dot(op, p01) + diff0
                                if np.allclose(p10, p01pr, atol=self._atol) and np.allclose(
                                    pmid1, pmid0pr, atol=self._atol
                                ):
                                    equivline = True

                if equivline:
                    key_lines_inds_orbits[i].append(l1)
                    to_pop.append(ind_key)

            for key in to_pop:
                key_lines_copy.pop(key)
            i += 1

        return key_lines_inds_orbits

    def _get_little_groups(self, key_points, key_points_inds_orbits, key_lines_inds_orbits):

        little_groups_points = []  # elements are lists of indicies of recip_point_group. the
        # list little_groups_points[i] is the little group for the
        # orbit key_points_inds_orbits[i]
        for (i, orbit) in enumerate(key_points_inds_orbits):
            k0 = key_points[orbit[0]]
            little_groups_points.append([])
            for (j, op) in enumerate(self._rpg):
                gamma_to = np.dot(op, -1 * k0) + k0
                check_gamma = True
                if not self._all_ints(gamma_to, atol=self._atol):
                    check_gamma = False
                if check_gamma:
                    little_groups_points[i].append(j)

        # elements are lists of indicies of recip_point_group. the list
        # little_groups_lines[i] is
        little_groups_lines = []
        # the little group for the orbit key_points_inds_lines[i]

        for (i, orbit) in enumerate(key_lines_inds_orbits):
            l0 = orbit[0]
            v = key_points[l0[1]] - key_points[l0[0]]
            k0 = key_points[l0[0]] + np.e / pi * v
            little_groups_lines.append([])
            for (j, op) in enumerate(self._rpg):
                gamma_to = np.dot(op, -1 * k0) + k0
                check_gamma = True
                if not self._all_ints(gamma_to, atol=self._atol):
                    check_gamma = False
                if check_gamma:
                    little_groups_lines[i].append(j)

        return little_groups_points, little_groups_lines

    def _convert_all_magmoms_to_vectors(self, magmom_axis, axis_specified):
        struct = self._structure.copy()
        magmom_axis = np.array(magmom_axis)
        if "magmom" not in struct.site_properties:
            warn(
                "The 'magmom' property is not set in the structure's site properties."
                "All magnetic moments are being set to zero."
            )
            struct.add_site_property("magmom", [np.array([0, 0, 0]) for i in range(len(struct.sites))])

            return struct

        old_magmoms = struct.site_properties["magmom"]
        new_magmoms = []
        found_scalar = False

        for magmom in old_magmoms:
            if isinstance(magmom, np.ndarray):
                new_magmoms.append(magmom)
            elif isinstance(magmom, list):
                new_magmoms.append(np.array(magmom))
            else:
                found_scalar = True
                new_magmoms.append(magmom * magmom_axis)

        if found_scalar and not axis_specified:
            warn("At least one magmom had a scalar value and magmom_axis was not specified." "Defaulted to z+ spinor.")

        struct.remove_site_property("magmom")
        struct.add_site_property("magmom", new_magmoms)
        return struct

    def _get_magnetic_symmetry_operations(self, struct, grey_ops, atol):
        mag_ops = []
        magmoms = struct.site_properties["magmom"]
        nonzero_magmom_inds = [i for i in range(len(struct.sites)) if not (magmoms[i] == np.array([0, 0, 0])).all()]
        init_magmoms = [site.properties["magmom"] for (i, site) in enumerate(struct.sites) if i in nonzero_magmom_inds]
        sites = [site for (i, site) in enumerate(struct.sites) if i in nonzero_magmom_inds]
        init_site_coords = [site.frac_coords for site in sites]
        for op in grey_ops:
            r = op.rotation_matrix
            t = op.translation_vector
            xformed_magmoms = [self._apply_op_to_magmom(r, magmom) for magmom in init_magmoms]
            xformed_site_coords = [np.dot(r, site.frac_coords) + t for site in sites]
            permutation = ["a" for i in range(len(sites))]
            not_found = list(range(len(sites)))
            for i in range(len(sites)):
                xformed = xformed_site_coords[i]
                for k, j in enumerate(not_found):
                    init = init_site_coords[j]
                    diff = xformed - init
                    if self._all_ints(diff, atol=atol):
                        permutation[i] = j
                        not_found.pop(k)
                        break

            same = np.zeros(len(sites))
            flipped = np.zeros(len(sites))
            for i, magmom in enumerate(xformed_magmoms):
                if (magmom == init_magmoms[permutation[i]]).all():
                    same[i] = 1
                elif (magmom == -1 * init_magmoms[permutation[i]]).all():
                    flipped[i] = 1

            if same.all():  # add symm op without tr
                mag_ops.append(
                    MagSymmOp.from_rotation_and_translation_and_time_reversal(
                        rotation_matrix=op.rotation_matrix,
                        translation_vec=op.translation_vector,
                        time_reversal=1,
                    )
                )
            if flipped.all():  # add symm op with tr
                mag_ops.append(
                    MagSymmOp.from_rotation_and_translation_and_time_reversal(
                        rotation_matrix=op.rotation_matrix,
                        translation_vec=op.translation_vector,
                        time_reversal=-1,
                    )
                )

        return mag_ops

    @staticmethod
    def _get_reciprocal_point_group(ops, R, A):
        Ainv = np.linalg.inv(A)
        # convert to reciprocal primitive basis
        recip_point_group = [np.around(np.dot(A, np.dot(R, Ainv)), decimals=2)]
        for op in ops:
            op = np.around(np.dot(A, np.dot(op, Ainv)), decimals=2)
            new = True
            new_coset = True
            for thing in recip_point_group:
                if (thing == op).all():
                    new = False
                if (thing == np.dot(R, op)).all():
                    new_coset = False

            if new:
                recip_point_group.append(op)
            if new_coset:
                recip_point_group.append(np.dot(R, op))

        return recip_point_group

    @staticmethod
    def _closewrapped(pos1, pos2, tolerance):
        pos1 = pos1 % 1.0
        pos2 = pos2 % 1.0

        if len(pos1) != len(pos2):
            return False
        for i, v in enumerate(pos1):
            if abs(pos1[i] - pos2[i]) > tolerance[i] and abs(pos1[i] - pos2[i]) < 1.0 - tolerance[i]:
                return False
        return True

    def _get_coset_factor(self, G, H):
        # finds g for left coset decomposition G = H + gH (H must be subgroup of G with index two.)
        # in this implementation, G and H are lists of objects of type
        # SymmOp
        gH = []
        for i, op1 in enumerate(G):
            in_H = False
            for op2 in H:
                if np.allclose(op1.rotation_matrix, op2.rotation_matrix, atol=self._atol) and self._closewrapped(
                    op1.translation_vector,
                    op2.translation_vector,
                    np.ones(3) * self._atol,
                ):
                    in_H = True
                    break
            if not in_H:
                gH.append(op1)

        for op in gH:
            opH = [op.__mul__(h) for h in H]
            is_coset_factor = True
            for op1 in opH:
                for op2 in H:
                    if np.allclose(op1.rotation_matrix, op2.rotation_matrix, atol=self._atol) and self._closewrapped(
                        op1.translation_vector,
                        op2.translation_vector,
                        np.ones(3) * self._atol,
                    ):

                        is_coset_factor = False
                        break
                if not is_coset_factor:
                    break
            if is_coset_factor:
                return op

        return "No coset factor found."

    @staticmethod
    def _apply_op_to_magmom(r, magmom):
        if np.linalg.det(r) == 1:
            return np.dot(r, magmom)
        return -1 * np.dot(r, magmom)

    @staticmethod
    def _all_ints(arr, atol):
        rounded_arr = np.around(arr, decimals=0)
        return np.allclose(rounded_arr, arr, atol=atol)

    def _get_IRBZ(self, recip_point_group, W, key_points, face_center_inds, atol):
        rpgdict = self._get_reciprocal_point_group_dict(recip_point_group, atol)

        g = np.dot(W.T, W)  # just using change of basis matrix rather than
        # Lattice.get_cartesian_coordinates for conciseness
        ginv = np.linalg.inv(g)
        D = np.linalg.det(W)

        primary_orientation = None
        secondary_orientation = None
        tertiary_orientation = None

        planar_boundaries = []
        IRBZ_points = list(enumerate(key_points))

        for sigma in rpgdict["reflections"]:
            norm = sigma["normal"]
            if primary_orientation is None:
                primary_orientation = norm
                planar_boundaries.append(norm)
            elif np.isclose(np.dot(primary_orientation, np.dot(g, norm)), 0, atol=atol):
                if secondary_orientation is None:
                    secondary_orientation = norm
                    planar_boundaries.append(norm)
                elif np.isclose(np.dot(secondary_orientation, np.dot(g, norm)), 0, atol=atol):
                    if tertiary_orientation is None:
                        tertiary_orientation = norm
                        planar_boundaries.append(norm)
                    elif np.allclose(norm, -1 * tertiary_orientation, atol=atol):
                        pass
                elif np.dot(secondary_orientation, np.dot(g, norm)) < 0:
                    planar_boundaries.append(-1 * norm)
                else:
                    planar_boundaries.append(norm)
            elif np.dot(primary_orientation, np.dot(g, norm)) < 0:
                planar_boundaries.append(-1 * norm)
            else:
                planar_boundaries.append(norm)

        IRBZ_points = self._reduce_IRBZ(IRBZ_points, planar_boundaries, g, atol)

        used_axes = []

        # six-fold rotoinversion always comes with horizontal mirror so don't
        # need to check
        for rotn in rpgdict["rotations"]["six-fold"]:
            ax = rotn["axis"]
            op = rotn["op"]
            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
                    face_center_found = False
                    for point in IRBZ_points:
                        if point[0] in face_center_inds:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                face_center_found = True
                                used_axes.append(ax)
                                break
                    if not face_center_found:
                        print("face center not found")
                        for point in IRBZ_points:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                used_axes.append(ax)
                                break
                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)

        for rotn in rpgdict["rotations"]["rotoinv-four-fold"]:
            ax = rotn["axis"]
            op = rotn["op"]
            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
                    face_center_found = False
                    for point in IRBZ_points:
                        if point[0] in face_center_inds:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, np.dot(op, cross)]
                                face_center_found = True
                                used_axes.append(ax)
                                break
                    if not face_center_found:
                        print("face center not found")
                        for point in IRBZ_points:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                used_axes.append(ax)
                                break
                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)

        for rotn in rpgdict["rotations"]["four-fold"]:
            ax = rotn["axis"]
            op = rotn["op"]
            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
                    face_center_found = False
                    for point in IRBZ_points:
                        if point[0] in face_center_inds:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                face_center_found = True
                                used_axes.append(ax)
                                break
                    if not face_center_found:
                        print("face center not found")
                        for point in IRBZ_points:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                used_axes.append(ax)
                                break
                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)

        for rotn in rpgdict["rotations"]["rotoinv-three-fold"]:
            ax = rotn["axis"]
            op = rotn["op"]
            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
                    face_center_found = False
                    for point in IRBZ_points:
                        if point[0] in face_center_inds:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [
                                    cross,
                                    -1 * np.dot(sqrtm(-1 * op), cross),
                                ]
                                face_center_found = True
                                used_axes.append(ax)
                                break
                    if not face_center_found:
                        print("face center not found")
                        for point in IRBZ_points:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                used_axes.append(ax)
                                break
                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)

        for rotn in rpgdict["rotations"]["three-fold"]:
            ax = rotn["axis"]
            op = rotn["op"]
            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
                    face_center_found = False
                    for point in IRBZ_points:
                        if point[0] in face_center_inds:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                face_center_found = True
                                used_axes.append(ax)
                                break
                    if not face_center_found:
                        print("face center not found")
                        for point in IRBZ_points:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                used_axes.append(ax)
                                break
                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)

        for rotn in rpgdict["rotations"]["two-fold"]:
            ax = rotn["axis"]
            op = rotn["op"]
            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
                    face_center_found = False
                    for point in IRBZ_points:
                        if point[0] in face_center_inds:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                face_center_found = True
                                used_axes.append(ax)
                                break
                    if not face_center_found:
                        print("face center not found")
                        for point in IRBZ_points:
                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
                            if not np.allclose(cross, 0, atol=atol):
                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
                                used_axes.append(ax)
                                break
                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)

        return [point[0] for point in IRBZ_points]

    @staticmethod
    def _get_reciprocal_point_group_dict(recip_point_group, atol):
        PAR = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, -1]])

        d = {
            "reflections": [],
            "rotations": {
                "two-fold": [],
                "three-fold": [],
                "four-fold": [],
                "six-fold": [],
                "rotoinv-three-fold": [],
                "rotoinv-four-fold": [],
                "rotoinv-six-fold": [],
            },
            "inversion": [],
        }

        for i, op in enumerate(recip_point_group):
            evals, evects = np.linalg.eig(op)

            tr = np.trace(op)
            det = np.linalg.det(op)

            # Proper rotations
            if np.isclose(det, 1, atol=atol):
                if np.isclose(tr, 3, atol=atol):
                    continue
                if np.isclose(tr, -1, atol=atol):  # two-fold rotation
                    for j in range(3):
                        if np.isclose(evals[j], 1, atol=atol):
                            ax = evects[:, j]
                    d["rotations"]["two-fold"].append({"ind": i, "axis": ax, "op": op})
                elif np.isclose(tr, 0, atol=atol):  # three-fold rotation
                    for j in range(3):
                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
                            ax = evects[:, j]
                    d["rotations"]["three-fold"].append({"ind": i, "axis": ax, "op": op})
                # four-fold rotation
                elif np.isclose(tr, 1, atol=atol):
                    for j in range(3):
                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
                            ax = evects[:, j]
                    d["rotations"]["four-fold"].append({"ind": i, "axis": ax, "op": op})
                elif np.isclose(tr, 2, atol=atol):  # six-fold rotation
                    for j in range(3):
                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
                            ax = evects[:, j]
                    d["rotations"]["six-fold"].append({"ind": i, "axis": ax, "op": op})

            # Improper rotations
            if np.isclose(det, -1, atol=atol):
                if np.isclose(tr, -3, atol=atol):
                    d["inversion"].append({"ind": i, "op": PAR})
                elif np.isclose(tr, 1, atol=atol):  # two-fold rotation
                    for j in range(3):
                        if np.isclose(evals[j], -1, atol=atol):
                            norm = evects[:, j]
                    d["reflections"].append({"ind": i, "normal": norm, "op": op})
                elif np.isclose(tr, 0, atol=atol):  # three-fold rotoinversion
                    for j in range(3):
                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
                            ax = evects[:, j]
                    d["rotations"]["rotoinv-three-fold"].append({"ind": i, "axis": ax, "op": op})
                # four-fold rotoinversion
                elif np.isclose(tr, -1, atol=atol):
                    for j in range(3):
                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
                            ax = evects[:, j]
                    d["rotations"]["rotoinv-four-fold"].append({"ind": i, "axis": ax, "op": op})
                # six-fold rotoinversion
                elif np.isclose(tr, -2, atol=atol):
                    for j in range(3):
                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
                            ax = evects[:, j]
                    d["rotations"]["rotoinv-six-fold"].append({"ind": i, "axis": ax, "op": op})

        return d

    @staticmethod
    def _op_maps_IRBZ_to_self(op, IRBZ_points, atol):
        point_coords = [point[1] for point in IRBZ_points]
        for point in point_coords:
            point_prime = np.dot(op, point)
            mapped_back = False
            for checkpoint in point_coords:
                if np.allclose(point_prime, checkpoint, atol):
                    mapped_back = True
                    break
            if not mapped_back:
                return False

        return True

    @staticmethod
    def _reduce_IRBZ(IRBZ_points, boundaries, g, atol):
        in_reduced_section = []
        for point in IRBZ_points:
            in_reduced_section.append(
                np.all(
                    [
                        (
                            np.dot(point[1], np.dot(g, boundary)) >= 0
                            or np.isclose(np.dot(point[1], np.dot(g, boundary)), 0, atol=atol)
                        )
                        for boundary in boundaries
                    ]
                )
            )

        return [IRBZ_points[i] for i in range(len(IRBZ_points)) if in_reduced_section[i]]

    def _get_orbit_labels(self, orbit_cosines_orig, key_points_inds_orbits, atol):

        orbit_cosines_copy = orbit_cosines_orig.copy()
        orbit_labels_unsorted = [(len(key_points_inds_orbits) - 1, 26)]
        orbit_inds_remaining = range(len(key_points_inds_orbits) - 1)
        pop_orbits = []
        pop_labels = []

        for i, orb_cos in enumerate(orbit_cosines_copy):
            if np.isclose(orb_cos[0][1], 1.0, atol=atol):
                # (point orbit index, label index)
                orbit_labels_unsorted.append((i, orb_cos[0][0]))
                pop_orbits.append(i)
                pop_labels.append(orb_cos[0][0])

        orbit_cosines_copy = self._reduce_cosines_array(orbit_cosines_copy, pop_orbits, pop_labels)
        orbit_inds_remaining = [i for i in orbit_inds_remaining if i not in pop_orbits]

        # orbit_labels_unsorted already contains gamma orbit
        while len(orbit_labels_unsorted) < len(orbit_cosines_orig) + 1:
            pop_orbits = []
            pop_labels = []
            max_cosine_value = max([orb_cos[0][1] for orb_cos in orbit_cosines_copy])
            max_cosine_value_inds = [
                j for j in range(len(orbit_cosines_copy)) if orbit_cosines_copy[j][0][1] == max_cosine_value
            ]
            max_cosine_label_inds = self._get_max_cosine_labels(
                [orbit_cosines_copy[j] for j in max_cosine_value_inds],
                key_points_inds_orbits,
                atol,
            )

            for j, label_ind in enumerate(max_cosine_label_inds):
                orbit_labels_unsorted.append((orbit_inds_remaining[max_cosine_value_inds[j]], label_ind))
                pop_orbits.append(max_cosine_value_inds[j])
                pop_labels.append(label_ind)
            orbit_cosines_copy = self._reduce_cosines_array(orbit_cosines_copy, pop_orbits, pop_labels)
            orbit_inds_remaining = [
                orbit_inds_remaining[j] for j in range(len(orbit_inds_remaining)) if j not in pop_orbits
            ]

        orbit_labels = np.zeros(len(key_points_inds_orbits))
        for tup in orbit_labels_unsorted:
            orbit_labels[tup[0]] = tup[1]

        return orbit_labels

    @staticmethod
    def _reduce_cosines_array(orbit_cosines, pop_orbits, pop_labels):
        return [
            [orb_cos[i] for i in range(len(orb_cos)) if orb_cos[i][0] not in pop_labels]
            for j, orb_cos in enumerate(orbit_cosines)
            if j not in pop_orbits
        ]

    def _get_max_cosine_labels(self, max_cosine_orbits_orig, key_points_inds_orbits, atol):
        max_cosine_orbits_copy = max_cosine_orbits_orig.copy()
        max_cosine_label_inds = np.zeros(len(max_cosine_orbits_copy))
        initial_max_cosine_label_inds = [max_cos_orb[0][0] for max_cos_orb in max_cosine_orbits_copy]
        u, inds, counts = np.unique(initial_max_cosine_label_inds, return_index=True, return_counts=True)
        grouped_inds = [
            [
                i
                for i in range(len(initial_max_cosine_label_inds))
                if max_cosine_orbits_copy[i][0][0] == max_cosine_orbits_copy[ind][0][0]
            ]
            for ind in inds
        ]
        pop_orbits = []
        pop_labels = []
        unassigned_orbits = []
        for i, ind in enumerate(inds):
            if counts[i] == 1:
                max_cosine_label_inds[ind] = initial_max_cosine_label_inds[ind]
                pop_orbits.append(ind)
                pop_labels.append(initial_max_cosine_label_inds[ind])
            else:
                next_choices = []
                for grouped_ind in grouped_inds[i]:
                    j = 1
                    while True:
                        if max_cosine_orbits_copy[grouped_ind][j][0] not in initial_max_cosine_label_inds:
                            next_choices.append(max_cosine_orbits_copy[grouped_ind][j][1])
                            break
                        j += 1
                worst_next_choice = next_choices.index(min(next_choices))
                for grouped_ind in grouped_inds[i]:
                    if grouped_ind != worst_next_choice:
                        unassigned_orbits.append(grouped_ind)
                max_cosine_label_inds[grouped_inds[i][worst_next_choice]] = initial_max_cosine_label_inds[
                    grouped_inds[i][worst_next_choice]
                ]
                pop_orbits.append(grouped_inds[i][worst_next_choice])
                pop_labels.append(initial_max_cosine_label_inds[grouped_inds[i][worst_next_choice]])

        if len(unassigned_orbits) != 0:
            max_cosine_orbits_copy = self._reduce_cosines_array(max_cosine_orbits_copy, pop_orbits, pop_labels)
            unassigned_orbits_labels = self._get_orbit_labels(max_cosine_orbits_copy, key_points_inds_orbits, atol)
            for i, unassigned_orbit in enumerate(unassigned_orbits):
                max_cosine_label_inds[unassigned_orbit] = unassigned_orbits_labels[i]

        return max_cosine_label_inds

    @staticmethod
    def LabelPoints(index):
        """
        Axes used in generating labels for Latimer-Munro convention
        """

        points = [
            [1, 0, 0],
            [0, 1, 0],
            [0, 0, 1],
            [1, 1, 0],
            [1, 0, 1],
            [0, 1, 1],
            [1, 1, 1],
            [1, 2, 0],
            [1, 0, 2],
            [1, 2, 2],
            [2, 1, 0],
            [0, 1, 2],
            [2, 1, 2],
            [2, 0, 1],
            [0, 2, 1],
            [2, 2, 1],
            [1, 1, 2],
            [1, 2, 1],
            [2, 1, 1],
            [3, 3, 2],
            [3, 2, 3],
            [2, 3, 3],
            [2, 2, 2],
            [3, 2, 2],
            [2, 3, 2],
            [1e-10, 1e-10, 1e-10],
        ]

        return points[index]

    @staticmethod
    def LabelSymbol(index):
        """
        Letters used in generating labels for Latimer-Munro convention
        """

        symbols = [
            "a",
            "b",
            "c",
            "d",
            "e",
            "f",
            "g",
            "h",
            "i",
            "j",
            "k",
            "l",
            "m",
            "n",
            "o",
            "p",
            "q",
            "r",
            "s",
            "t",
            "u",
            "v",
            "w",
            "x",
            "y",
            "z",
            "Γ",
        ]
        return symbols[index]