1# coding: utf-8
2# Copyright (c) Pymatgen Development Team.
3# Distributed under the terms of the MIT License.
4
5"""
6Provides classes for generating high-symmetry k-paths using different conventions.
7"""
8
9import abc
10import itertools
11import operator
12from math import ceil, cos, e, pi, sin, tan
13from warnings import warn
14
15import networkx as nx
16import numpy as np
17import spglib
18from monty.dev import requires
19from scipy.linalg import sqrtm
20
21from pymatgen.core.lattice import Lattice
22from pymatgen.core.operations import MagSymmOp, SymmOp
23from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
24
25try:
26    from seekpath import get_path  # type: ignore
27except ImportError:
28    get_path = None
29
30__author__ = "Geoffroy Hautier, Katherine Latimer, Jason Munro"
31__copyright__ = "Copyright 2020, The Materials Project"
32__version__ = "0.1"
33__maintainer__ = "Jason Munro"
34__email__ = "jmunro@lbl.gov"
35__status__ = "Development"
36__date__ = "March 2020"
37
38
39class KPathBase(metaclass=abc.ABCMeta):
40    """
41    This is the base class for classes used to generate high-symmetry
42    paths in reciprocal space (k-paths) for band structure calculations.
43    """
44
45    @abc.abstractmethod
46    def __init__(self, structure, symprec=0.01, angle_tolerance=5, atol=1e-5, *args, **kwargs):
47        """
48        Args:
49        structure (Structure): Structure object
50        symprec (float): Tolerance for symmetry finding
51        angle_tolerance (float): Angle tolerance for symmetry finding.
52        atol (float): Absolute tolerance used to compare structures
53            and determine symmetric equivalence of points and lines
54            in the BZ.
55        """
56
57        self._structure = structure
58        self._latt = self._structure.lattice
59        self._rec_lattice = self._structure.lattice.reciprocal_lattice
60        self._kpath = None
61
62        self._symprec = symprec
63        self._atol = atol
64        self._angle_tolerance = angle_tolerance
65
66    @property
67    def structure(self):
68        """
69        Returns:
70        The input structure
71        """
72        return self._structure
73
74    @property
75    def lattice(self):
76        """
77        Returns:
78        The real space lattice
79        """
80        return self._latt
81
82    @property
83    def rec_lattice(self):
84        """
85        Returns:
86        The reciprocal space lattice
87        """
88        return self._rec_lattice
89
90    @property
91    def kpath(self):
92        """
93        Returns:
94        The symmetry line path in reciprocal space
95        """
96        return self._kpath
97
98    def get_kpoints(self, line_density=20, coords_are_cartesian=True):
99        """
100        Returns:
101        the kpoints along the paths in cartesian coordinates
102        together with the labels for symmetry points -Wei.
103        """
104        list_k_points = []
105        sym_point_labels = []
106        for b in self.kpath["path"]:
107            for i in range(1, len(b)):
108                start = np.array(self.kpath["kpoints"][b[i - 1]])
109                end = np.array(self.kpath["kpoints"][b[i]])
110                distance = np.linalg.norm(
111                    self._rec_lattice.get_cartesian_coords(start) - self._rec_lattice.get_cartesian_coords(end)
112                )
113                nb = int(ceil(distance * line_density))
114                if nb == 0:
115                    continue
116                sym_point_labels.extend([b[i - 1]] + [""] * (nb - 1) + [b[i]])
117                list_k_points.extend(
118                    [
119                        self._rec_lattice.get_cartesian_coords(start)
120                        + float(i)
121                        / float(nb)
122                        * (self._rec_lattice.get_cartesian_coords(end) - self._rec_lattice.get_cartesian_coords(start))
123                        for i in range(0, nb + 1)
124                    ]
125                )
126        if coords_are_cartesian:
127            return list_k_points, sym_point_labels
128
129        frac_k_points = [self._rec_lattice.get_fractional_coords(k) for k in list_k_points]
130        return frac_k_points, sym_point_labels
131
132
133class KPathSetyawanCurtarolo(KPathBase):
134    """
135    This class looks for path along high symmetry lines in
136    the Brillouin Zone.
137    It is based on Setyawan, W., & Curtarolo, S. (2010).
138    High-throughput electronic band structure calculations:
139    Challenges and tools. Computational Materials Science,
140    49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010
141    It should be used with primitive structures that
142    comply with the definition from the paper.
143    The symmetry is determined by spglib through the
144    SpacegroupAnalyzer class. The analyzer can be used to
145    produce the correct primitive structure (method
146    get_primitive_standard_structure(international_monoclinic=False)).
147    A warning will signal possible compatibility problems
148    with the given structure. KPoints from get_kpoints() method
149    are returned in the reciprocal cell basis defined in the paper.
150    """
151
152    def __init__(self, structure, symprec=0.01, angle_tolerance=5, atol=1e-5):
153        """
154        Args:
155        structure (Structure): Structure object
156        symprec (float): Tolerance for symmetry finding
157        angle_tolerance (float): Angle tolerance for symmetry finding.
158        atol (float): Absolute tolerance used to compare the input
159            structure with the one expected as primitive standard.
160            A warning will be issued if the lattices don't match.
161        """
162        if "magmom" in structure.site_properties.keys():
163            warn(
164                "'magmom' entry found in site properties but will be ignored \
165                  for the Setyawan and Curtarolo convention."
166            )
167
168        super().__init__(structure, symprec=symprec, angle_tolerance=angle_tolerance, atol=atol)
169
170        self._sym = SpacegroupAnalyzer(structure, symprec=symprec, angle_tolerance=angle_tolerance)
171        self._prim = self._sym.get_primitive_standard_structure(international_monoclinic=False)
172        self._conv = self._sym.get_conventional_standard_structure(international_monoclinic=False)
173        self._rec_lattice = self._prim.lattice.reciprocal_lattice
174
175        # Note: this warning will be issued for space groups 38-41, since the primitive cell must be
176        # reformatted to match Setyawan/Curtarolo convention in order to work with the current k-path
177        # generation scheme.
178        if not np.allclose(self._structure.lattice.matrix, self._prim.lattice.matrix, atol=atol):
179            warn(
180                "The input structure does not match the expected standard primitive! "
181                "The path can be incorrect. Use at your own risk."
182            )
183
184        lattice_type = self._sym.get_lattice_type()
185        spg_symbol = self._sym.get_space_group_symbol()
186
187        if lattice_type == "cubic":
188            if "P" in spg_symbol:
189                self._kpath = self.cubic()
190            elif "F" in spg_symbol:
191                self._kpath = self.fcc()
192            elif "I" in spg_symbol:
193                self._kpath = self.bcc()
194            else:
195                warn("Unexpected value for spg_symbol: %s" % spg_symbol)
196
197        elif lattice_type == "tetragonal":
198            if "P" in spg_symbol:
199                self._kpath = self.tet()
200            elif "I" in spg_symbol:
201                a = self._conv.lattice.abc[0]
202                c = self._conv.lattice.abc[2]
203                if c < a:
204                    self._kpath = self.bctet1(c, a)
205                else:
206                    self._kpath = self.bctet2(c, a)
207            else:
208                warn("Unexpected value for spg_symbol: %s" % spg_symbol)
209
210        elif lattice_type == "orthorhombic":
211            a = self._conv.lattice.abc[0]
212            b = self._conv.lattice.abc[1]
213            c = self._conv.lattice.abc[2]
214
215            if "P" in spg_symbol:
216                self._kpath = self.orc()
217
218            elif "F" in spg_symbol:
219                if 1 / a ** 2 > 1 / b ** 2 + 1 / c ** 2:
220                    self._kpath = self.orcf1(a, b, c)
221                elif 1 / a ** 2 < 1 / b ** 2 + 1 / c ** 2:
222                    self._kpath = self.orcf2(a, b, c)
223                else:
224                    self._kpath = self.orcf3(a, b, c)
225
226            elif "I" in spg_symbol:
227                self._kpath = self.orci(a, b, c)
228
229            elif "C" in spg_symbol or "A" in spg_symbol:
230                self._kpath = self.orcc(a, b, c)
231            else:
232                warn("Unexpected value for spg_symbol: %s" % spg_symbol)
233
234        elif lattice_type == "hexagonal":
235            self._kpath = self.hex()
236
237        elif lattice_type == "rhombohedral":
238            alpha = self._prim.lattice.parameters[3]
239            if alpha < 90:
240                self._kpath = self.rhl1(alpha * pi / 180)
241            else:
242                self._kpath = self.rhl2(alpha * pi / 180)
243
244        elif lattice_type == "monoclinic":
245            a, b, c = self._conv.lattice.abc
246            alpha = self._conv.lattice.parameters[3]
247            # beta = self._conv.lattice.parameters[4]
248
249            if "P" in spg_symbol:
250                self._kpath = self.mcl(b, c, alpha * pi / 180)
251
252            elif "C" in spg_symbol:
253                kgamma = self._rec_lattice.parameters[5]
254                if kgamma > 90:
255                    self._kpath = self.mclc1(a, b, c, alpha * pi / 180)
256                if kgamma == 90:
257                    self._kpath = self.mclc2(a, b, c, alpha * pi / 180)
258                if kgamma < 90:
259                    if b * cos(alpha * pi / 180) / c + b ** 2 * sin(alpha * pi / 180) ** 2 / a ** 2 < 1:
260                        self._kpath = self.mclc3(a, b, c, alpha * pi / 180)
261                    if b * cos(alpha * pi / 180) / c + b ** 2 * sin(alpha * pi / 180) ** 2 / a ** 2 == 1:
262                        self._kpath = self.mclc4(a, b, c, alpha * pi / 180)
263                    if b * cos(alpha * pi / 180) / c + b ** 2 * sin(alpha * pi / 180) ** 2 / a ** 2 > 1:
264                        self._kpath = self.mclc5(a, b, c, alpha * pi / 180)
265            else:
266                warn("Unexpected value for spg_symbol: %s" % spg_symbol)
267
268        elif lattice_type == "triclinic":
269            kalpha = self._rec_lattice.parameters[3]
270            kbeta = self._rec_lattice.parameters[4]
271            kgamma = self._rec_lattice.parameters[5]
272            if kalpha > 90 and kbeta > 90 and kgamma > 90:
273                self._kpath = self.tria()
274            if kalpha < 90 and kbeta < 90 and kgamma < 90:
275                self._kpath = self.trib()
276            if kalpha > 90 and kbeta > 90 and kgamma == 90:
277                self._kpath = self.tria()
278            if kalpha < 90 and kbeta < 90 and kgamma == 90:
279                self._kpath = self.trib()
280
281        else:
282            warn("Unknown lattice type %s" % lattice_type)
283
284    @property
285    def conventional(self):
286        """
287        Returns:
288            The conventional cell structure
289        """
290        return self._conv
291
292    @property
293    def prim(self):
294        """
295        Returns:
296            The primitive cell structure
297        """
298        return self._prim
299
300    @property
301    def prim_rec(self):
302        """
303        Returns:
304            The primitive reciprocal cell structure
305        """
306        return self._rec_lattice
307
308    def cubic(self):
309        """
310        CUB Path
311        """
312        self.name = "CUB"
313        kpoints = {
314            "\\Gamma": np.array([0.0, 0.0, 0.0]),
315            "X": np.array([0.0, 0.5, 0.0]),
316            "R": np.array([0.5, 0.5, 0.5]),
317            "M": np.array([0.5, 0.5, 0.0]),
318        }
319        path = [["\\Gamma", "X", "M", "\\Gamma", "R", "X"], ["M", "R"]]
320        return {"kpoints": kpoints, "path": path}
321
322    def fcc(self):
323        """
324        FCC Path
325        """
326        self.name = "FCC"
327        kpoints = {
328            "\\Gamma": np.array([0.0, 0.0, 0.0]),
329            "K": np.array([3.0 / 8.0, 3.0 / 8.0, 3.0 / 4.0]),
330            "L": np.array([0.5, 0.5, 0.5]),
331            "U": np.array([5.0 / 8.0, 1.0 / 4.0, 5.0 / 8.0]),
332            "W": np.array([0.5, 1.0 / 4.0, 3.0 / 4.0]),
333            "X": np.array([0.5, 0.0, 0.5]),
334        }
335        path = [
336            ["\\Gamma", "X", "W", "K", "\\Gamma", "L", "U", "W", "L", "K"],
337            ["U", "X"],
338        ]
339        return {"kpoints": kpoints, "path": path}
340
341    def bcc(self):
342        """
343        BCC Path
344        """
345        self.name = "BCC"
346        kpoints = {
347            "\\Gamma": np.array([0.0, 0.0, 0.0]),
348            "H": np.array([0.5, -0.5, 0.5]),
349            "P": np.array([0.25, 0.25, 0.25]),
350            "N": np.array([0.0, 0.0, 0.5]),
351        }
352        path = [["\\Gamma", "H", "N", "\\Gamma", "P", "H"], ["P", "N"]]
353        return {"kpoints": kpoints, "path": path}
354
355    def tet(self):
356        """
357        TET Path
358        """
359        self.name = "TET"
360        kpoints = {
361            "\\Gamma": np.array([0.0, 0.0, 0.0]),
362            "A": np.array([0.5, 0.5, 0.5]),
363            "M": np.array([0.5, 0.5, 0.0]),
364            "R": np.array([0.0, 0.5, 0.5]),
365            "X": np.array([0.0, 0.5, 0.0]),
366            "Z": np.array([0.0, 0.0, 0.5]),
367        }
368        path = [
369            ["\\Gamma", "X", "M", "\\Gamma", "Z", "R", "A", "Z"],
370            ["X", "R"],
371            ["M", "A"],
372        ]
373        return {"kpoints": kpoints, "path": path}
374
375    def bctet1(self, c, a):
376        """
377        BCT1 Path
378        """
379        self.name = "BCT1"
380        eta = (1 + c ** 2 / a ** 2) / 4.0
381        kpoints = {
382            "\\Gamma": np.array([0.0, 0.0, 0.0]),
383            "M": np.array([-0.5, 0.5, 0.5]),
384            "N": np.array([0.0, 0.5, 0.0]),
385            "P": np.array([0.25, 0.25, 0.25]),
386            "X": np.array([0.0, 0.0, 0.5]),
387            "Z": np.array([eta, eta, -eta]),
388            "Z_1": np.array([-eta, 1 - eta, eta]),
389        }
390        path = [["\\Gamma", "X", "M", "\\Gamma", "Z", "P", "N", "Z_1", "M"], ["X", "P"]]
391        return {"kpoints": kpoints, "path": path}
392
393    def bctet2(self, c, a):
394        """
395        BCT2 Path
396        """
397        self.name = "BCT2"
398        eta = (1 + a ** 2 / c ** 2) / 4.0
399        zeta = a ** 2 / (2 * c ** 2)
400        kpoints = {
401            "\\Gamma": np.array([0.0, 0.0, 0.0]),
402            "N": np.array([0.0, 0.5, 0.0]),
403            "P": np.array([0.25, 0.25, 0.25]),
404            "\\Sigma": np.array([-eta, eta, eta]),
405            "\\Sigma_1": np.array([eta, 1 - eta, -eta]),
406            "X": np.array([0.0, 0.0, 0.5]),
407            "Y": np.array([-zeta, zeta, 0.5]),
408            "Y_1": np.array([0.5, 0.5, -zeta]),
409            "Z": np.array([0.5, 0.5, -0.5]),
410        }
411        path = [
412            [
413                "\\Gamma",
414                "X",
415                "Y",
416                "\\Sigma",
417                "\\Gamma",
418                "Z",
419                "\\Sigma_1",
420                "N",
421                "P",
422                "Y_1",
423                "Z",
424            ],
425            ["X", "P"],
426        ]
427        return {"kpoints": kpoints, "path": path}
428
429    def orc(self):
430        """
431        ORC Path
432        """
433        self.name = "ORC"
434        kpoints = {
435            "\\Gamma": np.array([0.0, 0.0, 0.0]),
436            "R": np.array([0.5, 0.5, 0.5]),
437            "S": np.array([0.5, 0.5, 0.0]),
438            "T": np.array([0.0, 0.5, 0.5]),
439            "U": np.array([0.5, 0.0, 0.5]),
440            "X": np.array([0.5, 0.0, 0.0]),
441            "Y": np.array([0.0, 0.5, 0.0]),
442            "Z": np.array([0.0, 0.0, 0.5]),
443        }
444        path = [
445            ["\\Gamma", "X", "S", "Y", "\\Gamma", "Z", "U", "R", "T", "Z"],
446            ["Y", "T"],
447            ["U", "X"],
448            ["S", "R"],
449        ]
450        return {"kpoints": kpoints, "path": path}
451
452    def orcf1(self, a, b, c):
453        """
454        ORFC1 Path
455        """
456        self.name = "ORCF1"
457        zeta = (1 + a ** 2 / b ** 2 - a ** 2 / c ** 2) / 4
458        eta = (1 + a ** 2 / b ** 2 + a ** 2 / c ** 2) / 4
459
460        kpoints = {
461            "\\Gamma": np.array([0.0, 0.0, 0.0]),
462            "A": np.array([0.5, 0.5 + zeta, zeta]),
463            "A_1": np.array([0.5, 0.5 - zeta, 1 - zeta]),
464            "L": np.array([0.5, 0.5, 0.5]),
465            "T": np.array([1, 0.5, 0.5]),
466            "X": np.array([0.0, eta, eta]),
467            "X_1": np.array([1, 1 - eta, 1 - eta]),
468            "Y": np.array([0.5, 0.0, 0.5]),
469            "Z": np.array([0.5, 0.5, 0.0]),
470        }
471        path = [
472            ["\\Gamma", "Y", "T", "Z", "\\Gamma", "X", "A_1", "Y"],
473            ["T", "X_1"],
474            ["X", "A", "Z"],
475            ["L", "\\Gamma"],
476        ]
477        return {"kpoints": kpoints, "path": path}
478
479    def orcf2(self, a, b, c):
480        """
481        ORFC2 Path
482        """
483        self.name = "ORCF2"
484        phi = (1 + c ** 2 / b ** 2 - c ** 2 / a ** 2) / 4
485        eta = (1 + a ** 2 / b ** 2 - a ** 2 / c ** 2) / 4
486        delta = (1 + b ** 2 / a ** 2 - b ** 2 / c ** 2) / 4
487        kpoints = {
488            "\\Gamma": np.array([0.0, 0.0, 0.0]),
489            "C": np.array([0.5, 0.5 - eta, 1 - eta]),
490            "C_1": np.array([0.5, 0.5 + eta, eta]),
491            "D": np.array([0.5 - delta, 0.5, 1 - delta]),
492            "D_1": np.array([0.5 + delta, 0.5, delta]),
493            "L": np.array([0.5, 0.5, 0.5]),
494            "H": np.array([1 - phi, 0.5 - phi, 0.5]),
495            "H_1": np.array([phi, 0.5 + phi, 0.5]),
496            "X": np.array([0.0, 0.5, 0.5]),
497            "Y": np.array([0.5, 0.0, 0.5]),
498            "Z": np.array([0.5, 0.5, 0.0]),
499        }
500        path = [
501            ["\\Gamma", "Y", "C", "D", "X", "\\Gamma", "Z", "D_1", "H", "C"],
502            ["C_1", "Z"],
503            ["X", "H_1"],
504            ["H", "Y"],
505            ["L", "\\Gamma"],
506        ]
507        return {"kpoints": kpoints, "path": path}
508
509    def orcf3(self, a, b, c):
510        """
511        ORFC3 Path
512        """
513        self.name = "ORCF3"
514        zeta = (1 + a ** 2 / b ** 2 - a ** 2 / c ** 2) / 4
515        eta = (1 + a ** 2 / b ** 2 + a ** 2 / c ** 2) / 4
516        kpoints = {
517            "\\Gamma": np.array([0.0, 0.0, 0.0]),
518            "A": np.array([0.5, 0.5 + zeta, zeta]),
519            "A_1": np.array([0.5, 0.5 - zeta, 1 - zeta]),
520            "L": np.array([0.5, 0.5, 0.5]),
521            "T": np.array([1, 0.5, 0.5]),
522            "X": np.array([0.0, eta, eta]),
523            "X_1": np.array([1, 1 - eta, 1 - eta]),
524            "Y": np.array([0.5, 0.0, 0.5]),
525            "Z": np.array([0.5, 0.5, 0.0]),
526        }
527        path = [
528            ["\\Gamma", "Y", "T", "Z", "\\Gamma", "X", "A_1", "Y"],
529            ["X", "A", "Z"],
530            ["L", "\\Gamma"],
531        ]
532        return {"kpoints": kpoints, "path": path}
533
534    def orci(self, a, b, c):
535        """
536        ORCI Path
537        """
538        self.name = "ORCI"
539        zeta = (1 + a ** 2 / c ** 2) / 4
540        eta = (1 + b ** 2 / c ** 2) / 4
541        delta = (b ** 2 - a ** 2) / (4 * c ** 2)
542        mu = (a ** 2 + b ** 2) / (4 * c ** 2)
543        kpoints = {
544            "\\Gamma": np.array([0.0, 0.0, 0.0]),
545            "L": np.array([-mu, mu, 0.5 - delta]),
546            "L_1": np.array([mu, -mu, 0.5 + delta]),
547            "L_2": np.array([0.5 - delta, 0.5 + delta, -mu]),
548            "R": np.array([0.0, 0.5, 0.0]),
549            "S": np.array([0.5, 0.0, 0.0]),
550            "T": np.array([0.0, 0.0, 0.5]),
551            "W": np.array([0.25, 0.25, 0.25]),
552            "X": np.array([-zeta, zeta, zeta]),
553            "X_1": np.array([zeta, 1 - zeta, -zeta]),
554            "Y": np.array([eta, -eta, eta]),
555            "Y_1": np.array([1 - eta, eta, -eta]),
556            "Z": np.array([0.5, 0.5, -0.5]),
557        }
558        path = [
559            ["\\Gamma", "X", "L", "T", "W", "R", "X_1", "Z", "\\Gamma", "Y", "S", "W"],
560            ["L_1", "Y"],
561            ["Y_1", "Z"],
562        ]
563        return {"kpoints": kpoints, "path": path}
564
565    def orcc(self, a, b, c):
566        """
567        ORCC Path
568        """
569        self.name = "ORCC"
570        zeta = (1 + a ** 2 / b ** 2) / 4
571        kpoints = {
572            "\\Gamma": np.array([0.0, 0.0, 0.0]),
573            "A": np.array([zeta, zeta, 0.5]),
574            "A_1": np.array([-zeta, 1 - zeta, 0.5]),
575            "R": np.array([0.0, 0.5, 0.5]),
576            "S": np.array([0.0, 0.5, 0.0]),
577            "T": np.array([-0.5, 0.5, 0.5]),
578            "X": np.array([zeta, zeta, 0.0]),
579            "X_1": np.array([-zeta, 1 - zeta, 0.0]),
580            "Y": np.array([-0.5, 0.5, 0]),
581            "Z": np.array([0.0, 0.0, 0.5]),
582        }
583        path = [
584            [
585                "\\Gamma",
586                "X",
587                "S",
588                "R",
589                "A",
590                "Z",
591                "\\Gamma",
592                "Y",
593                "X_1",
594                "A_1",
595                "T",
596                "Y",
597            ],
598            ["Z", "T"],
599        ]
600        return {"kpoints": kpoints, "path": path}
601
602    def hex(self):
603        """
604        HEX Path
605        """
606        self.name = "HEX"
607        kpoints = {
608            "\\Gamma": np.array([0.0, 0.0, 0.0]),
609            "A": np.array([0.0, 0.0, 0.5]),
610            "H": np.array([1.0 / 3.0, 1.0 / 3.0, 0.5]),
611            "K": np.array([1.0 / 3.0, 1.0 / 3.0, 0.0]),
612            "L": np.array([0.5, 0.0, 0.5]),
613            "M": np.array([0.5, 0.0, 0.0]),
614        }
615        path = [
616            ["\\Gamma", "M", "K", "\\Gamma", "A", "L", "H", "A"],
617            ["L", "M"],
618            ["K", "H"],
619        ]
620        return {"kpoints": kpoints, "path": path}
621
622    def rhl1(self, alpha):
623        """
624        RHL1 Path
625        """
626        self.name = "RHL1"
627        eta = (1 + 4 * cos(alpha)) / (2 + 4 * cos(alpha))
628        nu = 3.0 / 4.0 - eta / 2.0
629        kpoints = {
630            "\\Gamma": np.array([0.0, 0.0, 0.0]),
631            "B": np.array([eta, 0.5, 1.0 - eta]),
632            "B_1": np.array([1.0 / 2.0, 1.0 - eta, eta - 1.0]),
633            "F": np.array([0.5, 0.5, 0.0]),
634            "L": np.array([0.5, 0.0, 0.0]),
635            "L_1": np.array([0.0, 0.0, -0.5]),
636            "P": np.array([eta, nu, nu]),
637            "P_1": np.array([1.0 - nu, 1.0 - nu, 1.0 - eta]),
638            "P_2": np.array([nu, nu, eta - 1.0]),
639            "Q": np.array([1.0 - nu, nu, 0.0]),
640            "X": np.array([nu, 0.0, -nu]),
641            "Z": np.array([0.5, 0.5, 0.5]),
642        }
643        path = [
644            ["\\Gamma", "L", "B_1"],
645            ["B", "Z", "\\Gamma", "X"],
646            ["Q", "F", "P_1", "Z"],
647            ["L", "P"],
648        ]
649        return {"kpoints": kpoints, "path": path}
650
651    def rhl2(self, alpha):
652        """
653        RHL2 Path
654        """
655        self.name = "RHL2"
656        eta = 1 / (2 * tan(alpha / 2.0) ** 2)
657        nu = 3.0 / 4.0 - eta / 2.0
658        kpoints = {
659            "\\Gamma": np.array([0.0, 0.0, 0.0]),
660            "F": np.array([0.5, -0.5, 0.0]),
661            "L": np.array([0.5, 0.0, 0.0]),
662            "P": np.array([1 - nu, -nu, 1 - nu]),
663            "P_1": np.array([nu, nu - 1.0, nu - 1.0]),
664            "Q": np.array([eta, eta, eta]),
665            "Q_1": np.array([1.0 - eta, -eta, -eta]),
666            "Z": np.array([0.5, -0.5, 0.5]),
667        }
668        path = [["\\Gamma", "P", "Z", "Q", "\\Gamma", "F", "P_1", "Q_1", "L", "Z"]]
669        return {"kpoints": kpoints, "path": path}
670
671    def mcl(self, b, c, beta):
672        """
673        MCL Path
674        """
675        self.name = "MCL"
676        eta = (1 - b * cos(beta) / c) / (2 * sin(beta) ** 2)
677        nu = 0.5 - eta * c * cos(beta) / b
678        kpoints = {
679            "\\Gamma": np.array([0.0, 0.0, 0.0]),
680            "A": np.array([0.5, 0.5, 0.0]),
681            "C": np.array([0.0, 0.5, 0.5]),
682            "D": np.array([0.5, 0.0, 0.5]),
683            "D_1": np.array([0.5, 0.5, -0.5]),
684            "E": np.array([0.5, 0.5, 0.5]),
685            "H": np.array([0.0, eta, 1.0 - nu]),
686            "H_1": np.array([0.0, 1.0 - eta, nu]),
687            "H_2": np.array([0.0, eta, -nu]),
688            "M": np.array([0.5, eta, 1.0 - nu]),
689            "M_1": np.array([0.5, 1 - eta, nu]),
690            "M_2": np.array([0.5, 1 - eta, nu]),
691            "X": np.array([0.0, 0.5, 0.0]),
692            "Y": np.array([0.0, 0.0, 0.5]),
693            "Y_1": np.array([0.0, 0.0, -0.5]),
694            "Z": np.array([0.5, 0.0, 0.0]),
695        }
696        path = [
697            ["\\Gamma", "Y", "H", "C", "E", "M_1", "A", "X", "H_1"],
698            ["M", "D", "Z"],
699            ["Y", "D"],
700        ]
701        return {"kpoints": kpoints, "path": path}
702
703    def mclc1(self, a, b, c, alpha):
704        """
705        MCLC1 Path
706        """
707        self.name = "MCLC1"
708        zeta = (2 - b * cos(alpha) / c) / (4 * sin(alpha) ** 2)
709        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
710        psi = 0.75 - a ** 2 / (4 * b ** 2 * sin(alpha) ** 2)
711        phi = psi + (0.75 - psi) * b * cos(alpha) / c
712        kpoints = {
713            "\\Gamma": np.array([0.0, 0.0, 0.0]),
714            "N": np.array([0.5, 0.0, 0.0]),
715            "N_1": np.array([0.0, -0.5, 0.0]),
716            "F": np.array([1 - zeta, 1 - zeta, 1 - eta]),
717            "F_1": np.array([zeta, zeta, eta]),
718            "F_2": np.array([-zeta, -zeta, 1 - eta]),
719            "I": np.array([phi, 1 - phi, 0.5]),
720            "I_1": np.array([1 - phi, phi - 1, 0.5]),
721            "L": np.array([0.5, 0.5, 0.5]),
722            "M": np.array([0.5, 0.0, 0.5]),
723            "X": np.array([1 - psi, psi - 1, 0.0]),
724            "X_1": np.array([psi, 1 - psi, 0.0]),
725            "X_2": np.array([psi - 1, -psi, 0.0]),
726            "Y": np.array([0.5, 0.5, 0.0]),
727            "Y_1": np.array([-0.5, -0.5, 0.0]),
728            "Z": np.array([0.0, 0.0, 0.5]),
729        }
730        path = [
731            ["\\Gamma", "Y", "F", "L", "I"],
732            ["I_1", "Z", "F_1"],
733            ["Y", "X_1"],
734            ["X", "\\Gamma", "N"],
735            ["M", "\\Gamma"],
736        ]
737        return {"kpoints": kpoints, "path": path}
738
739    def mclc2(self, a, b, c, alpha):
740        """
741        MCLC2 Path
742        """
743        self.name = "MCLC2"
744        zeta = (2 - b * cos(alpha) / c) / (4 * sin(alpha) ** 2)
745        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
746        psi = 0.75 - a ** 2 / (4 * b ** 2 * sin(alpha) ** 2)
747        phi = psi + (0.75 - psi) * b * cos(alpha) / c
748        kpoints = {
749            "\\Gamma": np.array([0.0, 0.0, 0.0]),
750            "N": np.array([0.5, 0.0, 0.0]),
751            "N_1": np.array([0.0, -0.5, 0.0]),
752            "F": np.array([1 - zeta, 1 - zeta, 1 - eta]),
753            "F_1": np.array([zeta, zeta, eta]),
754            "F_2": np.array([-zeta, -zeta, 1 - eta]),
755            "F_3": np.array([1 - zeta, -zeta, 1 - eta]),
756            "I": np.array([phi, 1 - phi, 0.5]),
757            "I_1": np.array([1 - phi, phi - 1, 0.5]),
758            "L": np.array([0.5, 0.5, 0.5]),
759            "M": np.array([0.5, 0.0, 0.5]),
760            "X": np.array([1 - psi, psi - 1, 0.0]),
761            "X_1": np.array([psi, 1 - psi, 0.0]),
762            "X_2": np.array([psi - 1, -psi, 0.0]),
763            "Y": np.array([0.5, 0.5, 0.0]),
764            "Y_1": np.array([-0.5, -0.5, 0.0]),
765            "Z": np.array([0.0, 0.0, 0.5]),
766        }
767        path = [
768            ["\\Gamma", "Y", "F", "L", "I"],
769            ["I_1", "Z", "F_1"],
770            ["N", "\\Gamma", "M"],
771        ]
772        return {"kpoints": kpoints, "path": path}
773
774    def mclc3(self, a, b, c, alpha):
775        """
776        MCLC3 Path
777        """
778        self.name = "MCLC3"
779        mu = (1 + b ** 2 / a ** 2) / 4.0
780        delta = b * c * cos(alpha) / (2 * a ** 2)
781        zeta = mu - 0.25 + (1 - b * cos(alpha) / c) / (4 * sin(alpha) ** 2)
782        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
783        phi = 1 + zeta - 2 * mu
784        psi = eta - 2 * delta
785        kpoints = {
786            "\\Gamma": np.array([0.0, 0.0, 0.0]),
787            "F": np.array([1 - phi, 1 - phi, 1 - psi]),
788            "F_1": np.array([phi, phi - 1, psi]),
789            "F_2": np.array([1 - phi, -phi, 1 - psi]),
790            "H": np.array([zeta, zeta, eta]),
791            "H_1": np.array([1 - zeta, -zeta, 1 - eta]),
792            "H_2": np.array([-zeta, -zeta, 1 - eta]),
793            "I": np.array([0.5, -0.5, 0.5]),
794            "M": np.array([0.5, 0.0, 0.5]),
795            "N": np.array([0.5, 0.0, 0.0]),
796            "N_1": np.array([0.0, -0.5, 0.0]),
797            "X": np.array([0.5, -0.5, 0.0]),
798            "Y": np.array([mu, mu, delta]),
799            "Y_1": np.array([1 - mu, -mu, -delta]),
800            "Y_2": np.array([-mu, -mu, -delta]),
801            "Y_3": np.array([mu, mu - 1, delta]),
802            "Z": np.array([0.0, 0.0, 0.5]),
803        }
804        path = [
805            ["\\Gamma", "Y", "F", "H", "Z", "I", "F_1"],
806            ["H_1", "Y_1", "X", "\\Gamma", "N"],
807            ["M", "\\Gamma"],
808        ]
809        return {"kpoints": kpoints, "path": path}
810
811    def mclc4(self, a, b, c, alpha):
812        """
813        MCLC4 Path
814        """
815        self.name = "MCLC4"
816        mu = (1 + b ** 2 / a ** 2) / 4.0
817        delta = b * c * cos(alpha) / (2 * a ** 2)
818        zeta = mu - 0.25 + (1 - b * cos(alpha) / c) / (4 * sin(alpha) ** 2)
819        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
820        phi = 1 + zeta - 2 * mu
821        psi = eta - 2 * delta
822        kpoints = {
823            "\\Gamma": np.array([0.0, 0.0, 0.0]),
824            "F": np.array([1 - phi, 1 - phi, 1 - psi]),
825            "F_1": np.array([phi, phi - 1, psi]),
826            "F_2": np.array([1 - phi, -phi, 1 - psi]),
827            "H": np.array([zeta, zeta, eta]),
828            "H_1": np.array([1 - zeta, -zeta, 1 - eta]),
829            "H_2": np.array([-zeta, -zeta, 1 - eta]),
830            "I": np.array([0.5, -0.5, 0.5]),
831            "M": np.array([0.5, 0.0, 0.5]),
832            "N": np.array([0.5, 0.0, 0.0]),
833            "N_1": np.array([0.0, -0.5, 0.0]),
834            "X": np.array([0.5, -0.5, 0.0]),
835            "Y": np.array([mu, mu, delta]),
836            "Y_1": np.array([1 - mu, -mu, -delta]),
837            "Y_2": np.array([-mu, -mu, -delta]),
838            "Y_3": np.array([mu, mu - 1, delta]),
839            "Z": np.array([0.0, 0.0, 0.5]),
840        }
841        path = [
842            ["\\Gamma", "Y", "F", "H", "Z", "I"],
843            ["H_1", "Y_1", "X", "\\Gamma", "N"],
844            ["M", "\\Gamma"],
845        ]
846        return {"kpoints": kpoints, "path": path}
847
848    def mclc5(self, a, b, c, alpha):
849        """
850        MCLC5 Path
851        """
852        self.name = "MCLC5"
853        zeta = (b ** 2 / a ** 2 + (1 - b * cos(alpha) / c) / sin(alpha) ** 2) / 4
854        eta = 0.5 + 2 * zeta * c * cos(alpha) / b
855        mu = eta / 2 + b ** 2 / (4 * a ** 2) - b * c * cos(alpha) / (2 * a ** 2)
856        nu = 2 * mu - zeta
857        rho = 1 - zeta * a ** 2 / b ** 2
858        omega = (4 * nu - 1 - b ** 2 * sin(alpha) ** 2 / a ** 2) * c / (2 * b * cos(alpha))
859        delta = zeta * c * cos(alpha) / b + omega / 2 - 0.25
860        kpoints = {
861            "\\Gamma": np.array([0.0, 0.0, 0.0]),
862            "F": np.array([nu, nu, omega]),
863            "F_1": np.array([1 - nu, 1 - nu, 1 - omega]),
864            "F_2": np.array([nu, nu - 1, omega]),
865            "H": np.array([zeta, zeta, eta]),
866            "H_1": np.array([1 - zeta, -zeta, 1 - eta]),
867            "H_2": np.array([-zeta, -zeta, 1 - eta]),
868            "I": np.array([rho, 1 - rho, 0.5]),
869            "I_1": np.array([1 - rho, rho - 1, 0.5]),
870            "L": np.array([0.5, 0.5, 0.5]),
871            "M": np.array([0.5, 0.0, 0.5]),
872            "N": np.array([0.5, 0.0, 0.0]),
873            "N_1": np.array([0.0, -0.5, 0.0]),
874            "X": np.array([0.5, -0.5, 0.0]),
875            "Y": np.array([mu, mu, delta]),
876            "Y_1": np.array([1 - mu, -mu, -delta]),
877            "Y_2": np.array([-mu, -mu, -delta]),
878            "Y_3": np.array([mu, mu - 1, delta]),
879            "Z": np.array([0.0, 0.0, 0.5]),
880        }
881        path = [
882            ["\\Gamma", "Y", "F", "L", "I"],
883            ["I_1", "Z", "H", "F_1"],
884            ["H_1", "Y_1", "X", "\\Gamma", "N"],
885            ["M", "\\Gamma"],
886        ]
887        return {"kpoints": kpoints, "path": path}
888
889    def tria(self):
890        """
891        TRI1a Path
892        """
893        self.name = "TRI1a"
894        kpoints = {
895            "\\Gamma": np.array([0.0, 0.0, 0.0]),
896            "L": np.array([0.5, 0.5, 0.0]),
897            "M": np.array([0.0, 0.5, 0.5]),
898            "N": np.array([0.5, 0.0, 0.5]),
899            "R": np.array([0.5, 0.5, 0.5]),
900            "X": np.array([0.5, 0.0, 0.0]),
901            "Y": np.array([0.0, 0.5, 0.0]),
902            "Z": np.array([0.0, 0.0, 0.5]),
903        }
904        path = [
905            ["X", "\\Gamma", "Y"],
906            ["L", "\\Gamma", "Z"],
907            ["N", "\\Gamma", "M"],
908            ["R", "\\Gamma"],
909        ]
910        return {"kpoints": kpoints, "path": path}
911
912    def trib(self):
913        """
914        TRI1b Path
915        """
916        self.name = "TRI1b"
917        kpoints = {
918            "\\Gamma": np.array([0.0, 0.0, 0.0]),
919            "L": np.array([0.5, -0.5, 0.0]),
920            "M": np.array([0.0, 0.0, 0.5]),
921            "N": np.array([-0.5, -0.5, 0.5]),
922            "R": np.array([0.0, -0.5, 0.5]),
923            "X": np.array([0.0, -0.5, 0.0]),
924            "Y": np.array([0.5, 0.0, 0.0]),
925            "Z": np.array([-0.5, 0.0, 0.5]),
926        }
927        path = [
928            ["X", "\\Gamma", "Y"],
929            ["L", "\\Gamma", "Z"],
930            ["N", "\\Gamma", "M"],
931            ["R", "\\Gamma"],
932        ]
933        return {"kpoints": kpoints, "path": path}
934
935
936class KPathSeek(KPathBase):
937    """
938    This class looks for path along high symmetry lines in
939    the Brillouin Zone.
940    It is based on Hinuma, Y., Pizzi, G., Kumagai, Y., Oba, F.,
941    & Tanaka, I. (2017). Band structure diagram paths based on
942    crystallography. Computational Materials Science, 128, 140–184.
943    https://doi.org/10.1016/j.commatsci.2016.10.015
944    It should be used with primitive structures that
945    comply with the definition from the paper.
946    The symmetry is determined by spglib through the
947    SpacegroupAnalyzer class. KPoints from get_kpoints() method
948    are returned in the reciprocal cell basis defined in the paper.
949    """
950
951    @requires(
952        get_path is not None,
953        "SeeK-path is required to use the convention by Hinuma et al.",
954    )
955    def __init__(self, structure, symprec=0.01, angle_tolerance=5, atol=1e-5, system_is_tri=True):
956        """
957        Args:
958            structure (Structure): Structure object
959            symprec (float): Tolerance for symmetry finding
960            angle_tolerance (float): Angle tolerance for symmetry finding.
961            atol (float): Absolute tolerance used to determine edge cases
962                for settings of structures.
963            system_is_tri (boolean): Indicates if the system is time-reversal
964                invariant.
965        """
966
967        super().__init__(structure, symprec=symprec, angle_tolerance=angle_tolerance, atol=atol)
968
969        positions = structure.frac_coords
970
971        sp = structure.site_properties
972        species = [site.species for site in structure]
973        site_data = species
974
975        if not system_is_tri:
976            warn("Non-zero 'magmom' data will be used to define unique atoms in the cell.")
977            site_data = zip(species, [tuple(vec) for vec in sp["magmom"]])
978
979        unique_species = []
980        numbers = []
981
982        for species, g in itertools.groupby(site_data):
983            if species in unique_species:
984                ind = unique_species.index(species)
985                numbers.extend([ind + 1] * len(tuple(g)))
986            else:
987                unique_species.append(species)
988                numbers.extend([len(unique_species)] * len(tuple(g)))
989
990        cell = (self._latt.matrix, positions, numbers)
991
992        lattice, scale_pos, atom_num = spglib.standardize_cell(
993            cell, to_primitive=False, no_idealize=True, symprec=symprec
994        )
995
996        spg_struct = (lattice, scale_pos, atom_num)
997        spath_dat = get_path(spg_struct, system_is_tri, "hpkot", atol, symprec, angle_tolerance)
998
999        self._tmat = self._trans_sc_to_Hin(spath_dat["bravais_lattice_extended"])
1000        self._rec_lattice = Lattice(spath_dat["reciprocal_primitive_lattice"])
1001
1002        spath_data_formatted = [[spath_dat["path"][0][0]]]
1003        count = 0
1004        for pnum in range(len(spath_dat["path"]) - 1):
1005            if spath_dat["path"][pnum][1] == spath_dat["path"][pnum + 1][0]:
1006                spath_data_formatted[count].append(spath_dat["path"][pnum][1])
1007            else:
1008                spath_data_formatted[count].append(spath_dat["path"][pnum][1])
1009                spath_data_formatted.append([])
1010                count += 1
1011                spath_data_formatted[count].append(spath_dat["path"][pnum + 1][0])
1012
1013        spath_data_formatted[-1].append(spath_dat["path"][-1][1])
1014
1015        self._kpath = {
1016            "kpoints": spath_dat["point_coords"],
1017            "path": spath_data_formatted,
1018        }
1019
1020    @staticmethod
1021    def _trans_sc_to_Hin(sub_class):
1022
1023        if sub_class in [
1024            "cP1",
1025            "cP2",
1026            "cF1",
1027            "cF2",
1028            "cI1",
1029            "tP1",
1030            "oP1",
1031            "hP1",
1032            "hP2",
1033            "tI1",
1034            "tI2",
1035            "oF1",
1036            "oF3",
1037            "oI1",
1038            "oI3",
1039            "oC1",
1040            "hR1",
1041            "hR2",
1042            "aP1",
1043            "aP2",
1044            "aP3",
1045            "oA1",
1046        ]:
1047            return np.eye(3)
1048        if sub_class == "oF2":
1049            return np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]])
1050        if sub_class == "oI2":
1051            return np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]])
1052        if sub_class == "oI3":
1053            return np.array([[0, 0, 1], [1, 0, 0], [0, 1, 0]])
1054        if sub_class == "oA2":
1055            return np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]])
1056        if sub_class == "oC2":
1057            return np.array([[-1, 0, 0], [0, 1, 0], [0, 0, -1]])
1058        if sub_class in ["mP1", "mC1", "mC2", "mC3"]:
1059            return np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]])
1060        raise RuntimeError("Sub-classification of crystal not found!")
1061
1062
1063class KPathLatimerMunro(KPathBase):
1064    """
1065    This class looks for a path along high symmetry lines in the
1066    Brillouin zone. It is based on the method outlined in:
1067
1068    npj Comput Mater 6, 112 (2020). 10.1038/s41524-020-00383-7
1069
1070    The user should ensure that the lattice of the input structure
1071    is as reduced as possible, i.e. that there is no linear
1072    combination of lattice vectors which can produce a vector of
1073    lesser magnitude than the given set (this is required to
1074    obtain the correct Brillouin zone within the current
1075    implementaiton). This is checked during initialization and a
1076    warning is issued if the condition is not fulfilled.
1077    In the case of magnetic structures, care must also be taken to
1078    provide the magnetic primitive cell (i.e. that which reproduces
1079    the entire crystal, including the correct magnetic ordering,
1080    upon application of lattice translations). There is no way to
1081    programatically check for this, so if the input structure is
1082    incorrect, the class will output the incorrect kpath without
1083    any warning being issued.
1084    """
1085
1086    def __init__(
1087        self,
1088        structure,
1089        has_magmoms=False,
1090        magmom_axis=None,
1091        symprec=0.01,
1092        angle_tolerance=5,
1093        atol=1e-5,
1094    ):
1095        """
1096        Args:
1097            structure (Structure): Structure object
1098            has_magmoms (boolean): Whether the input structure contains
1099                magnetic moments as site properties with the key 'magmom.'
1100                Values may be in the form of 3-component vectors given in
1101                the basis of the input lattice vectors, or as scalars, in
1102                which case the spin axis will default to a_3, the third
1103                real-space lattice vector (this triggers a warning).
1104            magmom_axis (list or numpy array): 3-component vector specifying
1105                direction along which magnetic moments given as scalars
1106                should point. If all magnetic moments are provided as
1107                vectors then this argument is not used.
1108            symprec (float): Tolerance for symmetry finding
1109            angle_tolerance (float): Angle tolerance for symmetry finding.
1110            atol (float): Absolute tolerance used to determine symmetric
1111                equivalence of points and lines on the BZ.
1112        """
1113        super().__init__(structure, symprec=symprec, angle_tolerance=angle_tolerance, atol=atol)
1114
1115        # Check to see if input lattice is reducible. Ref: B Gruber in Acta. Cryst. Vol. A29,
1116        # pp. 433-440 ('The Relationship between Reduced Cells in a General Bravais lattice').
1117        # The correct BZ will still be obtained if the lattice vectors are reducible by any
1118        # linear combination of themselves with coefficients of absolute value less than 2,
1119        # hence a missing factor of 2 as compared to the reference.
1120        reducible = []
1121        for i in range(3):
1122            for j in range(3):
1123                if i != j:
1124                    if (
1125                        np.absolute(np.dot(self._latt.matrix[i], self._latt.matrix[j]))
1126                        > np.dot(self._latt.matrix[i], self._latt.matrix[i])
1127                        and np.absolute(
1128                            np.dot(self._latt.matrix[i], self._latt.matrix[j])
1129                            - np.dot(self._latt.matrix[i], self._latt.matrix[i])
1130                        )
1131                        > atol
1132                    ):
1133                        reducible.append(True)
1134                    else:
1135                        reducible.append(False)
1136        if np.any(reducible):
1137            print("reducible")
1138            warn(
1139                "The lattice of the input structure is not fully reduced!"
1140                "The path can be incorrect. Use at your own risk."
1141            )
1142
1143        if magmom_axis is None:
1144            magmom_axis = np.array([0, 0, 1])
1145            axis_specified = False
1146        else:
1147            axis_specified = True
1148
1149        self._kpath = self._get_ksymm_kpath(has_magmoms, magmom_axis, axis_specified, symprec, angle_tolerance, atol)
1150
1151    @property
1152    def mag_type(self):
1153        """
1154        Returns:
1155        The type of magnetic space group as a string.
1156        Current implementation does not distinguish
1157        between types 3 and 4, so return value is '3/4'.
1158        If has_magmoms is False, returns '0'.
1159        """
1160        return self._mag_type
1161
1162    def _get_ksymm_kpath(self, has_magmoms, magmom_axis, axis_specified, symprec, angle_tolerance, atol):
1163        ID = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
1164        # parity, aka the inversion operation (not calling it
1165        PAR = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, -1]])
1166        # INV to avoid confusion with np.linalg.inv() function)
1167
1168        # 1: Get lattices of real and reciprocal structures, and reciprocal
1169        # point group, and Brillouin zone (BZ)
1170
1171        V = self._latt.matrix.T  # fractional real space to cartesian real space
1172        # fractional reciprocal space to cartesian reciprocal space
1173        W = self._rec_lattice.matrix.T
1174        # fractional real space to fractional reciprocal space
1175        A = np.dot(np.linalg.inv(W), V)
1176
1177        if has_magmoms:
1178            grey_struct = self._structure.copy()
1179            grey_struct.remove_site_property("magmom")
1180            sga = SpacegroupAnalyzer(grey_struct, symprec=symprec, angle_tolerance=angle_tolerance)
1181            grey_ops = sga.get_symmetry_operations()
1182            self._structure = self._convert_all_magmoms_to_vectors(magmom_axis, axis_specified)
1183            mag_ops = self._get_magnetic_symmetry_operations(self._structure, grey_ops, atol)
1184
1185            D = [
1186                SymmOp.from_rotation_and_translation(
1187                    rotation_matrix=op.rotation_matrix,
1188                    translation_vec=op.translation_vector,
1189                )
1190                for op in mag_ops
1191                if op.time_reversal == 1
1192            ]
1193
1194            fD = [
1195                SymmOp.from_rotation_and_translation(
1196                    rotation_matrix=op.rotation_matrix,
1197                    translation_vec=op.translation_vector,
1198                )
1199                for op in mag_ops
1200                if op.time_reversal == -1
1201            ]
1202
1203            if np.array([m == np.array([0, 0, 0]) for m in self._structure.site_properties["magmom"]]).all():
1204                fD = D
1205                D = []
1206
1207            if len(fD) == 0:  # no operations contain time reversal; type 1
1208                self._mag_type = "1"
1209                isomorphic_point_group = [d.rotation_matrix for d in D]
1210                recip_point_group = self._get_reciprocal_point_group(isomorphic_point_group, ID, A)
1211            elif len(D) == 0:  # all operations contain time reversal / all magmoms zero; type 2
1212                self._mag_type = "2"
1213                isomorphic_point_group = [d.rotation_matrix for d in fD]
1214                recip_point_group = self._get_reciprocal_point_group(isomorphic_point_group, PAR, A)
1215            else:  # half and half; type 3 or 4
1216                self._mag_type = "3/4"
1217                f = self._get_coset_factor(D + fD, D)
1218                isomorphic_point_group = [d.rotation_matrix for d in D]
1219                recip_point_group = self._get_reciprocal_point_group(
1220                    isomorphic_point_group, np.dot(PAR, f.rotation_matrix), A
1221                )
1222        else:
1223            self._mag_type = "0"
1224            if "magmom" in self._structure.site_properties:
1225                warn(
1226                    "The parameter has_magmoms is False, but site_properties contains the key magmom."
1227                    "This property will be removed and could result in different symmetry operations."
1228                )
1229                self._structure.remove_site_property("magmom")
1230            sga = SpacegroupAnalyzer(self._structure)
1231            ops = sga.get_symmetry_operations()
1232            isomorphic_point_group = [op.rotation_matrix for op in ops]
1233
1234            recip_point_group = self._get_reciprocal_point_group(isomorphic_point_group, PAR, A)
1235
1236        self._rpg = recip_point_group
1237
1238        # 2: Get all vertices, edge- and face- center points of BZ ("key points")
1239
1240        key_points, bz_as_key_point_inds, face_center_inds = self._get_key_points()
1241
1242        # 3: Find symmetry-equivalent points, which can be mapped to each other by a combination of point group
1243        # operations and integer translations by lattice vectors. The integers will only be -1, 0, or 1, since
1244        # we are restricting to the BZ.
1245
1246        key_points_inds_orbits = self._get_key_point_orbits(key_points=key_points)
1247
1248        # 4: Get all lines on BZ between adjacent key points and between gamma
1249        # and key points ("key lines")
1250
1251        key_lines = self._get_key_lines(key_points=key_points, bz_as_key_point_inds=bz_as_key_point_inds)
1252
1253        # 5: Find symmetry-equivalent key lines, defined as endpoints of first line being equivalent
1254        # to end points of second line, and a random point in between being equivalent to the mapped
1255        # random point.
1256
1257        key_lines_inds_orbits = self._get_key_line_orbits(
1258            key_points=key_points,
1259            key_lines=key_lines,
1260            key_points_inds_orbits=key_points_inds_orbits,
1261        )
1262
1263        # 6 & 7: Get little groups for key points (group of symmetry elements present at that point).
1264        # Get little groups for key lines (group of symmetry elements present at every point
1265        # along the line). This is implemented by testing the symmetry at a point e/pi of the
1266        # way between the two endpoints.
1267
1268        little_groups_points, little_groups_lines = self._get_little_groups(
1269            key_points=key_points,
1270            key_points_inds_orbits=key_points_inds_orbits,
1271            key_lines_inds_orbits=key_lines_inds_orbits,
1272        )
1273
1274        # 8: Choose key lines for k-path. Loose criteria set: choose any points / segments
1275        # with spatial symmetry greater than the general position (AKA more symmetry operations
1276        # than just the identity or identity * TR in the little group).
1277        # This function can be edited to alter high-symmetry criteria for choosing points and lines
1278
1279        point_orbits_in_path, line_orbits_in_path = self._choose_path(
1280            key_points=key_points,
1281            key_points_inds_orbits=key_points_inds_orbits,
1282            key_lines_inds_orbits=key_lines_inds_orbits,
1283            little_groups_points=little_groups_points,
1284            little_groups_lines=little_groups_lines,
1285        )
1286
1287        # 10: Consolidate selected segments into a single irreducible section of the Brilouin zone (as determined
1288        # by the reciprocal point and lattice symmetries). This is accomplished by identifying the boundary
1289        # planes of the IRBZ. Also, get labels for points according to distance away from axes.
1290
1291        IRBZ_points_inds = self._get_IRBZ(recip_point_group, W, key_points, face_center_inds, atol)
1292        lines_in_path_inds = []
1293        for ind in line_orbits_in_path:
1294            for tup in key_lines_inds_orbits[ind]:
1295                if tup[0] in IRBZ_points_inds and tup[1] in IRBZ_points_inds:
1296                    lines_in_path_inds.append(tup)
1297                    break
1298        G = nx.Graph(lines_in_path_inds)
1299        lines_in_path_inds = list(nx.edge_dfs(G))
1300        points_in_path_inds = [ind for tup in lines_in_path_inds for ind in tup]
1301        points_in_path_inds_unique = list(set(points_in_path_inds))
1302
1303        orbit_cosines = []
1304        for i, orbit in enumerate(key_points_inds_orbits[:-1]):
1305
1306            orbit_cosines.append(
1307                sorted(
1308                    sorted(
1309                        [
1310                            (
1311                                j,
1312                                np.round(
1313                                    np.dot(key_points[k], self.LabelPoints(j))
1314                                    / (np.linalg.norm(key_points[k]) * np.linalg.norm(self.LabelPoints(j))),
1315                                    decimals=3,
1316                                ),
1317                            )
1318                            for k in orbit
1319                            for j in range(26)
1320                        ],
1321                        key=operator.itemgetter(0),
1322                    ),
1323                    key=operator.itemgetter(1),
1324                    reverse=True,
1325                )
1326            )
1327
1328        orbit_labels = self._get_orbit_labels(orbit_cosines, key_points_inds_orbits, atol)
1329        key_points_labels = ["" for i in range(len(key_points))]
1330        for i, orbit in enumerate(key_points_inds_orbits):
1331            for point_ind in orbit:
1332                key_points_labels[point_ind] = self.LabelSymbol(int(orbit_labels[i]))
1333
1334        kpoints = {}
1335        reverse_kpoints = {}
1336        for point_ind in points_in_path_inds_unique:
1337            point_label = key_points_labels[point_ind]
1338            if point_label not in kpoints.keys():
1339                kpoints[point_label] = key_points[point_ind]
1340                reverse_kpoints[point_ind] = point_label
1341            else:
1342                existing_labels = [key for key in kpoints.keys() if point_label in key]
1343                if "'" not in point_label:
1344                    existing_labels[:] = [label for label in existing_labels if "'" not in label]
1345                if len(existing_labels) == 1:
1346                    max_occurence = 0
1347                else:
1348                    if "'" not in point_label:
1349                        max_occurence = max([int(label[3:-1]) for label in existing_labels[1:]])
1350                    else:
1351                        max_occurence = max([int(label[4:-1]) for label in existing_labels[1:]])
1352                kpoints[point_label + "_{" + str(max_occurence + 1) + "}"] = key_points[point_ind]
1353                reverse_kpoints[point_ind] = point_label + "_{" + str(max_occurence + 1) + "}"
1354
1355        path = []
1356        i = 0
1357        start_of_subpath = True
1358        while i < len(points_in_path_inds):
1359            if start_of_subpath:
1360                path.append([reverse_kpoints[points_in_path_inds[i]]])
1361                i += 1
1362                start_of_subpath = False
1363            elif points_in_path_inds[i] == points_in_path_inds[i + 1]:
1364                path[-1].append(reverse_kpoints[points_in_path_inds[i]])
1365                i += 2
1366            else:
1367                path[-1].append(reverse_kpoints[points_in_path_inds[i]])
1368                i += 1
1369                start_of_subpath = True
1370            if i == len(points_in_path_inds) - 1:
1371                path[-1].append(reverse_kpoints[points_in_path_inds[i]])
1372                i += 1
1373
1374        return {"kpoints": kpoints, "path": path}
1375
1376    def _choose_path(
1377        self,
1378        key_points,
1379        key_points_inds_orbits,
1380        key_lines_inds_orbits,
1381        little_groups_points,
1382        little_groups_lines,
1383    ):
1384        #
1385        # This function can be edited to alter high-symmetry criteria for choosing points and lines
1386        #
1387
1388        ID = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
1389        PAR = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, -1]])
1390
1391        gamma_ind = len(key_points) - 1
1392
1393        line_orbits_in_path = []
1394        point_orbits_in_path = []
1395        for (i, little_group) in enumerate(little_groups_lines):
1396            add_rep = False
1397            nC2 = 0
1398            nC3 = 0
1399            nsig = 0
1400            for j, opind in enumerate(little_group):
1401                op = self._rpg[opind]
1402                if not (op == ID).all():
1403                    if (np.dot(op, op) == ID).all():
1404                        if np.linalg.det(op) == 1:
1405                            nC2 += 1
1406                            break
1407                        if not (op == PAR).all():
1408                            nsig += 1
1409                            break
1410                    elif (np.dot(op, np.dot(op, op)) == ID).all():
1411                        nC3 += 1
1412                        break
1413            if nC2 > 0 or nC3 > 0 or nsig > 0:
1414                add_rep = True
1415
1416            if add_rep:
1417                line_orbits_in_path.append(i)
1418                l = key_lines_inds_orbits[i][0]
1419                ind0 = l[0]
1420                ind1 = l[1]
1421                found0 = False
1422                found1 = False
1423                for (j, orbit) in enumerate(key_points_inds_orbits):
1424                    if ind0 in orbit:
1425                        point_orbits_in_path.append(j)
1426                        found0 = True
1427                    if ind1 in orbit:
1428                        point_orbits_in_path.append(j)
1429                        found1 = True
1430                    if found0 and found1:
1431                        break
1432        point_orbits_in_path = list(set(point_orbits_in_path))
1433
1434        # Choose remaining unconnected key points for k-path. The ones that remain are
1435        # those with inversion symmetry. Connect them to gamma.
1436
1437        unconnected = []
1438
1439        for i in range(len(key_points_inds_orbits)):
1440            if i not in point_orbits_in_path:
1441                unconnected.append(i)
1442
1443        for ind in unconnected:
1444            connect = False
1445            for op_ind in little_groups_points[ind]:
1446                op = self._rpg[op_ind]
1447                if (op == ID).all():
1448                    pass
1449                elif (op == PAR).all():
1450                    connect = True
1451                    break
1452                elif np.linalg.det(op) == 1:
1453                    if (np.dot(op, np.dot(op, op)) == ID).all():
1454                        pass
1455                    else:
1456                        connect = True
1457                        break
1458                else:
1459                    pass
1460            if connect:
1461                l = (key_points_inds_orbits[ind][0], gamma_ind)
1462                for (j, orbit) in enumerate(key_lines_inds_orbits):
1463                    if l in orbit:
1464                        line_orbits_in_path.append(j)
1465                        break
1466                if gamma_ind not in point_orbits_in_path:
1467                    point_orbits_in_path.append(gamma_ind)
1468                point_orbits_in_path.append(ind)
1469
1470        return point_orbits_in_path, line_orbits_in_path
1471
1472    def _get_key_points(self):
1473        decimals = ceil(-1 * np.log10(self._atol)) - 1
1474        bz = self._rec_lattice.get_wigner_seitz_cell()
1475
1476        key_points = []
1477        face_center_inds = []
1478        bz_as_key_point_inds = []
1479
1480        # pymatgen gives BZ in cartesian coordinates; convert to fractional in
1481        # the primitive basis for reciprocal space
1482        for (i, facet) in enumerate(bz):
1483            for (j, vert) in enumerate(facet):
1484                vert = self._rec_lattice.get_fractional_coords(vert)
1485                bz[i][j] = vert
1486        pop = []
1487        for i, facet in enumerate(bz):
1488            rounded_facet = np.around(facet, decimals=decimals)
1489            u, indices = np.unique(rounded_facet, axis=0, return_index=True)
1490            if len(u) in [1, 2]:
1491                pop.append(i)
1492            else:
1493                bz[i] = [facet[j] for j in np.sort(indices)]
1494        bz = [bz[i] for i in range(len(bz)) if i not in pop]
1495
1496        # use vertex points to calculate edge- and face- centers
1497        for (i, facet) in enumerate(bz):
1498            bz_as_key_point_inds.append([])
1499            for (j, vert) in enumerate(facet):
1500                edge_center = (vert + facet[j + 1]) / 2.0 if j != len(facet) - 1 else (vert + facet[0]) / 2.0
1501                duplicatevert = False
1502                duplicateedge = False
1503                for (k, point) in enumerate(key_points):
1504                    if np.allclose(vert, point, atol=self._atol):
1505                        bz_as_key_point_inds[i].append(k)
1506                        duplicatevert = True
1507                        break
1508                for (k, point) in enumerate(key_points):
1509                    if np.allclose(edge_center, point, atol=self._atol):
1510                        bz_as_key_point_inds[i].append(k)
1511                        duplicateedge = True
1512                        break
1513                if not duplicatevert:
1514                    key_points.append(vert)
1515                    bz_as_key_point_inds[i].append(len(key_points) - 1)
1516                if not duplicateedge:
1517                    key_points.append(edge_center)
1518                    bz_as_key_point_inds[i].append(len(key_points) - 1)
1519            if len(facet) == 4:  # parallelogram facet
1520                face_center = (facet[0] + facet[1] + facet[2] + facet[3]) / 4.0
1521                key_points.append(face_center)
1522                face_center_inds.append(len(key_points) - 1)
1523                bz_as_key_point_inds[i].append(len(key_points) - 1)
1524            else:  # hexagonal facet
1525                face_center = (facet[0] + facet[1] + facet[2] + facet[3] + facet[4] + facet[5]) / 6.0
1526                key_points.append(face_center)
1527                face_center_inds.append(len(key_points) - 1)
1528                bz_as_key_point_inds[i].append(len(key_points) - 1)
1529
1530        # add gamma point
1531        key_points.append(np.array([0, 0, 0]))
1532        return key_points, bz_as_key_point_inds, face_center_inds
1533
1534    def _get_key_point_orbits(self, key_points):
1535        key_points_copy = dict(zip(range(len(key_points) - 1), key_points[0 : len(key_points) - 1]))
1536        # gamma not equivalent to any on BZ and is last point added to
1537        # key_points
1538        key_points_inds_orbits = []
1539
1540        i = 0
1541        while len(key_points_copy) > 0:
1542            key_points_inds_orbits.append([])
1543            k0ind = list(key_points_copy.keys())[0]
1544            k0 = key_points_copy[k0ind]
1545            key_points_inds_orbits[i].append(k0ind)
1546            key_points_copy.pop(k0ind)
1547
1548            for op in self._rpg:
1549                to_pop = []
1550                k1 = np.dot(op, k0)
1551                for ind_key in key_points_copy:
1552                    diff = k1 - key_points_copy[ind_key]
1553                    if self._all_ints(diff, atol=self._atol):
1554                        key_points_inds_orbits[i].append(ind_key)
1555                        to_pop.append(ind_key)
1556
1557                for key in to_pop:
1558                    key_points_copy.pop(key)
1559            i += 1
1560
1561        key_points_inds_orbits.append([len(key_points) - 1])
1562
1563        return key_points_inds_orbits
1564
1565    @staticmethod
1566    def _get_key_lines(key_points, bz_as_key_point_inds):
1567        key_lines = []
1568        gamma_ind = len(key_points) - 1
1569
1570        for (i, facet_as_key_point_inds) in enumerate(bz_as_key_point_inds):
1571            facet_as_key_point_inds_bndy = facet_as_key_point_inds[: len(facet_as_key_point_inds) - 1]
1572            # not the face center point (don't need to check it since it's not
1573            # shared with other facets)
1574            face_center_ind = facet_as_key_point_inds[-1]
1575            for (j, ind) in enumerate(facet_as_key_point_inds_bndy):
1576                if (
1577                    min(ind, facet_as_key_point_inds_bndy[j - 1]),
1578                    max(ind, facet_as_key_point_inds_bndy[j - 1]),
1579                ) not in key_lines:
1580                    key_lines.append(
1581                        (
1582                            min(ind, facet_as_key_point_inds_bndy[j - 1]),
1583                            max(ind, facet_as_key_point_inds_bndy[j - 1]),
1584                        )
1585                    )
1586                k = j + 1 if j != len(facet_as_key_point_inds_bndy) - 1 else 0
1587                if (
1588                    min(ind, facet_as_key_point_inds_bndy[k]),
1589                    max(ind, facet_as_key_point_inds_bndy[k]),
1590                ) not in key_lines:
1591                    key_lines.append(
1592                        (
1593                            min(ind, facet_as_key_point_inds_bndy[k]),
1594                            max(ind, facet_as_key_point_inds_bndy[k]),
1595                        )
1596                    )
1597                if (ind, gamma_ind) not in key_lines:
1598                    key_lines.append((ind, gamma_ind))
1599                key_lines.append((min(ind, face_center_ind), max(ind, face_center_ind)))
1600            key_lines.append((face_center_ind, gamma_ind))
1601
1602        return key_lines
1603
1604    def _get_key_line_orbits(self, key_points, key_lines, key_points_inds_orbits):
1605        key_lines_copy = dict(zip(range(len(key_lines)), key_lines))
1606        key_lines_inds_orbits = []
1607
1608        i = 0
1609        while len(key_lines_copy) > 0:
1610            key_lines_inds_orbits.append([])
1611            l0ind = list(key_lines_copy.keys())[0]
1612            l0 = key_lines_copy[l0ind]
1613            key_lines_inds_orbits[i].append(l0)
1614            key_lines_copy.pop(l0ind)
1615            to_pop = []
1616            p00 = key_points[l0[0]]
1617            p01 = key_points[l0[1]]
1618            pmid0 = p00 + e / pi * (p01 - p00)
1619            for ind_key in key_lines_copy:
1620
1621                l1 = key_lines_copy[ind_key]
1622                p10 = key_points[l1[0]]
1623                p11 = key_points[l1[1]]
1624                equivptspar = False
1625                equivptsperp = False
1626                equivline = False
1627
1628                if (
1629                    np.array([l0[0] in orbit and l1[0] in orbit for orbit in key_points_inds_orbits]).any()
1630                    and np.array([l0[1] in orbit and l1[1] in orbit for orbit in key_points_inds_orbits]).any()
1631                ):
1632                    equivptspar = True
1633                elif (
1634                    np.array([l0[1] in orbit and l1[0] in orbit for orbit in key_points_inds_orbits]).any()
1635                    and np.array([l0[0] in orbit and l1[1] in orbit for orbit in key_points_inds_orbits]).any()
1636                ):
1637                    equivptsperp = True
1638
1639                if equivptspar:
1640                    pmid1 = p10 + e / pi * (p11 - p10)
1641                    for op in self._rpg:
1642                        if not equivline:
1643                            p00pr = np.dot(op, p00)
1644                            diff0 = p10 - p00pr
1645                            if self._all_ints(diff0, atol=self._atol):
1646                                pmid0pr = np.dot(op, pmid0) + diff0
1647                                p01pr = np.dot(op, p01) + diff0
1648                                if np.allclose(p11, p01pr, atol=self._atol) and np.allclose(
1649                                    pmid1, pmid0pr, atol=self._atol
1650                                ):
1651                                    equivline = True
1652
1653                elif equivptsperp:
1654                    pmid1 = p11 + e / pi * (p10 - p11)
1655                    for op in self._rpg:
1656                        if not equivline:
1657                            p00pr = np.dot(op, p00)
1658                            diff0 = p11 - p00pr
1659                            if self._all_ints(diff0, atol=self._atol):
1660                                pmid0pr = np.dot(op, pmid0) + diff0
1661                                p01pr = np.dot(op, p01) + diff0
1662                                if np.allclose(p10, p01pr, atol=self._atol) and np.allclose(
1663                                    pmid1, pmid0pr, atol=self._atol
1664                                ):
1665                                    equivline = True
1666
1667                if equivline:
1668                    key_lines_inds_orbits[i].append(l1)
1669                    to_pop.append(ind_key)
1670
1671            for key in to_pop:
1672                key_lines_copy.pop(key)
1673            i += 1
1674
1675        return key_lines_inds_orbits
1676
1677    def _get_little_groups(self, key_points, key_points_inds_orbits, key_lines_inds_orbits):
1678
1679        little_groups_points = []  # elements are lists of indicies of recip_point_group. the
1680        # list little_groups_points[i] is the little group for the
1681        # orbit key_points_inds_orbits[i]
1682        for (i, orbit) in enumerate(key_points_inds_orbits):
1683            k0 = key_points[orbit[0]]
1684            little_groups_points.append([])
1685            for (j, op) in enumerate(self._rpg):
1686                gamma_to = np.dot(op, -1 * k0) + k0
1687                check_gamma = True
1688                if not self._all_ints(gamma_to, atol=self._atol):
1689                    check_gamma = False
1690                if check_gamma:
1691                    little_groups_points[i].append(j)
1692
1693        # elements are lists of indicies of recip_point_group. the list
1694        # little_groups_lines[i] is
1695        little_groups_lines = []
1696        # the little group for the orbit key_points_inds_lines[i]
1697
1698        for (i, orbit) in enumerate(key_lines_inds_orbits):
1699            l0 = orbit[0]
1700            v = key_points[l0[1]] - key_points[l0[0]]
1701            k0 = key_points[l0[0]] + np.e / pi * v
1702            little_groups_lines.append([])
1703            for (j, op) in enumerate(self._rpg):
1704                gamma_to = np.dot(op, -1 * k0) + k0
1705                check_gamma = True
1706                if not self._all_ints(gamma_to, atol=self._atol):
1707                    check_gamma = False
1708                if check_gamma:
1709                    little_groups_lines[i].append(j)
1710
1711        return little_groups_points, little_groups_lines
1712
1713    def _convert_all_magmoms_to_vectors(self, magmom_axis, axis_specified):
1714        struct = self._structure.copy()
1715        magmom_axis = np.array(magmom_axis)
1716        if "magmom" not in struct.site_properties:
1717            warn(
1718                "The 'magmom' property is not set in the structure's site properties."
1719                "All magnetic moments are being set to zero."
1720            )
1721            struct.add_site_property("magmom", [np.array([0, 0, 0]) for i in range(len(struct.sites))])
1722
1723            return struct
1724
1725        old_magmoms = struct.site_properties["magmom"]
1726        new_magmoms = []
1727        found_scalar = False
1728
1729        for magmom in old_magmoms:
1730            if isinstance(magmom, np.ndarray):
1731                new_magmoms.append(magmom)
1732            elif isinstance(magmom, list):
1733                new_magmoms.append(np.array(magmom))
1734            else:
1735                found_scalar = True
1736                new_magmoms.append(magmom * magmom_axis)
1737
1738        if found_scalar and not axis_specified:
1739            warn("At least one magmom had a scalar value and magmom_axis was not specified." "Defaulted to z+ spinor.")
1740
1741        struct.remove_site_property("magmom")
1742        struct.add_site_property("magmom", new_magmoms)
1743        return struct
1744
1745    def _get_magnetic_symmetry_operations(self, struct, grey_ops, atol):
1746        mag_ops = []
1747        magmoms = struct.site_properties["magmom"]
1748        nonzero_magmom_inds = [i for i in range(len(struct.sites)) if not (magmoms[i] == np.array([0, 0, 0])).all()]
1749        init_magmoms = [site.properties["magmom"] for (i, site) in enumerate(struct.sites) if i in nonzero_magmom_inds]
1750        sites = [site for (i, site) in enumerate(struct.sites) if i in nonzero_magmom_inds]
1751        init_site_coords = [site.frac_coords for site in sites]
1752        for op in grey_ops:
1753            r = op.rotation_matrix
1754            t = op.translation_vector
1755            xformed_magmoms = [self._apply_op_to_magmom(r, magmom) for magmom in init_magmoms]
1756            xformed_site_coords = [np.dot(r, site.frac_coords) + t for site in sites]
1757            permutation = ["a" for i in range(len(sites))]
1758            not_found = list(range(len(sites)))
1759            for i in range(len(sites)):
1760                xformed = xformed_site_coords[i]
1761                for k, j in enumerate(not_found):
1762                    init = init_site_coords[j]
1763                    diff = xformed - init
1764                    if self._all_ints(diff, atol=atol):
1765                        permutation[i] = j
1766                        not_found.pop(k)
1767                        break
1768
1769            same = np.zeros(len(sites))
1770            flipped = np.zeros(len(sites))
1771            for i, magmom in enumerate(xformed_magmoms):
1772                if (magmom == init_magmoms[permutation[i]]).all():
1773                    same[i] = 1
1774                elif (magmom == -1 * init_magmoms[permutation[i]]).all():
1775                    flipped[i] = 1
1776
1777            if same.all():  # add symm op without tr
1778                mag_ops.append(
1779                    MagSymmOp.from_rotation_and_translation_and_time_reversal(
1780                        rotation_matrix=op.rotation_matrix,
1781                        translation_vec=op.translation_vector,
1782                        time_reversal=1,
1783                    )
1784                )
1785            if flipped.all():  # add symm op with tr
1786                mag_ops.append(
1787                    MagSymmOp.from_rotation_and_translation_and_time_reversal(
1788                        rotation_matrix=op.rotation_matrix,
1789                        translation_vec=op.translation_vector,
1790                        time_reversal=-1,
1791                    )
1792                )
1793
1794        return mag_ops
1795
1796    @staticmethod
1797    def _get_reciprocal_point_group(ops, R, A):
1798        Ainv = np.linalg.inv(A)
1799        # convert to reciprocal primitive basis
1800        recip_point_group = [np.around(np.dot(A, np.dot(R, Ainv)), decimals=2)]
1801        for op in ops:
1802            op = np.around(np.dot(A, np.dot(op, Ainv)), decimals=2)
1803            new = True
1804            new_coset = True
1805            for thing in recip_point_group:
1806                if (thing == op).all():
1807                    new = False
1808                if (thing == np.dot(R, op)).all():
1809                    new_coset = False
1810
1811            if new:
1812                recip_point_group.append(op)
1813            if new_coset:
1814                recip_point_group.append(np.dot(R, op))
1815
1816        return recip_point_group
1817
1818    @staticmethod
1819    def _closewrapped(pos1, pos2, tolerance):
1820        pos1 = pos1 % 1.0
1821        pos2 = pos2 % 1.0
1822
1823        if len(pos1) != len(pos2):
1824            return False
1825        for i, v in enumerate(pos1):
1826            if abs(pos1[i] - pos2[i]) > tolerance[i] and abs(pos1[i] - pos2[i]) < 1.0 - tolerance[i]:
1827                return False
1828        return True
1829
1830    def _get_coset_factor(self, G, H):
1831        # finds g for left coset decomposition G = H + gH (H must be subgroup of G with index two.)
1832        # in this implementation, G and H are lists of objects of type
1833        # SymmOp
1834        gH = []
1835        for i, op1 in enumerate(G):
1836            in_H = False
1837            for op2 in H:
1838                if np.allclose(op1.rotation_matrix, op2.rotation_matrix, atol=self._atol) and self._closewrapped(
1839                    op1.translation_vector,
1840                    op2.translation_vector,
1841                    np.ones(3) * self._atol,
1842                ):
1843                    in_H = True
1844                    break
1845            if not in_H:
1846                gH.append(op1)
1847
1848        for op in gH:
1849            opH = [op.__mul__(h) for h in H]
1850            is_coset_factor = True
1851            for op1 in opH:
1852                for op2 in H:
1853                    if np.allclose(op1.rotation_matrix, op2.rotation_matrix, atol=self._atol) and self._closewrapped(
1854                        op1.translation_vector,
1855                        op2.translation_vector,
1856                        np.ones(3) * self._atol,
1857                    ):
1858
1859                        is_coset_factor = False
1860                        break
1861                if not is_coset_factor:
1862                    break
1863            if is_coset_factor:
1864                return op
1865
1866        return "No coset factor found."
1867
1868    @staticmethod
1869    def _apply_op_to_magmom(r, magmom):
1870        if np.linalg.det(r) == 1:
1871            return np.dot(r, magmom)
1872        return -1 * np.dot(r, magmom)
1873
1874    @staticmethod
1875    def _all_ints(arr, atol):
1876        rounded_arr = np.around(arr, decimals=0)
1877        return np.allclose(rounded_arr, arr, atol=atol)
1878
1879    def _get_IRBZ(self, recip_point_group, W, key_points, face_center_inds, atol):
1880        rpgdict = self._get_reciprocal_point_group_dict(recip_point_group, atol)
1881
1882        g = np.dot(W.T, W)  # just using change of basis matrix rather than
1883        # Lattice.get_cartesian_coordinates for conciseness
1884        ginv = np.linalg.inv(g)
1885        D = np.linalg.det(W)
1886
1887        primary_orientation = None
1888        secondary_orientation = None
1889        tertiary_orientation = None
1890
1891        planar_boundaries = []
1892        IRBZ_points = list(enumerate(key_points))
1893
1894        for sigma in rpgdict["reflections"]:
1895            norm = sigma["normal"]
1896            if primary_orientation is None:
1897                primary_orientation = norm
1898                planar_boundaries.append(norm)
1899            elif np.isclose(np.dot(primary_orientation, np.dot(g, norm)), 0, atol=atol):
1900                if secondary_orientation is None:
1901                    secondary_orientation = norm
1902                    planar_boundaries.append(norm)
1903                elif np.isclose(np.dot(secondary_orientation, np.dot(g, norm)), 0, atol=atol):
1904                    if tertiary_orientation is None:
1905                        tertiary_orientation = norm
1906                        planar_boundaries.append(norm)
1907                    elif np.allclose(norm, -1 * tertiary_orientation, atol=atol):
1908                        pass
1909                elif np.dot(secondary_orientation, np.dot(g, norm)) < 0:
1910                    planar_boundaries.append(-1 * norm)
1911                else:
1912                    planar_boundaries.append(norm)
1913            elif np.dot(primary_orientation, np.dot(g, norm)) < 0:
1914                planar_boundaries.append(-1 * norm)
1915            else:
1916                planar_boundaries.append(norm)
1917
1918        IRBZ_points = self._reduce_IRBZ(IRBZ_points, planar_boundaries, g, atol)
1919
1920        used_axes = []
1921
1922        # six-fold rotoinversion always comes with horizontal mirror so don't
1923        # need to check
1924        for rotn in rpgdict["rotations"]["six-fold"]:
1925            ax = rotn["axis"]
1926            op = rotn["op"]
1927            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
1928                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
1929                    face_center_found = False
1930                    for point in IRBZ_points:
1931                        if point[0] in face_center_inds:
1932                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
1933                            if not np.allclose(cross, 0, atol=atol):
1934                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
1935                                face_center_found = True
1936                                used_axes.append(ax)
1937                                break
1938                    if not face_center_found:
1939                        print("face center not found")
1940                        for point in IRBZ_points:
1941                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
1942                            if not np.allclose(cross, 0, atol=atol):
1943                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
1944                                used_axes.append(ax)
1945                                break
1946                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)
1947
1948        for rotn in rpgdict["rotations"]["rotoinv-four-fold"]:
1949            ax = rotn["axis"]
1950            op = rotn["op"]
1951            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
1952                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
1953                    face_center_found = False
1954                    for point in IRBZ_points:
1955                        if point[0] in face_center_inds:
1956                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
1957                            if not np.allclose(cross, 0, atol=atol):
1958                                rot_boundaries = [cross, np.dot(op, cross)]
1959                                face_center_found = True
1960                                used_axes.append(ax)
1961                                break
1962                    if not face_center_found:
1963                        print("face center not found")
1964                        for point in IRBZ_points:
1965                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
1966                            if not np.allclose(cross, 0, atol=atol):
1967                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
1968                                used_axes.append(ax)
1969                                break
1970                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)
1971
1972        for rotn in rpgdict["rotations"]["four-fold"]:
1973            ax = rotn["axis"]
1974            op = rotn["op"]
1975            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
1976                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
1977                    face_center_found = False
1978                    for point in IRBZ_points:
1979                        if point[0] in face_center_inds:
1980                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
1981                            if not np.allclose(cross, 0, atol=atol):
1982                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
1983                                face_center_found = True
1984                                used_axes.append(ax)
1985                                break
1986                    if not face_center_found:
1987                        print("face center not found")
1988                        for point in IRBZ_points:
1989                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
1990                            if not np.allclose(cross, 0, atol=atol):
1991                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
1992                                used_axes.append(ax)
1993                                break
1994                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)
1995
1996        for rotn in rpgdict["rotations"]["rotoinv-three-fold"]:
1997            ax = rotn["axis"]
1998            op = rotn["op"]
1999            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
2000                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
2001                    face_center_found = False
2002                    for point in IRBZ_points:
2003                        if point[0] in face_center_inds:
2004                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
2005                            if not np.allclose(cross, 0, atol=atol):
2006                                rot_boundaries = [
2007                                    cross,
2008                                    -1 * np.dot(sqrtm(-1 * op), cross),
2009                                ]
2010                                face_center_found = True
2011                                used_axes.append(ax)
2012                                break
2013                    if not face_center_found:
2014                        print("face center not found")
2015                        for point in IRBZ_points:
2016                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
2017                            if not np.allclose(cross, 0, atol=atol):
2018                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
2019                                used_axes.append(ax)
2020                                break
2021                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)
2022
2023        for rotn in rpgdict["rotations"]["three-fold"]:
2024            ax = rotn["axis"]
2025            op = rotn["op"]
2026            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
2027                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
2028                    face_center_found = False
2029                    for point in IRBZ_points:
2030                        if point[0] in face_center_inds:
2031                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
2032                            if not np.allclose(cross, 0, atol=atol):
2033                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
2034                                face_center_found = True
2035                                used_axes.append(ax)
2036                                break
2037                    if not face_center_found:
2038                        print("face center not found")
2039                        for point in IRBZ_points:
2040                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
2041                            if not np.allclose(cross, 0, atol=atol):
2042                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
2043                                used_axes.append(ax)
2044                                break
2045                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)
2046
2047        for rotn in rpgdict["rotations"]["two-fold"]:
2048            ax = rotn["axis"]
2049            op = rotn["op"]
2050            if not np.any([np.allclose(ax, usedax, atol) for usedax in used_axes]):
2051                if self._op_maps_IRBZ_to_self(op, IRBZ_points, atol):
2052                    face_center_found = False
2053                    for point in IRBZ_points:
2054                        if point[0] in face_center_inds:
2055                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
2056                            if not np.allclose(cross, 0, atol=atol):
2057                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
2058                                face_center_found = True
2059                                used_axes.append(ax)
2060                                break
2061                    if not face_center_found:
2062                        print("face center not found")
2063                        for point in IRBZ_points:
2064                            cross = D * np.dot(ginv, np.cross(ax, point[1]))
2065                            if not np.allclose(cross, 0, atol=atol):
2066                                rot_boundaries = [cross, -1 * np.dot(op, cross)]
2067                                used_axes.append(ax)
2068                                break
2069                    IRBZ_points = self._reduce_IRBZ(IRBZ_points, rot_boundaries, g, atol)
2070
2071        return [point[0] for point in IRBZ_points]
2072
2073    @staticmethod
2074    def _get_reciprocal_point_group_dict(recip_point_group, atol):
2075        PAR = np.array([[-1, 0, 0], [0, -1, 0], [0, 0, -1]])
2076
2077        d = {
2078            "reflections": [],
2079            "rotations": {
2080                "two-fold": [],
2081                "three-fold": [],
2082                "four-fold": [],
2083                "six-fold": [],
2084                "rotoinv-three-fold": [],
2085                "rotoinv-four-fold": [],
2086                "rotoinv-six-fold": [],
2087            },
2088            "inversion": [],
2089        }
2090
2091        for i, op in enumerate(recip_point_group):
2092            evals, evects = np.linalg.eig(op)
2093
2094            tr = np.trace(op)
2095            det = np.linalg.det(op)
2096
2097            # Proper rotations
2098            if np.isclose(det, 1, atol=atol):
2099                if np.isclose(tr, 3, atol=atol):
2100                    continue
2101                if np.isclose(tr, -1, atol=atol):  # two-fold rotation
2102                    for j in range(3):
2103                        if np.isclose(evals[j], 1, atol=atol):
2104                            ax = evects[:, j]
2105                    d["rotations"]["two-fold"].append({"ind": i, "axis": ax, "op": op})
2106                elif np.isclose(tr, 0, atol=atol):  # three-fold rotation
2107                    for j in range(3):
2108                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
2109                            ax = evects[:, j]
2110                    d["rotations"]["three-fold"].append({"ind": i, "axis": ax, "op": op})
2111                # four-fold rotation
2112                elif np.isclose(tr, 1, atol=atol):
2113                    for j in range(3):
2114                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
2115                            ax = evects[:, j]
2116                    d["rotations"]["four-fold"].append({"ind": i, "axis": ax, "op": op})
2117                elif np.isclose(tr, 2, atol=atol):  # six-fold rotation
2118                    for j in range(3):
2119                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
2120                            ax = evects[:, j]
2121                    d["rotations"]["six-fold"].append({"ind": i, "axis": ax, "op": op})
2122
2123            # Improper rotations
2124            if np.isclose(det, -1, atol=atol):
2125                if np.isclose(tr, -3, atol=atol):
2126                    d["inversion"].append({"ind": i, "op": PAR})
2127                elif np.isclose(tr, 1, atol=atol):  # two-fold rotation
2128                    for j in range(3):
2129                        if np.isclose(evals[j], -1, atol=atol):
2130                            norm = evects[:, j]
2131                    d["reflections"].append({"ind": i, "normal": norm, "op": op})
2132                elif np.isclose(tr, 0, atol=atol):  # three-fold rotoinversion
2133                    for j in range(3):
2134                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
2135                            ax = evects[:, j]
2136                    d["rotations"]["rotoinv-three-fold"].append({"ind": i, "axis": ax, "op": op})
2137                # four-fold rotoinversion
2138                elif np.isclose(tr, -1, atol=atol):
2139                    for j in range(3):
2140                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
2141                            ax = evects[:, j]
2142                    d["rotations"]["rotoinv-four-fold"].append({"ind": i, "axis": ax, "op": op})
2143                # six-fold rotoinversion
2144                elif np.isclose(tr, -2, atol=atol):
2145                    for j in range(3):
2146                        if np.isreal(evals[j]) and np.isclose(np.absolute(evals[j]), 1, atol=atol):
2147                            ax = evects[:, j]
2148                    d["rotations"]["rotoinv-six-fold"].append({"ind": i, "axis": ax, "op": op})
2149
2150        return d
2151
2152    @staticmethod
2153    def _op_maps_IRBZ_to_self(op, IRBZ_points, atol):
2154        point_coords = [point[1] for point in IRBZ_points]
2155        for point in point_coords:
2156            point_prime = np.dot(op, point)
2157            mapped_back = False
2158            for checkpoint in point_coords:
2159                if np.allclose(point_prime, checkpoint, atol):
2160                    mapped_back = True
2161                    break
2162            if not mapped_back:
2163                return False
2164
2165        return True
2166
2167    @staticmethod
2168    def _reduce_IRBZ(IRBZ_points, boundaries, g, atol):
2169        in_reduced_section = []
2170        for point in IRBZ_points:
2171            in_reduced_section.append(
2172                np.all(
2173                    [
2174                        (
2175                            np.dot(point[1], np.dot(g, boundary)) >= 0
2176                            or np.isclose(np.dot(point[1], np.dot(g, boundary)), 0, atol=atol)
2177                        )
2178                        for boundary in boundaries
2179                    ]
2180                )
2181            )
2182
2183        return [IRBZ_points[i] for i in range(len(IRBZ_points)) if in_reduced_section[i]]
2184
2185    def _get_orbit_labels(self, orbit_cosines_orig, key_points_inds_orbits, atol):
2186
2187        orbit_cosines_copy = orbit_cosines_orig.copy()
2188        orbit_labels_unsorted = [(len(key_points_inds_orbits) - 1, 26)]
2189        orbit_inds_remaining = range(len(key_points_inds_orbits) - 1)
2190        pop_orbits = []
2191        pop_labels = []
2192
2193        for i, orb_cos in enumerate(orbit_cosines_copy):
2194            if np.isclose(orb_cos[0][1], 1.0, atol=atol):
2195                # (point orbit index, label index)
2196                orbit_labels_unsorted.append((i, orb_cos[0][0]))
2197                pop_orbits.append(i)
2198                pop_labels.append(orb_cos[0][0])
2199
2200        orbit_cosines_copy = self._reduce_cosines_array(orbit_cosines_copy, pop_orbits, pop_labels)
2201        orbit_inds_remaining = [i for i in orbit_inds_remaining if i not in pop_orbits]
2202
2203        # orbit_labels_unsorted already contains gamma orbit
2204        while len(orbit_labels_unsorted) < len(orbit_cosines_orig) + 1:
2205            pop_orbits = []
2206            pop_labels = []
2207            max_cosine_value = max([orb_cos[0][1] for orb_cos in orbit_cosines_copy])
2208            max_cosine_value_inds = [
2209                j for j in range(len(orbit_cosines_copy)) if orbit_cosines_copy[j][0][1] == max_cosine_value
2210            ]
2211            max_cosine_label_inds = self._get_max_cosine_labels(
2212                [orbit_cosines_copy[j] for j in max_cosine_value_inds],
2213                key_points_inds_orbits,
2214                atol,
2215            )
2216
2217            for j, label_ind in enumerate(max_cosine_label_inds):
2218                orbit_labels_unsorted.append((orbit_inds_remaining[max_cosine_value_inds[j]], label_ind))
2219                pop_orbits.append(max_cosine_value_inds[j])
2220                pop_labels.append(label_ind)
2221            orbit_cosines_copy = self._reduce_cosines_array(orbit_cosines_copy, pop_orbits, pop_labels)
2222            orbit_inds_remaining = [
2223                orbit_inds_remaining[j] for j in range(len(orbit_inds_remaining)) if j not in pop_orbits
2224            ]
2225
2226        orbit_labels = np.zeros(len(key_points_inds_orbits))
2227        for tup in orbit_labels_unsorted:
2228            orbit_labels[tup[0]] = tup[1]
2229
2230        return orbit_labels
2231
2232    @staticmethod
2233    def _reduce_cosines_array(orbit_cosines, pop_orbits, pop_labels):
2234        return [
2235            [orb_cos[i] for i in range(len(orb_cos)) if orb_cos[i][0] not in pop_labels]
2236            for j, orb_cos in enumerate(orbit_cosines)
2237            if j not in pop_orbits
2238        ]
2239
2240    def _get_max_cosine_labels(self, max_cosine_orbits_orig, key_points_inds_orbits, atol):
2241        max_cosine_orbits_copy = max_cosine_orbits_orig.copy()
2242        max_cosine_label_inds = np.zeros(len(max_cosine_orbits_copy))
2243        initial_max_cosine_label_inds = [max_cos_orb[0][0] for max_cos_orb in max_cosine_orbits_copy]
2244        u, inds, counts = np.unique(initial_max_cosine_label_inds, return_index=True, return_counts=True)
2245        grouped_inds = [
2246            [
2247                i
2248                for i in range(len(initial_max_cosine_label_inds))
2249                if max_cosine_orbits_copy[i][0][0] == max_cosine_orbits_copy[ind][0][0]
2250            ]
2251            for ind in inds
2252        ]
2253        pop_orbits = []
2254        pop_labels = []
2255        unassigned_orbits = []
2256        for i, ind in enumerate(inds):
2257            if counts[i] == 1:
2258                max_cosine_label_inds[ind] = initial_max_cosine_label_inds[ind]
2259                pop_orbits.append(ind)
2260                pop_labels.append(initial_max_cosine_label_inds[ind])
2261            else:
2262                next_choices = []
2263                for grouped_ind in grouped_inds[i]:
2264                    j = 1
2265                    while True:
2266                        if max_cosine_orbits_copy[grouped_ind][j][0] not in initial_max_cosine_label_inds:
2267                            next_choices.append(max_cosine_orbits_copy[grouped_ind][j][1])
2268                            break
2269                        j += 1
2270                worst_next_choice = next_choices.index(min(next_choices))
2271                for grouped_ind in grouped_inds[i]:
2272                    if grouped_ind != worst_next_choice:
2273                        unassigned_orbits.append(grouped_ind)
2274                max_cosine_label_inds[grouped_inds[i][worst_next_choice]] = initial_max_cosine_label_inds[
2275                    grouped_inds[i][worst_next_choice]
2276                ]
2277                pop_orbits.append(grouped_inds[i][worst_next_choice])
2278                pop_labels.append(initial_max_cosine_label_inds[grouped_inds[i][worst_next_choice]])
2279
2280        if len(unassigned_orbits) != 0:
2281            max_cosine_orbits_copy = self._reduce_cosines_array(max_cosine_orbits_copy, pop_orbits, pop_labels)
2282            unassigned_orbits_labels = self._get_orbit_labels(max_cosine_orbits_copy, key_points_inds_orbits, atol)
2283            for i, unassigned_orbit in enumerate(unassigned_orbits):
2284                max_cosine_label_inds[unassigned_orbit] = unassigned_orbits_labels[i]
2285
2286        return max_cosine_label_inds
2287
2288    @staticmethod
2289    def LabelPoints(index):
2290        """
2291        Axes used in generating labels for Latimer-Munro convention
2292        """
2293
2294        points = [
2295            [1, 0, 0],
2296            [0, 1, 0],
2297            [0, 0, 1],
2298            [1, 1, 0],
2299            [1, 0, 1],
2300            [0, 1, 1],
2301            [1, 1, 1],
2302            [1, 2, 0],
2303            [1, 0, 2],
2304            [1, 2, 2],
2305            [2, 1, 0],
2306            [0, 1, 2],
2307            [2, 1, 2],
2308            [2, 0, 1],
2309            [0, 2, 1],
2310            [2, 2, 1],
2311            [1, 1, 2],
2312            [1, 2, 1],
2313            [2, 1, 1],
2314            [3, 3, 2],
2315            [3, 2, 3],
2316            [2, 3, 3],
2317            [2, 2, 2],
2318            [3, 2, 2],
2319            [2, 3, 2],
2320            [1e-10, 1e-10, 1e-10],
2321        ]
2322
2323        return points[index]
2324
2325    @staticmethod
2326    def LabelSymbol(index):
2327        """
2328        Letters used in generating labels for Latimer-Munro convention
2329        """
2330
2331        symbols = [
2332            "a",
2333            "b",
2334            "c",
2335            "d",
2336            "e",
2337            "f",
2338            "g",
2339            "h",
2340            "i",
2341            "j",
2342            "k",
2343            "l",
2344            "m",
2345            "n",
2346            "o",
2347            "p",
2348            "q",
2349            "r",
2350            "s",
2351            "t",
2352            "u",
2353            "v",
2354            "w",
2355            "x",
2356            "y",
2357            "z",
2358            "Γ",
2359        ]
2360        return symbols[index]
2361