diff --git a/README.md b/README.md index 407a1dc..d5b37e1 100644 --- a/README.md +++ b/README.md @@ -20,6 +20,7 @@ If you use this tool, please cite the following work: - Y. Hinuma, G. Pizzi, Y. Kumagai, F. Oba, I. Tanaka, *Band structure diagram paths based on crystallography*, Comp. Mat. Sci. 128, 140 (2017) ([JOURNAL LINK](https://dx.doi.org/10.1016/j.commatsci.2016.10.015), [arXiv link](https://arxiv.org/abs/1602.06402)). - You should also cite [spglib](https://atztogo.github.io/spglib/) that is an essential library used in the implementation: A. Togo, I. Tanaka, "Spglib: a software library for crystal symmetry search", arXiv:1808.01590 (2018) ([spglib arXiv link](https://arxiv.org/abs/1808.01590)). +- If you run `SeeK-path` with the optional `moyopy` backend (`backend='moyopy'`), you should cite [moyo](https://github.com/spglib/moyo) in place of spglib for the symmetry analysis: K. Shinohara, *moyo: A fast and robust crystal symmetry finder, written in Rust* (2026) ([figshare link](https://figshare.com/articles/software/moyo_A_fast_and_robust_crystal_symmetry_finder_written_in_Rust_/31081162), doi:10.6084/m9.figshare.31081162.v1). ## How to install and how to use diff --git a/docs/source/maindoc.rst b/docs/source/maindoc.rst index f207c40..6ab15cb 100644 --- a/docs/source/maindoc.rst +++ b/docs/source/maindoc.rst @@ -42,7 +42,7 @@ How to use ========== The main interface of the code is the :py:func:`~seekpath.getpaths.get_path` python function:: - seekpath.get_path(structure, with_time_reversal, recipe, threshold, symprec, angle_tolerance) + seekpath.get_path(structure, with_time_reversal, recipe, threshold, symprec, angle_tolerance, backend) You need to pass a crystal structure, a boolean flag (``with_time_reversal``) to say if time-reversal symmetry is present or not, and optionally, a recipe (currently only the string ``hpkot`` is supported) and a numerical threshold. @@ -59,7 +59,112 @@ where (if ``N`` is the number of atoms): The output of the function is a dictionary containing, among other quantities, the k-vector coefficients, the suggested band path, whether the system has inversion symmetry, the crystallographic primitive lattice, the reciprocal primitive lattice. A detailed description of all output information and their format can be found in the function docstring. (Note that the ``threshold`` is the one used by seekpath to identify -e.g. the order of axes in an orthorhombic cell; instead ``symprec`` and ``angle_tolerance`` are just passed to spglib). +e.g. the order of axes in an orthorhombic cell; instead ``symprec`` and ``angle_tolerance`` are just passed to the symmetry backend). + +------------------ +Symmetry backends +------------------ + +The symmetry analysis that standardizes the input structure is done by +`spglib`_ by default. Passing ``backend='moyopy'`` uses `moyopy`_ instead, a +Rust reimplementation of spglib:: + + seekpath.get_path(structure, backend='moyopy') + +``moyopy`` is an optional dependency that requires Python >= 3.10, installed +with:: + + pip install seekpath[moyopy] + +What agrees +~~~~~~~~~~~ + +:py:func:`~seekpath.getpaths.get_path` and +:py:func:`~seekpath.getpaths.get_explicit_k_path` give the **same** result with +either backend. For all 60 reference structures shipped with seekpath, with and +without time reversal, the Bravais lattice classification, the extended Bravais +symbol, the space group, the path, the k-point labels and their coordinates are +identical, and the primitive lattices agree to within 1e-10. This also holds +under random rotations and unimodular changes of the input basis, for +left-handed input bases, and for supercells. Only the atom positions in the +standardized cells can differ (case 4). + +What can differ +~~~~~~~~~~~~~~~ + +Four cases are known. None of them makes one backend right and the other wrong, +but they are worth knowing about before switching an existing workflow. + +**1. K-points in the basis of the original cell.** +:py:func:`~seekpath.getpaths.get_path_orig_cell` and +:py:func:`~seekpath.getpaths.get_explicit_k_path_orig_cell` express the k-points +in the basis of the *input* cell. That depends on which of the several +symmetry-equivalent alignments of the standardized cell the backend picked, and +the two libraries do not always pick the same one. + +For 25 of the 60 reference structures the two backends give identical +k-points. For the other 35 the k-points differ by an integral operation in the +conventional basis, most often a 180 or 120 degree rotation, and the two paths +are symmetry-equivalent (allowing for time reversal) in every case: a band +structure computed along either is the same, but the numbers are not +identical. Code that compares ``get_path_orig_cell`` output against stored +k-points should therefore not switch backend without re-generating them. + +**2. Very loose** ``symprec`` **, right at a symmetry-promotion threshold.** +Where a structure sits on the boundary between two space groups, the two +libraries cross that boundary at slightly different tolerances, because they +apply ``symprec`` to different measures of the atomic mismatch: ``moyopy`` +promotes to the higher-symmetry group at the smaller tolerance of the two. +Displacing every atom of each reference structure by 0.02 Angstrom and raising +``symprec`` until the undistorted space group is recovered, ``spglib`` needs a +tolerance 1.1 to 2.0 times larger than ``moyopy`` does (median 1.7 over the 28 +structures where both recover it). Away from such a boundary the two agree over +the whole practical range. +If a structure's space group changes when ``symprec`` is nudged, that is a sign +the tolerance is too loose to be meaningful, whichever backend is used. + +**3. Triclinic cells with 90 degree reciprocal angles (** ``aP2`` **vs** ``aP3`` +**).** A triclinic crystal sitting on a higher-symmetry lattice - a defect cell +or a substituted supercell, say - can be classified either way, because the test +that separates the two looks at the sign of a quantity that is exactly zero. +This is **not** a backend difference: the spglib backend alone returns both +answers for the same crystal supplied in different orientations. seekpath raises +:py:class:`~seekpath.hpkot.EdgeCaseWarning` when it happens, and that warning +should be taken seriously whichever backend is in use. + +**4. Atom positions in the standardized cell.** The two libraries may choose +different, symmetry-equivalent origins, so ``conv_positions`` and +``primitive_positions`` can differ by a shift (and in atom order). This does not +affect the k-paths, which are origin-independent. + +.. note:: For a supercell whose lattice has lower symmetry than the crystal + itself, ``moyopy`` logs a warning that it omitted some symmetry operations. + seekpath does not use these operations, so the warning does not affect its + results. It can be silenced with + ``logging.getLogger('moyo').setLevel(logging.ERROR)``. + +Performance +~~~~~~~~~~~ + +``moyopy`` is the faster of the two across the board. Measured through +:py:func:`~seekpath.getpaths.get_path` on rocksalt supercells, either perfect or +with every atom displaced slightly to break the supercell symmetry: + +============================= ========== ========== +cell ``spglib`` ``moyopy`` +============================= ========== ========== +128 atoms, low symmetry 2.7 ms 1.6 ms +432 atoms, low symmetry 17.8 ms 8.4 ms +1024 atoms, low symmetry 87.3 ms 34.4 ms +128 atoms, perfect supercell 9.0 ms 1.0 ms +432 atoms, perfect supercell 10.0 ms 2.8 ms +1024 atoms, perfect supercell 14.7 ms 11.3 ms +============================= ========== ========== + +The speedup grows with the number of atoms for a low-symmetry cell, and shrinks +for a large perfect supercell, where ``spglib`` scales well because it reduces +to the primitive cell early. + ---------------------------------------- K-point path for non-standard unit cells @@ -127,6 +232,8 @@ https://aiida-core.readthedocs.io/en/latest/datatypes/kpoints.html .. _JOURNAL LINK: https://dx.doi.org/10.1016/j.commatsci.2016.10.015 .. _arXiv link: https://arxiv.org/abs/1602.06402 .. _spglib: https://atztogo.github.io/spglib/ + +.. _moyopy: https://spglib.github.io/moyo/python/ .. _Materials Cloud: https://www.materialscloud.org/tools/seekpath/ .. _docker hub: https://hub.docker.com/r/giovannipizzi/seekpath/ .. _AiiDA: https://www.aiida.net diff --git a/pyproject.toml b/pyproject.toml index a6dcdf9..423e8cd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,8 +43,12 @@ Downloads = "https://github.com/materialscloud-org/seekpath/archive/v2.2.2.tar.g bz = [ "scipy>=1", ] +moyopy = [ + "moyopy>=0.21; python_version >= '3.10'", +] dev = [ "black==23.3.0", + "moyopy>=0.21; python_version >= '3.10'", "pre-commit~=3.5", "prospector==1.11.0", "pytest==7.3.1", diff --git a/seekpath/getpaths.py b/seekpath/getpaths.py index f9a8ebc..3781e0d 100644 --- a/seekpath/getpaths.py +++ b/seekpath/getpaths.py @@ -5,6 +5,7 @@ import numpy as np import warnings from . import SupercellWarning +from .hpkot import DEFAULT_BACKEND def get_explicit_from_implicit(seekpath_output, reference_distance): @@ -82,6 +83,7 @@ def get_path( threshold=1.0e-7, symprec=1e-05, angle_tolerance=-1.0, + backend=DEFAULT_BACKEND, ): r""" Return the kpoint path information for band structure given a @@ -125,9 +127,16 @@ def get_path( Note that depending on the bravais lattice, the meaning of the threshold is different (angle, length, ...) - :param symprec: the symmetry precision used internally by SPGLIB + :param symprec: the symmetry precision used internally by the symmetry backend - :param angle_tolerance: the angle_tolerance used internally by SPGLIB + :param angle_tolerance: the angle_tolerance used internally by the symmetry + backend + + :param backend: the symmetry backend used to standardize the structure, + either ``'spglib'`` (the default) or ``'moyopy'``. The ``'moyopy'`` + backend requires the optional ``moyopy`` dependency and is usually + faster; see the documentation for where the results of the + two backends can differ. :return: a dictionary with the following @@ -185,6 +194,7 @@ def get_path( threshold=threshold, symprec=symprec, angle_tolerance=angle_tolerance, + backend=backend, ) else: @@ -203,6 +213,7 @@ def get_explicit_k_path( threshold=1.0e-7, symprec=1e-05, angle_tolerance=-1.0, + backend=DEFAULT_BACKEND, ): r""" Return the kpoint path for band structure (in scaled and absolute @@ -257,9 +268,16 @@ def get_explicit_k_path( Note that depending on the bravais lattice, the meaning of the threshold is different (angle, length, ...) - :param symprec: the symmetry precision used internally by SPGLIB + :param symprec: the symmetry precision used internally by the symmetry backend + + :param angle_tolerance: the angle_tolerance used internally by the symmetry + backend - :param angle_tolerance: the angle_tolerance used internally by SPGLIB + :param backend: the symmetry backend used to standardize the structure, + either ``'spglib'`` (the default) or ``'moyopy'``. The ``'moyopy'`` + backend requires the optional ``moyopy`` dependency and is usually + faster; see the documentation for where the results of the + two backends can differ. .. versionchanged:: 1.8 The key ``segments`` has been renamed ``explicit_segments`` @@ -313,6 +331,7 @@ def get_explicit_k_path( threshold=threshold, symprec=symprec, angle_tolerance=angle_tolerance, + backend=backend, ) else: @@ -336,6 +355,7 @@ def get_path_orig_cell( threshold=1.0e-7, symprec=1e-05, angle_tolerance=-1.0, + backend=DEFAULT_BACKEND, ): r""" Return the kpoint path information for band structure given a @@ -388,9 +408,16 @@ def get_path_orig_cell( Note that depending on the bravais lattice, the meaning of the threshold is different (angle, length, ...) - :param symprec: the symmetry precision used internally by SPGLIB + :param symprec: the symmetry precision used internally by the symmetry backend - :param angle_tolerance: the angle_tolerance used internally by SPGLIB + :param angle_tolerance: the angle_tolerance used internally by the symmetry + backend + + :param backend: the symmetry backend used to standardize the structure, + either ``'spglib'`` (the default) or ``'moyopy'``. The ``'moyopy'`` + backend requires the optional ``moyopy`` dependency and is usually + faster; see the documentation for where the results of the + two backends can differ. :return: a dictionary with the following @@ -420,6 +447,7 @@ def get_path_orig_cell( symprec=symprec, angle_tolerance=angle_tolerance, recipe=recipe, + backend=backend, ) # The volume ratio is negative for a left-handed input cell, so only its @@ -480,6 +508,7 @@ def get_explicit_k_path_orig_cell( threshold=1.0e-7, symprec=1e-05, angle_tolerance=-1.0, + backend=DEFAULT_BACKEND, ): r""" Return the kpoint path for band structure (in scaled and absolute @@ -537,9 +566,16 @@ def get_explicit_k_path_orig_cell( Note that depending on the bravais lattice, the meaning of the threshold is different (angle, length, ...) - :param symprec: the symmetry precision used internally by SPGLIB + :param symprec: the symmetry precision used internally by the symmetry backend + + :param angle_tolerance: the angle_tolerance used internally by the symmetry + backend - :param angle_tolerance: the angle_tolerance used internally by SPGLIB + :param backend: the symmetry backend used to standardize the structure, + either ``'spglib'`` (the default) or ``'moyopy'``. The ``'moyopy'`` + backend requires the optional ``moyopy`` dependency and is usually + faster; see the documentation for where the results of the + two backends can differ. .. versionchanged:: 1.8 The key ``segments`` has been renamed ``explicit_segments`` @@ -593,6 +629,7 @@ def get_explicit_k_path_orig_cell( symprec=symprec, angle_tolerance=angle_tolerance, recipe=recipe, + backend=backend, ) # Set reciprocal_primitive_lattice as the reciprocal lattice of the original diff --git a/seekpath/hpkot/__init__.py b/seekpath/hpkot/__init__.py index e6930c2..834c472 100644 --- a/seekpath/hpkot/__init__.py +++ b/seekpath/hpkot/__init__.py @@ -13,6 +13,11 @@ the Materials Project (https://materialsproject.org). """ +from .backends import ( # noqa: F401 + DEFAULT_BACKEND, + SymmetryDetectionError, +) + class EdgeCaseWarning(RuntimeWarning): """ @@ -21,18 +26,13 @@ class EdgeCaseWarning(RuntimeWarning): """ -class SymmetryDetectionError(Exception): - """ - Error raised if spglib could not detect the symmetry. - """ - - def get_path( structure, with_time_reversal=True, threshold=1.0e-7, symprec=1e-05, angle_tolerance=-1.0, + backend=DEFAULT_BACKEND, ): r""" Return the kpoint path information for band structure given a @@ -76,9 +76,16 @@ def get_path( Note that depending on the bravais lattice, the meaning of the threshold is different (angle, length, ...) - :param symprec: the symmetry precision used internally by SPGLIB + :param symprec: the symmetry precision used internally by the symmetry backend - :param angle_tolerance: the angle_tolerance used internally by SPGLIB + :param angle_tolerance: the angle_tolerance used internally by the symmetry + backend + + :param backend: the symmetry backend used to standardize the structure, + either ``'spglib'`` (the default) or ``'moyopy'``. The ``'moyopy'`` + backend requires the optional ``moyopy`` dependency and is usually + faster; see the documentation for where the results of the + two backends can differ. :return: a dictionary with the following @@ -134,51 +141,40 @@ def get_path( import warnings import numpy as np + import spglib from .tools import ( - check_spglib_version, extend_kparam, eval_expr, eval_expr_simple, get_cell_params, - get_dot_access_dataset, get_path_data, get_reciprocal_cell_rows, get_real_cell_from_reciprocal_rows, ) + from .backends import get_symmetry_dataset from .spg_mapping import get_spgroup_data, get_primitive - # I check if the SPGlib version is recent enough (raises ValueError) - # otherwise - spglib = check_spglib_version() - structure_internal = ( np.array(structure[0]), np.array(structure[1]), np.array(structure[2]), ) - # Symmetry analysis by SPGlib, get crystallographic lattice, + # Symmetry analysis by the chosen backend, get crystallographic lattice, # and cell parameters for this lattice - dataset = get_dot_access_dataset( - spglib.get_symmetry_dataset( - structure_internal, symprec=symprec, angle_tolerance=angle_tolerance - ) + dataset = get_symmetry_dataset( + structure_internal, + symprec=symprec, + angle_tolerance=angle_tolerance, + backend=backend, ) - if dataset is None: - raise SymmetryDetectionError( - 'Spglib could not detect the symmetry of the system' - ) conv_lattice = dataset.std_lattice conv_positions = dataset.std_positions conv_types = dataset.std_types a, b, c, cosalpha, cosbeta, cosgamma = get_cell_params(conv_lattice) spgrp_num = dataset.number - # This is the transformation from the original to the crystallographic - # conventional (called std in spglib) - # Lattice^{crystallographic_bravais} = L^{original} * transf_matrix - transf_matrix = dataset.transformation_matrix - volume_conv_wrt_original = np.linalg.det(transf_matrix) + volume_original_wrt_conv = dataset.volume_original_wrt_conv # Get the properties of the spacegroup, needed to get the bravais_lattice properties = get_spgroup_data()[spgrp_num] @@ -502,8 +498,8 @@ def get_path( # For the time being disabled, not valid for aP lattices # (for which we would need the transformation matrix from niggli) #'transformation_matrix': transf_matrix, - 'volume_original_wrt_conv': volume_conv_wrt_original, - 'volume_original_wrt_prim': volume_conv_wrt_original * np.linalg.det(invP), + 'volume_original_wrt_conv': volume_original_wrt_conv, + 'volume_original_wrt_prim': volume_original_wrt_conv * np.linalg.det(invP), 'spacegroup_number': dataset.number, 'spacegroup_international': dataset.international, 'rotation_matrix': dataset.std_rotation_matrix, diff --git a/seekpath/hpkot/backends.py b/seekpath/hpkot/backends.py new file mode 100644 index 0000000..d01141b --- /dev/null +++ b/seekpath/hpkot/backends.py @@ -0,0 +1,185 @@ +"""Symmetry backends used to standardize a crystal structure. + +seekpath needs only a small part of a full symmetry dataset: the conventional +(standardized) cell, the space-group number and international symbol, the +rotation from the input orientation to the standardized one, and the volume +ratio between the input and the conventional cell. This module normalizes what +the supported symmetry libraries return into a single :class:`SymmetryDataset`, +so that everything downstream in :mod:`seekpath.hpkot` stays backend agnostic. + +Two backends are supported: + +``spglib`` + The default, and a hard dependency of seekpath. + +``moyopy`` + An optional backend. moyopy is a Rust reimplementation of spglib and is + usually faster. It is queried with ``Setting.spglib()``, which selects the + same Hall-symbol representative that spglib uses. +""" + +from dataclasses import dataclass +from functools import lru_cache + +import numpy as np + +DEFAULT_BACKEND = 'spglib' +MOYOPY_MIN_VERSION = '0.21' + + +class SymmetryDetectionError(Exception): + """Error raised if the symmetry of the structure could not be detected.""" + + +@dataclass +class SymmetryDataset: + """The subset of a symmetry dataset that seekpath uses. + + :param number: the international space-group number. + :param international: the international (Hermann-Mauguin) symbol, without + spaces. + :param std_lattice: the conventional cell, as a 3x3 array with the lattice + vectors as *rows*. + :param std_positions: fractional coordinates of the atoms in the + conventional cell. + :param std_types: atomic numbers of the atoms in the conventional cell. + :param std_rotation_matrix: the rotation in Cartesian space bringing the + input structure onto the standardized one. + :param volume_original_wrt_conv: the volume of the input cell divided by + the volume of the conventional cell. + """ + + number: int + international: str + std_lattice: np.ndarray + std_positions: np.ndarray + std_types: np.ndarray + std_rotation_matrix: np.ndarray + volume_original_wrt_conv: float + + +def get_symmetry_dataset( + structure, symprec=1e-05, angle_tolerance=-1.0, backend=DEFAULT_BACKEND +): + """Standardize ``structure`` with the requested backend. + + :param structure: a tuple ``(cell, positions, numbers)``, with the lattice + vectors as rows of ``cell`` and ``positions`` in fractional + coordinates. + :param symprec: the symmetry precision. + :param angle_tolerance: the angle tolerance. The spglib convention of a + negative value meaning "use the default" is honoured by both backends. + :param backend: one of :data:`SUPPORTED_BACKENDS`. + + :return: a :class:`SymmetryDataset`. + + :raise ValueError: if ``backend`` is not a supported backend. + :raise SymmetryDetectionError: if the backend could not detect the + symmetry of the structure. + """ + try: + get_dataset = _BACKENDS[backend] + except KeyError: + raise ValueError( + f"Unknown symmetry backend '{backend}', should be one of {SUPPORTED_BACKENDS}" + ) from None + return get_dataset(structure, symprec, angle_tolerance) + + +def _get_spglib_dataset(structure, symprec, angle_tolerance): + """Build a :class:`SymmetryDataset` with spglib.""" + from .tools import check_spglib_version, get_dot_access_dataset + + spglib = check_spglib_version() + + dataset = get_dot_access_dataset( + spglib.get_symmetry_dataset( + structure, symprec=symprec, angle_tolerance=angle_tolerance + ) + ) + if dataset is None: + raise SymmetryDetectionError( + 'spglib could not detect the symmetry of the system' + ) + + return SymmetryDataset( + number=dataset.number, + international=dataset.international, + std_lattice=np.array(dataset.std_lattice), + std_positions=np.array(dataset.std_positions), + std_types=np.array(dataset.std_types), + std_rotation_matrix=np.array(dataset.std_rotation_matrix), + volume_original_wrt_conv=np.linalg.det(dataset.transformation_matrix), + ) + + +def _get_moyopy_dataset(structure, symprec, angle_tolerance): + """Build a :class:`SymmetryDataset` with moyopy.""" + moyopy = _check_moyopy_version() + + cell, positions, numbers = structure + moyo_cell = moyopy.Cell( + np.array(cell, dtype=float).tolist(), + np.array(positions, dtype=float).tolist(), + np.array(numbers).tolist(), + ) + + try: + dataset = moyopy.MoyoDataset( + moyo_cell, + symprec=symprec, + angle_tolerance=( + np.deg2rad(angle_tolerance) if angle_tolerance > 0 else None + ), + setting=moyopy.Setting.spglib(), + ) + except Exception as exc: + raise SymmetryDetectionError( + f'moyopy could not detect the symmetry of the system: {exc}' + ) from exc + + international = moyopy.SpaceGroupType(dataset.number).hm_short.replace(' ', '') + + volume_original_wrt_conv = 1 / np.linalg.det(np.array(dataset.std_linear)) + + return SymmetryDataset( + number=dataset.number, + international=international, + std_lattice=np.array(dataset.std_cell.basis), + std_positions=np.array(dataset.std_cell.positions), + std_types=np.array(dataset.std_cell.numbers), + std_rotation_matrix=np.array(dataset.std_rotation_matrix), + volume_original_wrt_conv=volume_original_wrt_conv, + ) + + +@lru_cache(maxsize=None) +def _check_moyopy_version(): + """Import moyopy, checking it is recent enough. + + :return: the moyopy module. + + :raise ValueError: if moyopy is missing or older than + :data:`MOYOPY_MIN_VERSION`. + """ + from importlib.metadata import version + + from packaging.version import Version + + try: + import moyopy + except ImportError as exc: + raise ValueError( + f'moyopy >= {MOYOPY_MIN_VERSION} is required for the ' + "'moyopy' backend, but it could not be imported. Install it " + 'with `pip install seekpath[moyopy]` (requires Python >= 3.10)' + ) from exc + + if Version(version('moyopy')) < Version(MOYOPY_MIN_VERSION): + raise ValueError(f'Invalid moyopy version, need >= {MOYOPY_MIN_VERSION}') + + return moyopy + + +_BACKENDS = {'spglib': _get_spglib_dataset, 'moyopy': _get_moyopy_dataset} +SUPPORTED_BACKENDS = tuple(_BACKENDS) diff --git a/tests/test_backends.py b/tests/test_backends.py new file mode 100644 index 0000000..f64251f --- /dev/null +++ b/tests/test_backends.py @@ -0,0 +1,250 @@ +"""Test the symmetry backends, and that they agree with each other. + +The ``moyopy`` backend is optional, so every test that needs it is skipped when +it is not installed. +""" + +import glob +import os + +import numpy as np +import pytest +from test_paths_hpkot import simple_read_poscar + +import seekpath +from seekpath import hpkot +from seekpath.hpkot.backends import ( + DEFAULT_BACKEND, + SymmetryDetectionError, + get_symmetry_dataset, +) + +BAND_PATH_DATA = os.path.join(os.path.dirname(hpkot.__file__), 'band_path_data') + +try: + import moyopy as _moyopy # noqa: F401 + + HAS_MOYOPY = True +except ImportError: + HAS_MOYOPY = False + +needs_moyopy = pytest.mark.skipif(not HAS_MOYOPY, reason='moyopy is not installed') + +REFERENCE_STRUCTURES = sorted( + (os.path.basename(folder), os.path.basename(poscar).replace('POSCAR_', '')) + for folder in glob.glob(os.path.join(BAND_PATH_DATA, '*')) + if os.path.isdir(folder) + for poscar in glob.glob(os.path.join(folder, 'POSCAR_*')) +) + + +def read_reference(ext_bravais, variant): + """Read the reference POSCAR for an extended Bravais symbol.""" + return simple_read_poscar( + os.path.join(BAND_PATH_DATA, ext_bravais, f'POSCAR_{variant}') + ) + + +def ids(param): + """Readable test ids for the (ext_bravais, variant) pairs.""" + return f'{param[0]}-{param[1]}' + + +class TestBackendSelection: + """Test how the backend is selected and validated.""" + + def test_default_is_spglib(self): + """The default backend must stay spglib, for backwards compatibility.""" + assert DEFAULT_BACKEND == 'spglib' + + def test_unknown_backend_raises(self): + """An unknown backend name is rejected with a helpful message.""" + structure = read_reference('cP1', 'inversion') + with pytest.raises(ValueError, match="Unknown symmetry backend 'nope'"): + hpkot.get_path(structure, backend='nope') + + @pytest.mark.parametrize( + 'backend', ['spglib', pytest.param('moyopy', marks=needs_moyopy)] + ) + def test_symmetry_detection_error(self, backend): + """Both backends report a failed detection the same way. + + A degenerate cell (two parallel lattice vectors, so zero volume) cannot + be analysed, and must raise ``SymmetryDetectionError`` whichever backend + is used, rather than the backend's own exception type. + """ + degenerate = ( + [[4.0, 0.0, 0.0], [4.0, 0.0, 0.0], [0.0, 0.0, 4.0]], + [[0.0, 0.0, 0.0]], + [6], + ) + with pytest.raises(SymmetryDetectionError): + get_symmetry_dataset(degenerate, backend=backend) + + +@needs_moyopy +class TestMoyopyDataset: + """Test the moyopy dataset against the spglib one, field by field.""" + + @pytest.mark.parametrize('case', REFERENCE_STRUCTURES, ids=ids) + def test_dataset_fields_agree(self, case): + """The standardized cell, space group and volume ratio must agree. + + The atom positions are deliberately not compared: moyo may pick a + different, symmetry-equivalent origin, which is a legitimate choice and + does not affect the k-path. + """ + structure = read_reference(*case) + spglib_ds = get_symmetry_dataset(structure, backend='spglib') + moyopy_ds = get_symmetry_dataset(structure, backend='moyopy') + + assert spglib_ds.number == moyopy_ds.number + assert spglib_ds.international == moyopy_ds.international + assert sorted(spglib_ds.std_types) == sorted(moyopy_ds.std_types) + np.testing.assert_allclose( + spglib_ds.std_lattice, moyopy_ds.std_lattice, atol=1e-8 + ) + np.testing.assert_allclose( + spglib_ds.volume_original_wrt_conv, + moyopy_ds.volume_original_wrt_conv, + atol=1e-8, + ) + + @pytest.mark.parametrize(('degrees', 'expected'), [(-1.0, None), (5.0, 0.0872665)]) + def test_angle_tolerance_is_converted(self, monkeypatch, degrees, expected): + """seekpath's angle tolerance is in degrees (as spglib's), moyopy's in radians.""" + import moyopy + + received = {} + original = moyopy.MoyoDataset + + def spy(*args, **kwargs): + received.update(kwargs) + return original(*args, **kwargs) + + monkeypatch.setattr(moyopy, 'MoyoDataset', spy) + get_symmetry_dataset( + read_reference('cP1', 'inversion'), + angle_tolerance=degrees, + backend='moyopy', + ) + if expected is None: + assert received['angle_tolerance'] is None + else: + assert received['angle_tolerance'] == pytest.approx(expected) + + +@needs_moyopy +class TestBackendsAgree: + """Test that both backends give the same k-path for the reference set.""" + + @pytest.mark.parametrize('with_time_reversal', [True, False]) + @pytest.mark.parametrize('case', REFERENCE_STRUCTURES, ids=ids) + def test_get_path_agrees(self, case, with_time_reversal): + """``get_path`` must give identical results with either backend.""" + structure = read_reference(*case) + kwargs = {'with_time_reversal': with_time_reversal} + res_spglib = hpkot.get_path(structure, backend='spglib', **kwargs) + res_moyopy = hpkot.get_path(structure, backend='moyopy', **kwargs) + + for key in ( + 'bravais_lattice', + 'bravais_lattice_extended', + 'spacegroup_number', + 'spacegroup_international', + 'has_inversion_symmetry', + 'augmented_path', + 'path', + ): + assert res_spglib[key] == res_moyopy[key], f'{key} differs' + + assert set(res_spglib['point_coords']) == set(res_moyopy['point_coords']) + for label, coords in res_spglib['point_coords'].items(): + np.testing.assert_allclose( + coords, res_moyopy['point_coords'][label], atol=1e-6, err_msg=label + ) + + for key in ('primitive_lattice', 'reciprocal_primitive_lattice'): + np.testing.assert_allclose( + res_spglib[key], res_moyopy[key], atol=1e-8, err_msg=key + ) + + @pytest.mark.parametrize('case', REFERENCE_STRUCTURES, ids=ids) + def test_reference_classification(self, case): + """moyopy must reproduce the HPKOT extended Bravais classification. + + This is the same assertion the spglib-only tests make, so it checks + moyopy against the reference data rather than only against spglib. + """ + ext_bravais, variant = case + structure = read_reference(ext_bravais, variant) + res = hpkot.get_path(structure, with_time_reversal=False, backend='moyopy') + + assert res['bravais_lattice_extended'] == ext_bravais + assert res['has_inversion_symmetry'] == variant.startswith('inversion') + + @pytest.mark.parametrize('case', REFERENCE_STRUCTURES, ids=ids) + def test_explicit_k_path_agrees(self, case): + """The explicit (interpolated) k-point list must agree as well.""" + structure = read_reference(*case) + res_spglib = seekpath.get_explicit_k_path(structure, backend='spglib') + res_moyopy = seekpath.get_explicit_k_path(structure, backend='moyopy') + + assert ( + res_spglib['explicit_kpoints_labels'] + == res_moyopy['explicit_kpoints_labels'] + ) + np.testing.assert_allclose( + res_spglib['explicit_kpoints_rel'], + res_moyopy['explicit_kpoints_rel'], + atol=1e-6, + ) + + +@needs_moyopy +class TestOrigCellPathsAreEquivalent: + """Test the k-path expressed in the basis of the *original* cell. + + Unlike ``get_path``, ``get_path_orig_cell`` depends on + ``std_rotation_matrix``, i.e. on which of the symmetry-equivalent + standardizations the backend picked. spglib and moyopy may pick different + ones, so the two paths need not be numerically identical - but they must be + related by a single symmetry operation of the crystal, which makes them + physically equivalent. + """ + + @staticmethod + def _symmetry_equivalent(points_a, points_b, rotations): + """Is there one operation mapping every k-point of a onto b?""" + labels = sorted(points_a) + k_a = np.array([points_a[label] for label in labels]) + k_b = np.array([points_b[label] for label in labels]) + + candidates = [np.array(rot, dtype=float) for rot in rotations] + # Time reversal (k -> -k) is a symmetry of the band structure too + candidates += [-candidate for candidate in candidates] + + for candidate in candidates: + for matrix in (candidate, candidate.T): + difference = k_a @ matrix - k_b + # Two k-points differing by a reciprocal lattice vector are equal + difference -= np.round(difference) + if np.abs(difference).max() < 1e-5: + return True + return False + + @pytest.mark.parametrize('case', REFERENCE_STRUCTURES, ids=ids) + def test_orig_cell_path_equivalent(self, case): + """The original-cell k-points must be symmetry-equivalent.""" + import spglib + + structure = read_reference(*case) + res_spglib = seekpath.get_path_orig_cell(structure, backend='spglib') + res_moyopy = seekpath.get_path_orig_cell(structure, backend='moyopy') + + assert res_spglib['path'] == res_moyopy['path'] + + rotations = spglib.get_symmetry_dataset(structure, symprec=1e-5).rotations + assert self._symmetry_equivalent( + res_spglib['point_coords'], res_moyopy['point_coords'], rotations + )