Repository: symfc/symfc (language: python)

Issue:
Preparing for release v1.0.0
This is PR toward release v1.0.0.

* Migration from setup.py and setup.cfg to pyproject.toml
* `src/symfc/version.py` contains `__version__ = "1.0.0"`


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 68bc17f..e6ef6f2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -158,3 +158,5 @@ cython_debug/
 #  and can be added to the global gitignore or merged into this file.  For a more nuclear
 #  option (not recommended) you can uncomment the following to ignore the entire idea folder.
 #.idea/
+
+.vscode
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 0000000..f7d3467
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,41 @@
+# See https://pre-commit.com for more informatio
+# See https://pre-commit.com/hooks.html for more hooks
+repos:
+- repo: https://github.com/pre-commit/pre-commit-hooks
+  rev: v4.5.0
+  hooks:
+  - id: trailing-whitespace
+  - id: end-of-file-fixer
+  - id: check-yaml
+    exclude: ^conda/
+  - id: check-added-large-files
+
+- repo: https://github.com/pycqa/flake8
+  rev: 7.0.0
+  hooks:
+  - id: flake8
+    args:
+    - "--max-line-length=88"
+    - "--ignore=E203,W503"
+  - id: flake8
+    args:
+    - "--max-line-length=88"
+    - "--ignore=E203,W503,E501"
+
+- repo: https://github.com/psf/black
+  rev: 24.1.1
+  hooks:
+  - id: black
+    args:
+    - --line-length=88
+
+- repo: https://github.com/pycqa/pydocstyle
+  rev: 6.3.0
+  hooks:
+  - id: pydocstyle
+
+- repo: https://github.com/pycqa/isort
+  rev: 5.13.2
+  hooks:
+    - id: isort
+      name: isort (python)
diff --git a/LICENSE b/LICENSE
index 0c36bd9..050c9fd 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 BSD 3-Clause License
 
-Copyright (c) 2023, phonopy
+Copyright (c) 2023, symfc project
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
diff --git a/README.md b/README.md
index 9e9236e..a6524d1 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,22 @@
 # symfc
-Generate symmetrized force constants
+
+## What does symfc do?
+
+Atomic vibrations in crystals are often conveniently described using the phonon
+model. In this model, the crystal potential is expanded into a Taylor series
+with respect to atomic displacements from their equilibrium positions, and the
+expansion coefficients are referred to as force constants.
+
+Predicting phonon properties through computer simulations is becoming
+increasingly popular, with the supercell approach being one of the techniques
+employed for phonon calculations. In this method, force constants are derived
+from datasets of atomic forces and displacements obtained from supercell
+snapshots, which feature various configurations of atomic displacements.
+
+While force constants possess specific symmetries, those computed from
+displacement-force datasets often do not adhere to these symmetries due to
+factors such as numerical noise or approximations used. Symfc is a software
+designed to compute force constants from displacement-force datasets in the
+supercell approach, ensuring they meet the required symmetry constraints.
+
+For more detail, please refer https://arxiv.org/abs/2403.03588.
diff --git a/examples/NaCl444-rd/README.md b/examples/NaCl444-rd/README.md
new file mode 100644
index 0000000..215e1e5
--- /dev/null
+++ b/examples/NaCl444-rd/README.md
@@ -0,0 +1,8 @@
+# NaCl example
+
+The python example is found in `calc-fc.py`.
+
+Force constants are calculated from 20 supercell snapshots. The crystal
+structure, displacement-force dataset, and parameters needed for non-analytical
+term correction are found in `phonopy_NaCl444_rd.yaml.xz` in the phonopy-yaml
+format.
diff --git a/examples/NaCl444-rd/calc-fc.py b/examples/NaCl444-rd/calc-fc.py
new file mode 100644
index 0000000..8338362
--- /dev/null
+++ b/examples/NaCl444-rd/calc-fc.py
@@ -0,0 +1,18 @@
+"""Calculate force constants using symfc."""
+
+import phonopy
+
+from symfc import Symfc
+
+ph = phonopy.load("phonopy_NaCl444_rd.yaml.xz", produce_fc=False)
+symfc = Symfc(
+    ph.supercell,
+    displacements=ph.dataset["displacements"],
+    forces=ph.dataset["forces"],
+    orders=[
+        2,
+    ],
+    log_level=1,
+)
+ph.force_constants = symfc.force_constants[2]
+ph.auto_band_structure(plot=True).show()
diff --git a/examples/NaCl444-rd/phonopy_NaCl444_rd.yaml.xz b/examples/NaCl444-rd/phonopy_NaCl444_rd.yaml.xz
new file mode 100644
index 0000000..be4506c
Binary files /dev/null and b/examples/NaCl444-rd/phonopy_NaCl444_rd.yaml.xz differ
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..9db6369
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,38 @@
+[build-system]
+requires = ["setuptools >= 61.0"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "symfc"
+dynamic = ["version"]
+readme = {file = "README.md", content-type = "text/markdown"}
+description = "This is the symfc module."
+authors = [
+    { name = "Atsuto Seko", email = "seko@cms.mtl.kyoto-u.ac.jp" }
+]
+maintainers = [
+    {name = "Atsushi Togo", email = "atztogo@gmail.com"}
+]
+requires-python = ">=3.8"
+dependencies = [
+    "numpy",
+    "scipy",
+    "spglib",
+]
+license = {file = "LICENSE"}
+
+[project.urls]
+"Homepage" = "https://github.com/symfc/symfc"
+
+[tool.setuptools.dynamic]
+version = { attr = "symfc.version.__version__" }
+
+[tool.flake8]
+max-line-length = 88
+extend-ignore = "E203,W503"
+
+[tool.black]
+line-length = 88
+
+[tool.isort]
+profile = "black"
diff --git a/setup.cfg.bak b/setup.cfg.bak
new file mode 100644
index 0000000..95e040d
--- /dev/null
+++ b/setup.cfg.bak
@@ -0,0 +1,10 @@
+[options]
+package_dir=
+    =src
+packages=find:
+
+[options.packages.find]
+where=.
+
+[metadata]
+version = attr: symfc.version.__version__
diff --git a/setup.py.bak b/setup.py.bak
new file mode 100644
index 0000000..9caa8a1
--- /dev/null
+++ b/setup.py.bak
@@ -0,0 +1,24 @@
+"""Setup script."""
+
+from setuptools import setup
+from pathlib import Path
+
+cwd = Path(__file__).parent
+
+# with open("src/symfc/version.py") as f:
+#     for line in f:
+#         if "__version__" in line:
+#             version = line.split("=")[1].strip().strip('"')
+
+setup(
+    name="symfc",
+    # version=version,
+    description="This is the symfc module.",
+    author="Atsuto Seko",
+    author_email="seko@cms.mtl.kyoto-u.ac.jp",
+    url="https://github.com/symfc/symfc",
+    install_requires=["numpy", "scipy", "spglib"],
+    provides=["symfc"],
+    platforms=["all"],
+    python_requires=">=3.8",
+)
diff --git a/src/symfc/__init__.py b/src/symfc/__init__.py
new file mode 100644
index 0000000..dd41795
--- /dev/null
+++ b/src/symfc/__init__.py
@@ -0,0 +1,4 @@
+"""Forcd constants calculation code: Symfc."""
+
+from .api_symfc import Symfc  # noqa F401
+from .version import __version__  # noqa F401
diff --git a/src/symfc/api_symfc.py b/src/symfc/api_symfc.py
new file mode 100644
index 0000000..3fb5117
--- /dev/null
+++ b/src/symfc/api_symfc.py
@@ -0,0 +1,180 @@
+"""Symfc API."""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+from typing import Optional, Union
+
+import numpy as np
+
+from symfc.basis_sets import FCBasisSetBase, FCBasisSetO2, FCBasisSetO3
+from symfc.solvers import FCSolverO2, FCSolverO2O3
+from symfc.utils.utils import SymfcAtoms
+
+
+class Symfc:
+    """Symfc API."""
+
+    def __init__(
+        self,
+        supercell: SymfcAtoms,
+        displacements: Optional[np.ndarray] = None,
+        forces: Optional[np.ndarray] = None,
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method."""
+        self._supercell: SymfcAtoms = supercell
+        self._displacements: Optional[np.ndarray] = displacements
+        self._forces: Optional[np.ndarray] = forces
+        self._use_mkl = use_mkl
+        self._log_level = log_level
+
+        self._basis_set: dict[FCBasisSetBase] = {}
+        self._force_constants: dict[np.ndarray] = {}
+
+    @property
+    def basis_set(self) -> dict[FCBasisSetBase]:
+        """Return basis set instance.
+
+        Returns
+        -------
+        dict[FCBasisSet]
+            The key is the order of basis set in int.
+
+        """
+        return self._basis_set
+
+    @property
+    def force_constants(self) -> dict[np.ndarray]:
+        """Return force constants.
+
+        Returns
+        -------
+        dict[np.ndarray]
+            The key is the order of force_constants in int.
+
+        """
+        return self._force_constants
+
+    @property
+    def displacements(self) -> np.ndarray:
+        """Setter and getter of supercell displacements.
+
+        ndarray
+            shape=(n_snapshot, natom, 3), dtype='double', order='C'
+
+        """
+        return self._displacements
+
+    @displacements.setter
+    def displacements(self, displacements: Union[np.ndarray, list, tuple]):
+        self._displacements = np.array(displacements, dtype="double", order="C")
+
+    @property
+    def forces(self) -> np.ndarray:
+        """Setter and getter of supercell forces.
+
+        ndarray
+            shape=(n_snapshot, natom, 3), dtype='double', order='C'
+
+        """
+        return self._forces
+
+    @forces.setter
+    def forces(self, forces: Union[np.ndarray, list, tuple]):
+        self._forces = np.array(forces, dtype="double", order="C")
+
+    def run(self, orders: Sequence[int], is_compact_fc=True) -> Symfc:
+        """Run basis set and force constants calculation."""
+        if (
+            orders is not None
+            and self._displacements is not None
+            and self._forces is not None
+        ):
+            for order in orders:
+                self._compute_basis_set(order)
+            self.solve(orders, is_compact_fc=is_compact_fc)
+        return self
+
+    def solve(self, orders: Sequence[int], is_compact_fc=True) -> Symfc:
+        """Calculate force constants.
+
+        orders : Sequence[int]
+            Sequence of fc orders.
+
+        """
+        self._check_dataset()
+        for order in orders:
+            if order not in (2, 3):
+                raise NotImplementedError("Only order-2 and order-3 are implemented.")
+            if order == 2:
+                basis_set: FCBasisSetO2 = self._basis_set[2]
+                solver_o2 = FCSolverO2(
+                    basis_set,
+                    use_mkl=self._use_mkl,
+                    log_level=self._log_level,
+                ).solve(self._displacements, self._forces)
+                if is_compact_fc:
+                    self._force_constants[2] = solver_o2.compact_fc
+                else:
+                    self._force_constants[2] = solver_o2.full_fc
+            if order == 3:
+                basis_set_o2: FCBasisSetO2 = self._basis_set[2]
+                basis_set_o3: FCBasisSetO3 = self._basis_set[3]
+                solver_o2o3 = FCSolverO2O3(
+                    [basis_set_o2, basis_set_o3],
+                    use_mkl=self._use_mkl,
+                    log_level=self._log_level,
+                ).solve(self._displacements, self._forces)
+                if is_compact_fc:
+                    fc2, fc3 = solver_o2o3.compact_fc
+                else:
+                    fc2, fc3 = solver_o2o3.full_fc
+                self._force_constants[2] = fc2
+                self._force_constants[3] = fc3
+        return self
+
+    def _compute_basis_set(self, order: int):
+        """Set order of force constants."""
+        if order not in (2, 3):
+            raise NotImplementedError("Only fc2 basis set is implemented.")
+
+        if order == 2:
+            basis_set_o2 = FCBasisSetO2(
+                self._supercell, use_mkl=self._use_mkl, log_level=self._log_level
+            ).run()
+            self._basis_set[2] = basis_set_o2
+        if order == 3:
+            basis_set_o3 = FCBasisSetO3(
+                self._supercell, use_mkl=self._use_mkl, log_level=self._log_level
+            ).run()
+            self._basis_set[3] = basis_set_o3
+
+    def _check_dataset(self):
+        if self._displacements is None:
+            raise RuntimeError("Dispalcements not found.")
+        if self._forces is None:
+            raise RuntimeError("Forces not found.")
+        if self._displacements.shape != self._forces.shape:
+            raise RuntimeError("Shape mismatch between dispalcements and forces.")
+        if self._displacements.shape != self._forces.shape:
+            raise RuntimeError("Shape mismatch between dispalcements and forces.")
+        if self._displacements.ndim != 3 or self._displacements.shape[1:] != (
+            len(self._supercell),
+            3,
+        ):
+            raise RuntimeError(
+                "Inconsistent array shape of displacements "
+                f"{self._displacements.shape} with respect to supercell "
+                f"{len(self._supercell)}."
+            )
+        if self._forces.ndim != 3 or self._forces.shape[1:] != (
+            len(self._supercell),
+            3,
+        ):
+            raise RuntimeError(
+                "Inconsistent array shape of forces "
+                f"{self._forces.shape} with respect to supercell "
+                f"{len(self._supercell)}."
+            )
diff --git a/src/symfc/basis_sets/__init__.py b/src/symfc/basis_sets/__init__.py
new file mode 100644
index 0000000..20f1c38
--- /dev/null
+++ b/src/symfc/basis_sets/__init__.py
@@ -0,0 +1,5 @@
+"""Symmetry adapted basis sets of force constants."""
+
+from .basis_sets_base import FCBasisSetBase  # noqa 401
+from .basis_sets_O2 import FCBasisSetO2  # noqa 401
+from .basis_sets_O3 import FCBasisSetO3  # noqa 401
diff --git a/src/symfc/basis_sets/basis_sets_O1.py b/src/symfc/basis_sets/basis_sets_O1.py
new file mode 100644
index 0000000..83633f8
--- /dev/null
+++ b/src/symfc/basis_sets/basis_sets_O1.py
@@ -0,0 +1,140 @@
+"""Symmetry adapted basis sets of 1st order force constants."""
+
+from __future__ import annotations
+
+from typing import Optional
+
+import numpy as np
+from scipy.sparse import csr_array
+
+from symfc.spg_reps import SpgRepsO1
+from symfc.utils.eig_tools import eigsh_projector
+from symfc.utils.matrix_tools_O1 import compressed_projector_sum_rules
+from symfc.utils.utils import SymfcAtoms
+from symfc.utils.utils_O1 import (
+    get_compr_coset_reps_sum,
+    get_lat_trans_compr_matrix,
+    get_lat_trans_decompr_indices,
+)
+
+from .basis_sets_base import FCBasisSetBase
+
+
+class FCBasisSetO1Base(FCBasisSetBase):
+    """Base class of FCBasisSetO1."""
+
+    def __init__(
+        self,
+        supercell: SymfcAtoms,
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+        use_mkl : bool
+            Use MKL or not. Default is False.
+        log_level : int, optional
+            Log level. Default is 0.
+
+        """
+        super().__init__(supercell, use_mkl=use_mkl, log_level=log_level)
+        self._spg_reps = SpgRepsO1(supercell)
+
+    def _get_c_trans(self) -> csr_array:
+        trans_perms = self._spg_reps.translation_permutations
+        n_lp, N = trans_perms.shape
+        decompr_idx = get_lat_trans_decompr_indices(trans_perms)
+        c_trans = get_lat_trans_compr_matrix(decompr_idx, N, n_lp)
+        return c_trans
+
+
+class FCBasisSetO1(FCBasisSetO1Base):
+    """Dense symmetry adapted basis set for 1st order force constants.
+
+    Attributes
+    ----------
+    basis_set : ndarray
+        Compressed force constants basis set. The first dimension n_x (< n_a) is
+        given as a result of compression, which is depends on the system.
+        shape=(n_x, n_bases), dtype='double'
+    full_basis_set : ndarray
+        Full (decompressed) force constants basis set. shape=(N * 3,
+        n_bases), dtype='double'
+    translation_permutations : ndarray
+        Atom indices after lattice translations. shape=(lattice_translations,
+        supercell_atoms), dtype=int.
+
+    """
+
+    def __init__(
+        self,
+        supercell: SymfcAtoms,
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+        log_level : int, optional
+            Log level. Default is 0.
+
+        """
+        super().__init__(supercell, use_mkl=use_mkl, log_level=log_level)
+        self._spg_reps = SpgRepsO1(supercell)
+        self._n_a_compression_matrix: Optional[csr_array] = None
+
+    @property
+    def basis_set(self) -> Optional[np.ndarray]:
+        """Return compressed basis set.
+
+        n_c = len(compressed_indices).
+
+        shape=(n_c, n_bases), dtype='double'.
+
+        """
+        return self._basis_set
+
+    @property
+    def full_basis_set(self) -> Optional[np.ndarray]:
+        """Return full (decompressed) basis set.
+
+        shape=(N*3, n_bases), dtype='double'.
+
+        Data in first dimension is ordered by (N,3).
+
+        """
+        return self._full_basis_set
+
+    @property
+    def compact_compression_matrix(self) -> Optional[np.ndarray]:
+        """Return compression matrix for compact basis set."""
+        pass
+
+    @property
+    def compression_matrix(self) -> Optional[np.ndarray]:
+        """Return compression matrix."""
+        pass
+
+    def run(self) -> FCBasisSetO1:
+        """Compute compressed force constants basis set."""
+        c_trans = self._get_c_trans()
+        coset_reps_sum = get_compr_coset_reps_sum(self._spg_reps)
+        proj_rt = coset_reps_sum
+
+        if len(proj_rt.data) == 0:
+            raise ValueError("No basis vectors exist.")
+
+        c_rt = eigsh_projector(proj_rt)
+        compress_mat = c_trans @ c_rt
+        proj = compressed_projector_sum_rules(compress_mat, self._natom)
+        self._basis_set = eigsh_projector(proj)
+        self._full_basis_set = compress_mat @ self._basis_set
+
+        return self
diff --git a/src/symfc/basis_sets/basis_sets_O2.py b/src/symfc/basis_sets/basis_sets_O2.py
new file mode 100644
index 0000000..a7e9a3e
--- /dev/null
+++ b/src/symfc/basis_sets/basis_sets_O2.py
@@ -0,0 +1,171 @@
+"""Symmetry adapted basis sets of 2nd order force constants."""
+
+from __future__ import annotations
+
+from typing import Optional
+
+import numpy as np
+from scipy.sparse import csr_array
+
+from symfc.spg_reps import SpgRepsO2
+from symfc.utils.eig_tools import (
+    dot_product_sparse,
+    eigsh_projector,
+    eigsh_projector_sumrule,
+)
+from symfc.utils.matrix_tools_O2 import compressed_projector_sum_rules
+from symfc.utils.utils import SymfcAtoms
+from symfc.utils.utils_O2 import (
+    get_compr_coset_reps_sum,
+    get_lat_trans_compr_matrix,
+    get_lat_trans_decompr_indices,
+    get_perm_compr_matrix,
+)
+
+from .basis_sets_base import FCBasisSetBase
+
+
+class FCBasisSetO2(FCBasisSetBase):
+    """Dense symmetry adapted basis set for 2nd order force constants.
+
+    Attributes
+    ----------
+    basis_set : ndarray
+        Compressed force constants basis set. The first dimension n_x (< n_a) is
+        given as a result of compression, which is depends on the system.
+        shape=(n_x * N * 9, n_bases), dtype='double'
+    full_basis_set : ndarray
+        Full (decompressed) force constants basis set. shape=(N * N * 9,
+        n_bases), dtype='double'
+    translation_permutations : ndarray
+        Atom indices after lattice translations. shape=(lattice_translations,
+        supercell_atoms), dtype=int.
+
+    """
+
+    def __init__(
+        self,
+        supercell: SymfcAtoms,
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+        use_mkl : bool
+            Use MKL or not. Default is False.
+        log_level : int, optional
+            Log level. Default is 0.
+
+        """
+        super().__init__(supercell, use_mkl=use_mkl, log_level=log_level)
+        self._spg_reps = SpgRepsO2(supercell)
+        self._n_a_compression_matrix: Optional[csr_array] = None
+
+    @property
+    def basis_set(self) -> Optional[np.ndarray]:
+        """Return compressed basis set.
+
+        n_c = len(compressed_indices).
+
+        shape=(n_c*N*3*3, n_bases), dtype='double'.
+
+        Data in first dimension is ordered by (n_c,N,3,3).
+
+        """
+        return self._basis_set
+
+    @property
+    def compression_matrix(self) -> Optional[csr_array]:
+        """Return compression matrix.
+
+        This expands fc basis_sets to (N*N*3*3, n_bases).
+
+        """
+        c_trans = self._get_c_trans()
+        return dot_product_sparse(
+            c_trans, self._n_a_compression_matrix, use_mkl=self._use_mkl
+        )
+
+    @property
+    def compact_compression_matrix(self) -> Optional[csr_array]:
+        """Return compact compression matrix.
+
+        This expands basis_sets to (n_a*N*3*3, n_bases).
+
+        """
+        n_lp = self.translation_permutations.shape[0]
+        return self._n_a_compression_matrix / np.sqrt(n_lp)
+
+    @property
+    def decompression_indices(self) -> np.ndarray:
+        """Return decompression indices in (N,N,3,3) order."""
+        trans_perms = self.translation_permutations
+        return get_lat_trans_decompr_indices(trans_perms)
+
+    def run(self) -> FCBasisSetO2:
+        """Compute compressed force constants basis set.
+
+        compression using C(pt)
+            = eigvecs of C(trans).T @ C(perm) @ C(perm).T @ C(trans)
+        proj_rpt = c_pt.transpose() @ coset_reps_sum @ c_pt
+
+        [C_pt @ C_trans.T] @ P_sum @ [C_trans @ C_pt] @ C_rpt
+         = C_rpt @ [C_rpt.T @ C_pt @ C_trans.T] @ P_sum @ [C_trans @ C_pt @ C_rpt]
+         = C_rpt @ compression_mat.T @ (I - P_sum^(c)) @ compression_mat
+         = C_rpt @ compression_mat.T @ (I - Csum(c) @ Csum(c).T) @ compression_mat
+         = C_rpt @ proj
+
+        compress_mat = c_trans @ c_pt @ c_rpt
+
+        """
+        N = self._natom
+        c_trans = self._get_c_trans()
+        n_a_compress_mat = self._get_n_a_compress_mat(c_trans)
+        compress_mat = dot_product_sparse(
+            c_trans, n_a_compress_mat, use_mkl=self._use_mkl
+        )
+        proj = compressed_projector_sum_rules(compress_mat, N, use_mkl=self._use_mkl)
+        eigvecs = eigsh_projector_sumrule(proj)
+
+        if self._log_level:
+            print(f"Final size of basis set: {eigvecs.shape}")
+
+        self._basis_set = eigvecs
+        self._n_a_compression_matrix = n_a_compress_mat
+        return self
+
+    def _get_n_a_compress_mat(self, c_trans: csr_array) -> csr_array:
+        """Return compression matrix without c_trans mutiplied.
+
+        This compression matrix is preserved as a class instance variable.
+        The full compression matrix is obtained by
+
+        c_trans @ n_a_compression_matrix.
+
+        The compact compression matrix is obtained by
+
+        n_a_compression_matrix / sqrt(n_lp).
+
+        """
+        N = self._natom
+        c_perm = get_perm_compr_matrix(N)
+        c_pt = dot_product_sparse(c_perm.T, c_trans, use_mkl=self._use_mkl)
+        proj_pt = dot_product_sparse(c_pt.T, c_pt, use_mkl=self._use_mkl)
+        coset_reps_sum = get_compr_coset_reps_sum(self._spg_reps)
+        c_pt = eigsh_projector(proj_pt)
+        proj_rpt = dot_product_sparse(coset_reps_sum, c_pt, use_mkl=self._use_mkl)
+        proj_rpt = dot_product_sparse(c_pt.T, proj_rpt)
+        c_rpt = eigsh_projector(proj_rpt)
+        n_a_compress_mat = dot_product_sparse(c_pt, c_rpt, use_mkl=self._use_mkl)
+        return n_a_compress_mat
+
+    def _get_c_trans(self) -> csr_array:
+        trans_perms = self._spg_reps.translation_permutations
+        n_lp, N = trans_perms.shape
+        decompr_idx = get_lat_trans_decompr_indices(trans_perms)
+        c_trans = get_lat_trans_compr_matrix(decompr_idx, N, n_lp)
+        return c_trans
diff --git a/src/symfc/basis_sets/basis_sets_O3.py b/src/symfc/basis_sets/basis_sets_O3.py
new file mode 100644
index 0000000..4b50b12
--- /dev/null
+++ b/src/symfc/basis_sets/basis_sets_O3.py
@@ -0,0 +1,218 @@
+"""Symmetry adapted basis sets of 3rd order force constants."""
+
+from __future__ import annotations
+
+import time
+from typing import Optional, Union
+
+import numpy as np
+from scipy.sparse import coo_array, csr_array
+
+from symfc.spg_reps import SpgRepsO3
+from symfc.utils.eig_tools import (
+    dot_product_sparse,
+    eigsh_projector,
+    eigsh_projector_sumrule,
+)
+from symfc.utils.matrix_tools_O3 import (
+    compressed_projector_sum_rules,
+    get_perm_compr_matrix_O3,
+)
+from symfc.utils.utils import SymfcAtoms
+from symfc.utils.utils_O3 import (
+    get_compr_coset_reps_sum_O3,
+    get_lat_trans_compr_matrix_O3,
+)
+
+from . import FCBasisSetBase
+
+
+def print_sp_matrix_size(c: Union[csr_array, coo_array], header: str):
+    """Show sparse matrix size."""
+    print(header, c.shape, len(c.data))
+
+
+class FCBasisSetO3(FCBasisSetBase):
+    """Symmetry adapted basis set for 3rd order force constants.
+
+    Attributes
+    ----------
+    basis_set : ndarray
+        Compressed force constants basis set.
+        shape=(n_a * N * 9, n_bases), dtype='double'
+    full_basis_set : ndarray
+        Full (decompressed) force constants basis set.
+        shape=(N * N * 9, n_bases), dtype='double'
+    decompression_indices : ndarray
+        Decompression indices in (N,N,3,3) order.
+        shape=(N^2*9,), dtype='int_'.
+    compresssion_indices : ndarray
+        Compression indices in (n_a,N,3,3) order.
+        shape=(n_a*N*9, n_lp), dtype='int_'.
+    translation_permutations : ndarray
+        Atom indices after lattice translations.
+        shape=(lattice_translations, supercell_atoms), dtype=int.
+
+    """
+
+    def __init__(
+        self,
+        supercell: SymfcAtoms,
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+        log_level : int, optional
+            Log level. Default is 0.
+
+        """
+        super().__init__(supercell, use_mkl=use_mkl, log_level=log_level)
+        self._spg_reps = SpgRepsO3(supercell)
+
+    @property
+    def basis_set(self) -> Optional[csr_array]:
+        """Return compressed basis set.
+
+        n_c = len(compressed_indices).
+
+        shape=(n_c*N*N*3*3*3, n_bases), dtype='double'.
+
+        Data in first dimension is ordered by (n_c,N,N,3,3,3).
+
+        """
+        return self._basis_set
+
+    @property
+    def compression_matrix(self) -> Optional[csr_array]:
+        """Return compression matrix.
+
+        This expands fc basis_sets to (N*N*N*3*3*3, n_bases).
+
+        """
+        trans_perms = self._spg_reps.translation_permutations
+        c_trans = get_lat_trans_compr_matrix_O3(trans_perms)
+        return dot_product_sparse(
+            c_trans, self._n_a_compression_matrix, use_mkl=self._use_mkl
+        )
+
+    @property
+    def compact_compression_matrix(self) -> Optional[csr_array]:
+        """Return compact compression matrix.
+
+        This expands basis_sets to (n_a*N*N*3*3*3, n_bases).
+
+        """
+        n_lp = self.translation_permutations.shape[0]
+        return self._n_a_compression_matrix / np.sqrt(n_lp)
+
+    def run(self) -> FCBasisSetO3:
+        """Compute compressed force constants basis set."""
+        trans_perms = self._spg_reps.translation_permutations
+        N = self._natom
+
+        tt0 = time.time()
+        c_trans = get_lat_trans_compr_matrix_O3(trans_perms)
+        tt1 = time.time()
+
+        coset_reps_sum = get_compr_coset_reps_sum_O3(self._spg_reps)
+        print_sp_matrix_size(coset_reps_sum, " R_(coset):")
+        tt2 = time.time()
+
+        c_pt = self._get_perm_trans_compr_matrix(c_trans, N)
+
+        tt4 = time.time()
+
+        n_a_compress_mat = self._get_n_a_compress_mat(c_pt, coset_reps_sum)
+        compress_mat = dot_product_sparse(
+            c_trans, n_a_compress_mat, use_mkl=self._use_mkl
+        )
+        print_sp_matrix_size(compress_mat, " compression matrix:")
+        tt5 = time.time()
+
+        proj = compressed_projector_sum_rules(compress_mat, N, use_mkl=self._use_mkl)
+        print_sp_matrix_size(proj, " P_(perm,trans,coset,sum):")
+        tt6 = time.time()
+
+        eigvecs = eigsh_projector_sumrule(proj)
+        print(" basis (size) =", eigvecs.shape)
+
+        tt7 = time.time()
+
+        print("  t (init., trans)        = ", tt1 - tt0)
+        print("  t (coset_reps_sum)      = ", tt2 - tt1)
+        print("  t (dot, trans, perm)    = ", tt4 - tt2)
+        print("  t (rot, trans, perm)    = ", tt5 - tt4)
+        print("  t (proj_st)             = ", tt6 - tt5)
+        print("  t (eigh(svd))           = ", tt7 - tt6)
+        # print('  t (reconstruction)      = ', tt8-tt7)
+
+        self._basis_set = eigvecs
+        self._n_a_compression_matrix = n_a_compress_mat
+        # self._compression_matrix = compress_mat
+
+        return self
+
+    def _get_perm_trans_compr_matrix(self, c_trans: csr_array, N: int):
+        """Return perm trans compression matrix.
+
+        compression using C(pt)
+            = eigvecs of C(trans).T @ C(perm) @ C(perm).T @ C(trans)
+        proj_rpt = c_pt.T @ coset_reps_sum @ c_pt
+
+        C(pt) = C(perm).T @ C(trans)
+
+        """
+        c_perm = get_perm_compr_matrix_O3(N)
+        print_sp_matrix_size(c_perm, " C_perm:")
+        c_pt = dot_product_sparse(c_perm.T, c_trans, use_mkl=self._use_mkl)
+        print_sp_matrix_size(c_pt, " C_(perm,trans):")
+        proj_pt = dot_product_sparse(c_pt.T, c_pt, use_mkl=self._use_mkl)
+        print_sp_matrix_size(proj_pt, " P_(perm,trans):")
+        c_pt = eigsh_projector(proj_pt)
+        print_sp_matrix_size(c_pt, " C_(perm,trans,compressed):")
+        return c_pt
+
+    def _get_n_a_compress_mat(
+        self, c_pt: csr_array, coset_reps_sum: csr_array
+    ) -> csr_array:
+        """Return compression matrix without c_trans mutiplied.
+
+        This compression matrix is preserved as a class instance variable.
+        The full compression matrix is obtained by
+
+        c_trans @ n_a_compression_matrix.
+
+        The compact compression matrix is obtained by
+
+        n_a_compression_matrix / sqrt(n_lp).
+
+        About compression matrix
+        ------------------------
+        compress_mat = c_trans @ c_pt @ c_rpt
+
+        c_trans : compression matrix by lattice translations
+        c_pt : basis vectors of projection matrix of perm & trans
+        c_rpt : basis vectors of projection matrix of  coset rots & perm & trans
+
+        More precisely,
+
+        [C_pt @ C_trans.T] @ P_sum @ [C_trans @ C_pt] @ C_rpt
+         = C_rpt @ [C_rpt.T @ C_pt @ C_trans.T] @ P_sum @ [C_trans @ C_pt @ C_rpt]
+         = C_rpt @ compression_mat.T @ (I - P_sum^(c)) @ compression_mat
+         = C_rpt @ compression_mat.T @ (I - Csum(c) @ Csum(c).T) @ compression_mat
+         = C_rpt @ proj
+
+        """
+        proj_rpt = dot_product_sparse(coset_reps_sum, c_pt, use_mkl=self._use_mkl)
+        proj_rpt = dot_product_sparse(c_pt.T, proj_rpt, use_mkl=self._use_mkl)
+        print_sp_matrix_size(proj_rpt, " P_(perm,trans,coset):")
+        c_rpt = eigsh_projector(proj_rpt)
+        print_sp_matrix_size(c_rpt, " C_(perm,trans,coset):")
+        n_a_compress_mat = dot_product_sparse(c_pt, c_rpt, use_mkl=self._use_mkl)
+        print_sp_matrix_size(n_a_compress_mat, " n_a_compression matrix:")
+        return n_a_compress_mat
diff --git a/src/symfc/basis_sets/basis_sets_base.py b/src/symfc/basis_sets/basis_sets_base.py
new file mode 100644
index 0000000..36d15ef
--- /dev/null
+++ b/src/symfc/basis_sets/basis_sets_base.py
@@ -0,0 +1,65 @@
+"""Symmetry adapted basis sets of force constants."""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from typing import Optional
+
+import numpy as np
+
+from symfc.spg_reps import SpgRepsBase
+from symfc.utils.utils import SymfcAtoms
+
+
+class FCBasisSetBase(ABC):
+    """Abstract base class of symmetry adapted basis set for force constants."""
+
+    def __init__(
+        self,
+        supercell: SymfcAtoms,
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+        log_level : int, optional
+            Log level. Default is 0.
+
+        """
+        self._natom = len(supercell)
+        self._use_mkl = use_mkl
+        self._log_level = log_level
+        self._basis_set: Optional[np.ndarray] = None
+        self._spg_reps: Optional[SpgRepsBase] = None
+
+    @property
+    @abstractmethod
+    def basis_set(self) -> Optional[np.ndarray]:
+        """Return (compressed) basis set."""
+        pass
+
+    @property
+    @abstractmethod
+    def compact_compression_matrix(self) -> Optional[np.ndarray]:
+        """Return compression matrix for compact basis set."""
+        pass
+
+    @property
+    @abstractmethod
+    def compression_matrix(self) -> Optional[np.ndarray]:
+        """Return compression matrix."""
+        pass
+
+    @property
+    def translation_permutations(self) -> np.ndarray:
+        """Return permutations by lattice translation."""
+        return self._spg_reps.translation_permutations
+
+    @abstractmethod
+    def run(self):
+        """Run basis set calculation."""
+        pass
diff --git a/src/symfc/solvers/__init__.py b/src/symfc/solvers/__init__.py
new file mode 100644
index 0000000..a3a2495
--- /dev/null
+++ b/src/symfc/solvers/__init__.py
@@ -0,0 +1,5 @@
+"""Force constants solvers."""
+
+from .solver_base import FCSolverBase  # noqa 401
+from .solver_O2 import FCSolverO2  # noqa 401
+from .solver_O2O3 import FCSolverO2O3  # noqa 401
diff --git a/src/symfc/solvers/solver_O2.py b/src/symfc/solvers/solver_O2.py
new file mode 100644
index 0000000..d73cccb
--- /dev/null
+++ b/src/symfc/solvers/solver_O2.py
@@ -0,0 +1,254 @@
+"""2nd order force constants solver."""
+
+from __future__ import annotations
+
+import time
+from typing import Optional
+
+import numpy as np
+from scipy.sparse import csr_array
+
+from symfc.basis_sets import FCBasisSetO2
+from symfc.utils.eig_tools import dot_product_sparse
+from symfc.utils.solver_funcs import fit, get_batch_slice, solve_linear_equation
+
+from .solver_base import FCSolverBase
+
+
+class FCSolverO2(FCSolverBase):
+    """Third order force constants solver."""
+
+    def __init__(
+        self,
+        basis_set: FCBasisSetO2,
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method."""
+        super().__init__(basis_set, use_mkl=use_mkl, log_level=log_level)
+
+    @property
+    def full_fc(self) -> Optional[np.ndarray]:
+        """Return full force constants.
+
+        Returns
+        -------
+        np.ndarray
+            shape=(N, N, 3, 3), dtype='double', order='C'
+
+        """
+        N = self._natom
+        fc = self._basis_set.basis_set @ self._coefs
+        return np.array(
+            (self._basis_set.compression_matrix @ fc).reshape((-1, N, 3, 3)),
+            dtype="double",
+            order="C",
+        )
+
+    @property
+    def compact_fc(self) -> Optional[np.ndarray]:
+        """Return full force constants.
+
+        Returns
+        -------
+        np.ndarray
+            shape=(n_a, N, 3, 3), dtype='double', order='C'
+
+        """
+        N = self._natom
+        fc = self._basis_set.basis_set @ self._coefs
+        return np.array(
+            (self._basis_set.compact_compression_matrix @ fc).reshape((-1, N, 3, 3)),
+            dtype="double",
+            order="C",
+        )
+
+    def solve(
+        self,
+        displacements: np.ndarray,
+        forces: np.ndarray,
+    ) -> FCSolverO2:
+        """Solve coefficients of basis set from displacements and forces.
+
+        Parameters
+        ----------
+        displacements : ndarray
+            Displacements of atoms in Cartesian coordinates.
+            shape=(n_snapshot, N, 3), dtype='double'
+        forces : ndarray
+            Forces of atoms in Cartesian coordinates.
+            shape=(n_snapshot, N, 3), dtype='double'
+
+        Returns
+        -------
+        self : FCSolverO2
+
+        """
+        n_data = forces.shape[0]
+        f = forces.reshape(n_data, -1)
+        d = displacements.reshape(n_data, -1)
+        self._coefs = run_solver_sparse_O2(
+            d, f, self._basis_set.compression_matrix, self._basis_set.basis_set
+        )
+        return self
+
+
+def get_training_from_full_basis(disps, forces, full_basis):
+    """Return training dataset (X, y).
+
+    The dataset is transformed from displacements, forces, and full basis-set.
+
+    disps: (n_samples, N3)
+    forces: (n_samples, N3)
+    full_basis: (n_basis, N, N, 3, 3)
+
+    X: features (calculated from displacements),(n_samples * N3, N_basis)
+    y: observations (forces), (n_samples * N3)
+
+    Note
+    ----
+    #algorithm 1
+    X = []
+    for disp_st, force_st in zip(disps, forces):
+        disp_reshape = disp_st.reshape(-1)
+        X_tmp = []
+        for basis in full_basis:
+            Bb = basis.transpose((0,2,1,3)).reshape(N3,N3)
+            X_tmp.append(-np.dot(Bb, disp_reshape))
+        X_tmp = np.array(X_tmp).T
+        X.append(X_tmp)
+    X = np.vstack(X)
+
+    #algorithm 2
+    disps_reshape = disps.reshape((n_samples, N3))
+    X = np.zeros((N3*n_samples, n_basis))
+    for i, basis in enumerate(full_basis):
+        Bb = basis.transpose((0,2,1,3)).reshape(N3,N3)
+        prod = - np.dot(disps_reshape, Bb)
+        X[:,i] = np.ravel(prod)
+
+    """
+    N3 = full_basis.shape[1] * 3
+    n_basis = full_basis.shape[0]
+    n_samples = disps.shape[0]
+
+    disps_reshape = disps.reshape((n_samples, N3))
+    full_basis = full_basis.transpose((1, 3, 2, 4, 0)).reshape((N3, -1))
+    X = -np.dot(disps_reshape, full_basis)
+    X = X.reshape((-1, n_basis))
+    y = forces.reshape(-1)
+
+    return X, y
+
+
+def run_solver_dense_O2(disps, forces, compress_mat, compress_eigvecs):
+    """Estimating coeffs. in X @ coeffs = y using least squares.
+
+    X: features (calculated from displacements), (n_samples * N3, N_basis)
+    y: observations (forces), (n_samples * N3)
+
+    """
+    N3 = disps.shape[1]
+    N = N3 // 3
+
+    full_basis = compress_mat @ compress_eigvecs
+    n_basis = full_basis.shape[1]
+    coefs = _solve(disps, forces, full_basis.T.reshape((n_basis, N, N, 3, 3)))
+    return coefs
+
+
+def run_solver_sparse_O2(disps, forces, compress_mat, compress_eigvecs, batch_size=200):
+    """Estimating coeffs. in X @ coeffs = y by solving normal equation.
+
+    (X.T @ X) @ coeffs = X.T @ y
+
+    X = displacements @ compress_mat @ compress_eigvecs
+    Matrix reshapings are appropriately applied.
+    X: features (n_samples * N3, N_basis)
+    y: observations (forces), (n_samples * N3)
+
+    """
+    XTX, XTy = _get_training(
+        disps, forces, compress_mat, compress_eigvecs, batch_size=batch_size
+    )
+    coefs = solve_linear_equation(XTX, XTy)
+    return coefs
+
+
+def _get_training(
+    disps, forces, compress_mat, compress_eigvecs, batch_size=200, use_mkl=False
+):
+    r"""Calculate X.T @ X and X.T @ y.
+
+    X = displacements @ compress_mat @ compress_eigvecs
+
+    displacements: (n_samples, N33)
+    compress_mat: (NN33, n_compr)
+    compress_eigvecs: (n_compr, n_basis)
+    Matrix reshapings are appropriately applied to compress_mat
+    and its products.
+
+    X.T @ X and X.T @ y are sequentially calculated using divided dataset.
+        X.T @ X = \sum_i X_i.T @ X_i
+        X.T @ y = \sum_i X_i.T @ y_i (i: batch index)
+
+    """
+    N3 = disps.shape[1]
+    N = N3 // 3
+    n_basis = compress_eigvecs.shape[1]
+    n_compr = compress_mat.shape[1]
+
+    t1 = time.time()
+    compress_mat = _csr_NN33_to_N3N3(compress_mat, N).reshape((N3, -1)).tocsr()
+    t2 = time.time()
+    print(" reshape(compr):   ", t2 - t1)
+
+    begin_batch, end_batch = get_batch_slice(disps.shape[0], batch_size)
+    XTX = np.zeros((n_basis, n_basis), dtype=float)
+    XTy = np.zeros(n_basis, dtype=float)
+    for begin, end in zip(begin_batch, end_batch):
+        t01 = time.time()
+        disps_batch = csr_array(disps[begin:end])
+        X = dot_product_sparse(
+            disps_batch, compress_mat, use_mkl=use_mkl, dense=True
+        ).reshape((-1, n_compr))
+        X = X @ compress_eigvecs
+        y = -forces[begin:end].reshape(-1)
+        XTX += X.T @ X
+        XTy += X.T @ y
+        t02 = time.time()
+        print(" solver_block:", end, ":, t =", t02 - t01)
+    t3 = time.time()
+    print(" (disp @ compr @ eigvecs).T @ (disp @ compr @ eigvecs):", t3 - t2)
+    return XTX, XTy
+
+
+def _solve(disps, forces, full_basis):
+    """Solve force constants."""
+    X, y = get_training_from_full_basis(disps, forces, full_basis)
+    coefs = fit(X, y)
+    return coefs
+
+
+def _NN33_to_N3N3(row, N):
+    """Reorder row indices in a sparse matrix (NN33->N3N3)."""
+    i, rem = np.divmod(row, 9 * N)
+    j, rem = np.divmod(rem, 9)
+    a, b = np.divmod(rem, 3)
+
+    vec = i * 9 * N
+    vec += j * 3
+    vec += a * 3 * N + b
+    return vec
+
+
+def _csr_NN33_to_N3N3(mat, N):
+    """Reorder row indices in a sparse matrix (NN33->N3N3).
+
+    Return reordered csr_array.
+    """
+    NN33, nx = mat.shape
+    mat = mat.tocoo()
+    row = _NN33_to_N3N3(mat.row, N)
+    mat = csr_array((mat.data, (row, mat.col)), shape=(NN33, nx))
+    return mat
diff --git a/src/symfc/solvers/solver_O2O3.py b/src/symfc/solvers/solver_O2O3.py
new file mode 100644
index 0000000..593c9e1
--- /dev/null
+++ b/src/symfc/solvers/solver_O2O3.py
@@ -0,0 +1,583 @@
+"""Solver of 2nd and 3rd order force constants simultaneously."""
+
+from __future__ import annotations
+
+import time
+from collections.abc import Sequence
+from typing import Literal, Optional
+
+import numpy as np
+from scipy.sparse import csr_array
+
+from symfc.basis_sets import FCBasisSetO2, FCBasisSetO3
+from symfc.utils.eig_tools import dot_product_sparse
+from symfc.utils.solver_funcs import get_batch_slice, solve_linear_equation
+from symfc.utils.utils_O2 import get_perm_compr_matrix
+
+from .solver_base import FCSolverBase
+from .solver_O2 import get_training_from_full_basis
+from .solver_O3 import csr_NNN333_to_NN33N3, set_2nd_disps
+
+
+class FCSolverO2O3(FCSolverBase):
+    """Simultaneous second and third order force constants solver."""
+
+    def __init__(
+        self,
+        basis_set: Sequence[FCBasisSetO2, FCBasisSetO3],
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method."""
+        super().__init__(basis_set, use_mkl=use_mkl, log_level=log_level)
+
+    def solve(
+        self, displacements: np.ndarray, forces: np.ndarray, batch_size: int = 200
+    ) -> FCSolverO2O3:
+        """Solve force constants.
+
+        Note
+        ----
+        self._coefs = (coefs_fc2, coefs_fc3)
+
+        Parameters
+        ----------
+        displacements : ndarray
+            Displacements of atoms in Cartesian coordinates. shape=(n_snapshot,
+            N, 3), dtype='double'
+        forces : ndarray
+            Forces of atoms in Cartesian coordinates. shape=(n_snapshot, N, 3),
+            dtype='double'
+
+        Returns
+        -------
+        ndarray
+            Force constants. shape=(n_a, N, 3, 3) or (N, N, 3, 3). See
+            `is_compact_fc` parameter. dtype='double', order='C'
+
+        """
+        n_data = forces.shape[0]
+        f = forces.reshape(n_data, -1)
+        d = displacements.reshape(n_data, -1)
+
+        fc2_basis: FCBasisSetO2 = self._basis_set[0]
+        fc3_basis: FCBasisSetO3 = self._basis_set[1]
+        compress_mat_fc2 = fc2_basis.compression_matrix
+        basis_set_fc2 = fc2_basis.basis_set
+        compress_mat_fc3 = fc3_basis.compression_matrix
+        basis_set_fc3 = fc3_basis.basis_set
+
+        self._coefs = run_solver_sparse_O2O3(
+            d,
+            f,
+            compress_mat_fc2,
+            compress_mat_fc3,
+            basis_set_fc2,
+            basis_set_fc3,
+            batch_size=batch_size,
+            use_mkl=self._use_mkl,
+        )
+        return self
+
+    @property
+    def full_fc(self) -> Optional[tuple[np.ndarray, np.ndarray]]:
+        """Return full force constants.
+
+        Returns
+        -------
+        tuple[np.ndarray, np.ndarray]
+            shape=(N, N, 3, 3), dtype='double', order='C'
+            shape=(N, N, N, 3, 3, 3), dtype='double', order='C'
+
+        """
+        return self._recover_fcs("full")
+
+    @property
+    def compact_fc(self) -> Optional[tuple[np.ndarray, np.ndarray]]:
+        """Return full force constants.
+
+        Returns
+        -------
+        tuple[np.ndarray, np.ndarray]
+            shape=(n_a, N, 3, 3), dtype='double', order='C'
+            shape=(n_a, N, N, 3, 3, 3), dtype='double', order='C'
+
+        """
+        return self._recover_fcs("compact")
+
+    def _recover_fcs(
+        self, comp_mat_type: str = Literal["full", "compact"]
+    ) -> Optional[tuple[np.ndarray, np.ndarray]]:
+        if self._coefs is None:
+            return None
+
+        fc2_basis: FCBasisSetO2 = self._basis_set[0]
+        fc3_basis: FCBasisSetO3 = self._basis_set[1]
+        if comp_mat_type == "full":
+            comp_mat_fc2 = fc2_basis.compression_matrix
+            comp_mat_fc3 = fc3_basis.compression_matrix
+        elif comp_mat_type == "compact":
+            comp_mat_fc2 = fc2_basis.compact_compression_matrix
+            comp_mat_fc3 = fc3_basis.compact_compression_matrix
+        else:
+            raise ValueError("Invalid comp_mat_type.")
+
+        N = self._natom
+        fc2 = fc2_basis.basis_set @ self._coefs[0]
+        fc2 = np.array(
+            (comp_mat_fc2 @ fc2).reshape((-1, N, 3, 3)), dtype="double", order="C"
+        )
+        fc3 = fc3_basis.basis_set @ self._coefs[1]
+        fc3 = np.array(
+            (comp_mat_fc3 @ fc3).reshape((-1, N, N, 3, 3, 3)),
+            dtype="double",
+            order="C",
+        )
+        return fc2, fc3
+
+
+def get_training_exact(
+    disps,
+    forces,
+    compress_mat_fc2,
+    compress_mat_fc3,
+    compress_eigvecs_fc2,
+    compress_eigvecs_fc3,
+    batch_size=200,
+    use_mkl=False,
+):
+    r"""Calculate X.T @ X and X.T @ y.
+
+    X = displacements @ compress_mat @ compress_eigvecs
+    X = np.hstack([X_fc2, X_fc3])
+
+    displacements (fc2): (n_samples, N3)
+    displacements (fc3): (n_samples, NN33)
+    compress_mat_fc2: (NN33, n_compr)
+    compress_mat_fc3: (NNN333, n_compr_fc3)
+    compress_eigvecs_fc2: (n_compr_fc2, n_basis_fc2)
+    compress_eigvecs_fc3: (n_compr_fc3, n_basis_fc3)
+    Matrix reshapings are appropriately applied to compress_mat
+    and its products.
+
+    X.T @ X and X.T @ y are sequentially calculated using divided dataset.
+        X.T @ X = \sum_i X_i.T @ X_i
+        X.T @ y = \sum_i X_i.T @ y_i (i: batch index)
+
+    """
+    N3 = disps.shape[1]
+    N = N3 // 3
+    NN33 = N3 * N3
+    n_basis_fc2 = compress_eigvecs_fc2.shape[1]
+    n_basis_fc3 = compress_eigvecs_fc3.shape[1]
+    n_compr_fc3 = compress_mat_fc3.shape[1]
+    n_basis = n_basis_fc2 + n_basis_fc3
+
+    t1 = time.time()
+    full_basis_fc2 = compress_mat_fc2 @ compress_eigvecs_fc2
+    X2, y_all = get_training_from_full_basis(
+        disps, forces, full_basis_fc2.T.reshape((n_basis_fc2, N, N, 3, 3))
+    )
+    t2 = time.time()
+    print(" training data (fc2):    ", t2 - t1)
+
+    t1 = time.time()
+    c_perm_fc2 = get_perm_compr_matrix(N)
+    compress_mat_fc3 = (
+        csr_NNN333_to_NN33N3(compress_mat_fc3, N).reshape((NN33, -1)).tocsr()
+    )
+    compress_mat_fc3 = -0.5 * (c_perm_fc2.T @ compress_mat_fc3)
+    t2 = time.time()
+    print(" precond. compress_mat (for fc3):", t2 - t1)
+
+    sparse_disps = True if use_mkl else False
+    mat33 = np.zeros((n_compr_fc3, n_compr_fc3), dtype=float)
+    mat23 = np.zeros((n_basis_fc2, n_compr_fc3), dtype=float)
+    mat3y = np.zeros(n_compr_fc3, dtype=float)
+    begin_batch, end_batch = get_batch_slice(disps.shape[0], batch_size)
+    for begin, end in zip(begin_batch, end_batch):
+        t01 = time.time()
+        disps_batch = set_2nd_disps(disps[begin:end], sparse=sparse_disps)
+        disps_batch = disps_batch @ c_perm_fc2
+        X3 = dot_product_sparse(
+            disps_batch, compress_mat_fc3, use_mkl=use_mkl, dense=True
+        ).reshape((-1, n_compr_fc3))
+        y_batch = forces[begin:end].reshape(-1)
+        mat23 += X2[begin * N3 : end * N3].T @ X3
+        mat33 += X3.T @ X3
+        mat3y += X3.T @ y_batch
+        t02 = time.time()
+        print(" solver_block:", end, ":, t =", t02 - t01)
+
+    XTX = np.zeros((n_basis, n_basis), dtype=float)
+    XTy = np.zeros(n_basis, dtype=float)
+    XTX[:n_basis_fc2, :n_basis_fc2] = X2.T @ X2
+    XTX[:n_basis_fc2, n_basis_fc2:] = mat23 @ compress_eigvecs_fc3
+    XTX[n_basis_fc2:, :n_basis_fc2] = XTX[:n_basis_fc2, n_basis_fc2:].T
+    XTX[n_basis_fc2:, n_basis_fc2:] = (
+        compress_eigvecs_fc3.T @ mat33 @ compress_eigvecs_fc3
+    )
+    XTy[:n_basis_fc2] = X2.T @ y_all
+    XTy[n_basis_fc2:] = compress_eigvecs_fc3.T @ mat3y
+
+    t3 = time.time()
+    print(" (disp @ compr @ eigvecs).T @ (disp @ compr @ eigvecs):", t3 - t2)
+    return XTX, XTy
+
+
+def run_solver_sparse_O2O3(
+    disps,
+    forces,
+    compress_mat_fc2,
+    compress_mat_fc3,
+    compress_eigvecs_fc2,
+    compress_eigvecs_fc3,
+    batch_size=200,
+    use_mkl=False,
+):
+    """Estimate coeffs. in X @ coeffs = y.
+
+    X_fc2 = displacements_fc2 @ compress_mat_fc2 @ compress_eigvecs_fc2
+    X_fc3 = displacements_fc3 @ compress_mat_fc3 @ compress_eigvecs_fc3
+    X = np.hstack([X_fc2, X_fc3])
+
+    Matrix reshapings are appropriately applied.
+    X: features (n_samples * N3, N_basis_fc2 + N_basis_fc3)
+    y: observations (forces), (n_samples * N3)
+
+    """
+    XTX, XTy = get_training_exact(
+        disps,
+        forces,
+        compress_mat_fc2,
+        compress_mat_fc3,
+        compress_eigvecs_fc2,
+        compress_eigvecs_fc3,
+        batch_size=batch_size,
+        use_mkl=use_mkl,
+    )
+    coefs = solve_linear_equation(XTX, XTy)
+    n_basis_fc2 = compress_eigvecs_fc2.shape[1]
+    coefs_fc2, coefs_fc3 = coefs[:n_basis_fc2], coefs[n_basis_fc2:]
+    return coefs_fc2, coefs_fc3
+
+
+def _NNN333_to_NN33N3_core(row, N):
+    """Reorder row indices in a sparse matrix (NNN333->NN33N3)."""
+    # i
+    div, rem = np.divmod(row, 27 * (N**2))
+    row_update = div * 27 * (N**2)
+    # j
+    div, rem = np.divmod(rem, 27 * N)
+    row_update += div * 27 * N
+    # k
+    div, rem = np.divmod(rem, 27)
+    row_update += div * 3
+    # a
+    div, rem = np.divmod(rem, 9)
+    row_update += div * 9 * N
+    # b, c
+    div, rem = np.divmod(rem, 3)
+    row_update += div * 3 * N + rem
+    return row_update
+
+
+def _NNN333_to_NN33N3(row, N, n_batch=10):
+    """Reorder row indices in a sparse matrix (NNN333->NN33N3).
+
+    This function uses divided row index sets.
+    """
+    batch_size = len(row) // n_batch
+    begin_batch, end_batch = get_batch_slice(len(row), batch_size)
+    for begin, end in zip(begin_batch, end_batch):
+        row[begin:end] = _NNN333_to_NN33N3_core(row[begin:end], N)
+    return row
+
+
+def reshape_compress_mat(mat, N, n_batch=10):
+    """Reorder row indices in a sparse matrix (NNN333->NN33N3).
+
+    This function involves reshaping the sparse matrix
+    into another sparse matrix (NN33N3,Nx) -> (NN33, N3Nx).
+
+    Return reordered csr_matrix.
+    """
+    NNN333, nx = mat.shape
+    mat = mat.tocoo()
+    mat.row = _NNN333_to_NN33N3(mat.row, N, n_batch=n_batch)
+
+    NN33 = (N**2) * 9
+    N3 = N * 3
+    """
+    # reshape: (NN33N3,Nx) -> (NN33, N3Nx)
+    mat.row, rem = np.divmod(mat.row, N3)
+    mat.col += rem * nx
+    """
+    batch_size = len(mat.row) // n_batch
+    begin_batch, end_batch = get_batch_slice(len(mat.row), batch_size)
+    for begin, end in zip(begin_batch, end_batch):
+        mat.row[begin:end], rem = np.divmod(mat.row[begin:end], N3)
+        mat.col[begin:end] += rem * nx
+
+    return csr_array((mat.data, (mat.row, mat.col)), shape=(NN33, N3 * nx))
+
+
+def get_training(
+    disps,
+    forces,
+    compress_mat_fc2,
+    compress_mat_fc3,
+    compress_eigvecs_fc2,
+    compress_eigvecs_fc3,
+    batch_size=100,
+    use_mkl=False,
+    compress_perm_fc2=False,
+):
+    r"""Calculate X.T @ X and X.T @ y.
+
+    X = displacements @ compress_mat @ compress_eigvecs
+    X = np.hstack([X_fc2, X_fc3])
+
+    displacements (fc2): (n_samples, N3)
+    displacements (fc3): (n_samples, NN33)
+    compress_mat_fc2: (NN33, n_compr)
+    compress_mat_fc3: (NNN333, n_compr_fc3)
+    compress_eigvecs_fc2: (n_compr_fc2, n_basis_fc2)
+    compress_eigvecs_fc3: (n_compr_fc3, n_basis_fc3)
+    Matrix reshapings are appropriately applied to compress_mat
+    and its products.
+
+    X.T @ X and X.T @ y are sequentially calculated using divided dataset.
+    X.T @ X = \sum_i X_i.T @ X_i
+    X.T @ y = \sum_i X_i.T @ y_i (i: batch index)
+    """
+    N3 = disps.shape[1]
+    N = N3 // 3
+    n_basis_fc2 = compress_eigvecs_fc2.shape[1]
+    n_basis_fc3 = compress_eigvecs_fc3.shape[1]
+    n_compr_fc3 = compress_mat_fc3.shape[1]
+    n_basis = n_basis_fc2 + n_basis_fc3
+
+    t1 = time.time()
+    full_basis_fc2 = compress_mat_fc2 @ compress_eigvecs_fc2
+    X2, y_all = get_training_from_full_basis(
+        disps, forces, full_basis_fc2.T.reshape((n_basis_fc2, N, N, 3, 3))
+    )
+    t2 = time.time()
+    print(" training data (fc2):    ", t2 - t1)
+
+    t1 = time.time()
+    """
+    compress_mat_fc3 = (
+        csr_NNN333_to_NN33N3(compress_mat_fc3, N).reshape((NN33, -1)).tocsr()
+    )
+    """
+    """peak memory part (when batch size is less than nearly 50)"""
+    compress_mat_fc3 = -0.5 * reshape_compress_mat(compress_mat_fc3, N)
+    if compress_perm_fc2:
+        c_perm_fc2 = get_perm_compr_matrix(N)
+        compress_mat_fc3 = dot_product_sparse(
+            c_perm_fc2.T, compress_mat_fc3, use_mkl=use_mkl
+        )
+
+    t2 = time.time()
+    print(" precond. compress_mat (for fc3):", t2 - t1)
+
+    sparse_disps = True if use_mkl else False
+    mat33 = np.zeros((n_compr_fc3, n_compr_fc3), dtype=float)
+    mat23 = np.zeros((n_basis_fc2, n_compr_fc3), dtype=float)
+    mat3y = np.zeros(n_compr_fc3, dtype=float)
+    begin_batch, end_batch = get_batch_slice(disps.shape[0], batch_size)
+    for begin, end in zip(begin_batch, end_batch):
+        t01 = time.time()
+        """peak memory part (when batch size is more than nearly 50)"""
+        disps_batch = set_2nd_disps(disps[begin:end], sparse=sparse_disps)
+        if compress_perm_fc2:
+            disps_batch = disps_batch @ c_perm_fc2
+
+        X3 = dot_product_sparse(
+            disps_batch, compress_mat_fc3, use_mkl=use_mkl, dense=True
+        ).reshape((-1, n_compr_fc3))
+        y_batch = forces[begin:end].reshape(-1)
+        mat23 += X2[begin * N3 : end * N3].T @ X3
+        mat33 += X3.T @ X3
+        mat3y += X3.T @ y_batch
+        t02 = time.time()
+        print(" solver_block:", end, ":, t =", t02 - t01)
+
+    XTX = np.zeros((n_basis, n_basis), dtype=float)
+    XTy = np.zeros(n_basis, dtype=float)
+    XTX[:n_basis_fc2, :n_basis_fc2] = X2.T @ X2
+    XTX[:n_basis_fc2, n_basis_fc2:] = mat23 @ compress_eigvecs_fc3
+    XTX[n_basis_fc2:, :n_basis_fc2] = XTX[:n_basis_fc2, n_basis_fc2:].T
+    XTX[n_basis_fc2:, n_basis_fc2:] = (
+        compress_eigvecs_fc3.T @ mat33 @ compress_eigvecs_fc3
+    )
+    XTy[:n_basis_fc2] = X2.T @ y_all
+    XTy[n_basis_fc2:] = compress_eigvecs_fc3.T @ mat3y
+
+    t3 = time.time()
+    print(" (disp @ compr @ eigvecs).T @ (disp @ compr @ eigvecs):", t3 - t2)
+    return XTX, XTy
+
+
+def get_training_no_sum_rule_basis(
+    disps,
+    forces,
+    compress_mat_fc2,
+    compress_mat_fc3,
+    compress_eigvecs_fc2,
+    batch_size=100,
+    use_mkl=False,
+    compress_perm_fc2=False,
+):
+    r"""Calculate X.T @ X and X.T @ y.
+
+    X = displacements @ compress_mat @ compress_eigvecs
+    X = np.hstack([X_fc2, X_fc3])
+
+    displacements (fc2): (n_samples, N3)
+    displacements (fc3): (n_samples, NN33)
+    compress_mat_fc2: (NN33, n_compr)
+    compress_mat_fc3: (NNN333, n_compr_fc3)
+    compress_eigvecs_fc2: (n_compr_fc2, n_basis_fc2)
+    compress_eigvecs_fc3: (n_compr_fc3, n_basis_fc3)
+    Matrix reshapings are appropriately applied to compress_mat
+    and its products.
+
+    X.T @ X and X.T @ y are sequentially calculated using divided dataset.
+    X.T @ X = \sum_i X_i.T @ X_i
+    X.T @ y = \sum_i X_i.T @ y_i (i: batch index)
+    """
+    N3 = disps.shape[1]
+    N = N3 // 3
+    n_basis_fc2 = compress_eigvecs_fc2.shape[1]
+    n_compr_fc3 = compress_mat_fc3.shape[1]
+    n_basis = n_basis_fc2 + n_compr_fc3
+
+    t1 = time.time()
+    full_basis_fc2 = compress_mat_fc2 @ compress_eigvecs_fc2
+    X2, y_all = get_training_from_full_basis(
+        disps, forces, full_basis_fc2.T.reshape((n_basis_fc2, N, N, 3, 3))
+    )
+    t2 = time.time()
+    print(" training data (fc2):    ", t2 - t1)
+
+    t1 = time.time()
+    """peak memory part (when batch size is less than nearly 50)"""
+    compress_mat_fc3 = -0.5 * reshape_compress_mat(compress_mat_fc3, N)
+    if compress_perm_fc2:
+        c_perm_fc2 = get_perm_compr_matrix(N)
+        compress_mat_fc3 = dot_product_sparse(
+            c_perm_fc2.T, compress_mat_fc3, use_mkl=use_mkl
+        )
+    t2 = time.time()
+    print(" precond. compress_mat (for fc3):", t2 - t1)
+
+    sparse_disps = True if use_mkl else False
+    mat33 = np.zeros((n_compr_fc3, n_compr_fc3), dtype=float)
+    mat23 = np.zeros((n_basis_fc2, n_compr_fc3), dtype=float)
+    mat3y = np.zeros(n_compr_fc3, dtype=float)
+    begin_batch, end_batch = get_batch_slice(disps.shape[0], batch_size)
+    for begin, end in zip(begin_batch, end_batch):
+        t01 = time.time()
+        """peak memory part (when batch size is more than 50)"""
+        disps_batch = set_2nd_disps(disps[begin:end], sparse=sparse_disps)
+        if compress_perm_fc2:
+            disps_batch = disps_batch @ c_perm_fc2
+
+        X3 = dot_product_sparse(
+            disps_batch, compress_mat_fc3, use_mkl=use_mkl, dense=True
+        ).reshape((-1, n_compr_fc3))
+        y_batch = forces[begin:end].reshape(-1)
+
+        mat23 += X2[begin * N3 : end * N3].T @ X3
+        mat33 += X3.T @ X3
+        mat3y += X3.T @ y_batch
+        t02 = time.time()
+        print(" solver_block:", end, ":, t =", t02 - t01)
+
+    XTX = np.zeros((n_basis, n_basis), dtype=float)
+    XTy = np.zeros(n_basis, dtype=float)
+    XTX[:n_basis_fc2, :n_basis_fc2] = X2.T @ X2
+    XTX[:n_basis_fc2, n_basis_fc2:] = mat23
+    XTX[n_basis_fc2:, :n_basis_fc2] = XTX[:n_basis_fc2, n_basis_fc2:].T
+    XTX[n_basis_fc2:, n_basis_fc2:] = mat33
+    XTy[:n_basis_fc2] = X2.T @ y_all
+    XTy[n_basis_fc2:] = mat3y
+
+    t3 = time.time()
+    print(" (disp @ compr @ eigvecs).T @ (disp @ compr @ eigvecs):", t3 - t2)
+    return XTX, XTy
+
+
+def run_solver_O2O3(
+    disps,
+    forces,
+    compress_mat_fc2,
+    compress_mat_fc3,
+    compress_eigvecs_fc2,
+    compress_eigvecs_fc3,
+    batch_size=100,
+    use_mkl=False,
+):
+    """Estimate coeffs. in X @ coeffs = y.
+
+    X_fc2 = displacements_fc2 @ compress_mat_fc2 @ compress_eigvecs_fc2
+    X_fc3 = displacements_fc3 @ compress_mat_fc3 @ compress_eigvecs_fc3
+    X = np.hstack([X_fc2, X_fc3])
+
+    Matrix reshapings are appropriately applied.
+    X: features (n_samples * N3, N_basis_fc2 + N_basis_fc3)
+    y: observations (forces), (n_samples * N3)
+
+    """
+    XTX, XTy = get_training(
+        disps,
+        forces,
+        compress_mat_fc2,
+        compress_mat_fc3,
+        compress_eigvecs_fc2,
+        compress_eigvecs_fc3,
+        batch_size=batch_size,
+        use_mkl=use_mkl,
+    )
+    coefs = solve_linear_equation(XTX, XTy)
+    n_basis_fc2 = compress_eigvecs_fc2.shape[1]
+    coefs_fc2, coefs_fc3 = coefs[:n_basis_fc2], coefs[n_basis_fc2:]
+    return coefs_fc2, coefs_fc3
+
+
+def run_solver_O2O3_no_sum_rule_basis(
+    disps,
+    forces,
+    compress_mat_fc2,
+    compress_mat_fc3,
+    compress_eigvecs_fc2,
+    batch_size=100,
+    use_mkl=False,
+):
+    """Estimate coeffs. in X @ coeffs = y.
+
+    X_fc2 = displacements_fc2 @ compress_mat_fc2 @ compress_eigvecs_fc2
+    X_fc3 = displacements_fc3 @ compress_mat_fc3 @ compress_eigvecs_fc3
+    X = np.hstack([X_fc2, X_fc3])
+
+    Matrix reshapings are appropriately applied.
+    X: features (n_samples * N3, N_basis_fc2 + N_basis_fc3)
+    y: observations (forces), (n_samples * N3)
+
+    """
+    XTX, XTy = get_training_no_sum_rule_basis(
+        disps,
+        forces,
+        compress_mat_fc2,
+        compress_mat_fc3,
+        compress_eigvecs_fc2,
+        batch_size=batch_size,
+        use_mkl=use_mkl,
+    )
+    coefs = solve_linear_equation(XTX, XTy)
+    n_basis_fc2 = compress_eigvecs_fc2.shape[1]
+    coefs_fc2, coefs_fc3 = coefs[:n_basis_fc2], coefs[n_basis_fc2:]
+    return coefs_fc2, coefs_fc3
diff --git a/src/symfc/solvers/solver_O3.py b/src/symfc/solvers/solver_O3.py
new file mode 100644
index 0000000..e5ae87b
--- /dev/null
+++ b/src/symfc/solvers/solver_O3.py
@@ -0,0 +1,236 @@
+"""3rd order force constants solver."""
+
+import time
+from typing import Optional
+
+import numpy as np
+from scipy.sparse import csr_array
+
+from symfc.utils.eig_tools import dot_product_sparse
+from symfc.utils.solver_funcs import fit, get_batch_slice, solve_linear_equation
+from symfc.utils.utils_O2 import get_perm_compr_matrix
+
+
+def run_solver_dense_O3(disps, forces, compress_mat, compress_eigvecs):
+    """Estimate coeffs. in X @ coeffs = y.
+
+    X: features (calculated from Kronecker products of displacements)
+                (n_samples * N3, N_basis)
+    y: observations (forces), (n_samples * N3)
+
+    """
+    N3 = disps.shape[1]
+    N = N3 // 3
+
+    t0 = time.time()
+    full_basis = compress_mat @ compress_eigvecs
+    n_basis = full_basis.shape[1]
+    full_basis = full_basis.T.reshape((n_basis, N, N, N, 3, 3, 3))
+    t1 = time.time()
+    print(" elapsed time (compute full basis) =", t1 - t0)
+
+    X, y = _get_training_from_full_basis_set(disps, forces, full_basis)
+    coefs = fit(X, y)
+    t2 = time.time()
+    print(" elapsed time (solve fc3)          =", t2 - t1)
+
+    return coefs
+
+
+def run_solver_sparse_O3(
+    disps, forces, compress_mat, compress_eigvecs, batch_size=200, use_mkl=False
+):
+    """Estimating coeffs. in X @ coeffs = y by solving normal equation.
+
+    (X.T @ X) @ coeffs = X.T @ y
+    X = displacements @ compress_mat @ compress_eigvecs
+    Matrix reshapings are appropriately applied.
+    X: features (n_samples * N3, N_basis)
+    y: observations (forces), (n_samples * N3)
+
+    """
+    XTX, XTy = _get_training_exact(
+        disps,
+        forces,
+        compress_mat,
+        compress_eigvecs,
+        batch_size=batch_size,
+        use_mkl=use_mkl,
+    )
+    coefs = solve_linear_equation(XTX, XTy)
+    return coefs
+
+
+def csr_NNN333_to_NN33N3(mat, N):
+    """Reorder row indices in a sparse matrix (NNN333->NN33N3).
+
+    Return reordered csr_matrix.
+
+    """
+    NNN333, nx = mat.shape
+    mat = mat.tocoo()
+    row = _NNN333_to_NN33N3(mat.row, N)
+    mat = csr_array((mat.data, (row, mat.col)), shape=(NNN333, nx))
+    return mat
+
+
+def csr_NNN333_to_NN33na3(compress_mat: csr_array, N: int, n_a: Optional[int] = None):
+    """Reorder row indices in a sparse matrix.
+
+    (n,N,N,3,3,3) -> (N,3,N,3,n,3)
+    (i,j,k,a,b,c) -> (j,b,k,c,i,a)
+
+    Return reordered csr_matrix.
+
+    n can be any integer, but it is supporsed to be either N or n_a, where n_a
+    is the number of atoms in primitive cell.
+
+    """
+    compress_mat_coo = compress_mat.tocoo()
+    data = compress_mat_coo.data
+    row = compress_mat_coo.row
+    col = compress_mat_coo.col
+
+    if n_a is None:
+        _N = N
+    else:
+        _N = n_a
+
+    NN333 = N * _N * 27
+    N333 = _N * 27
+    N33 = _N * 9
+    N3 = _N * 3
+    conversion_array = np.array(
+        [
+            j * NN333 + k * N333 + m * N33 + n * N3 + i * 3 + l
+            for i, j, k, l, m, n in np.ndindex((_N, N, N, 3, 3, 3))
+        ],
+        dtype=int,
+    )
+    new_row = conversion_array[row]
+    return csr_array((data, (new_row, col)), shape=compress_mat.shape).reshape(
+        N * N * 3 * 3, -1
+    )
+
+
+def set_2nd_disps(disps, sparse=True):
+    """Calculate Kronecker products of displacements.
+
+    disps: (n_samples, N3)
+    disps_2nd: (n_samples, NN33)
+
+    """
+    N = disps.shape[1] // 3
+    disps_2nd = np.zeros((disps.shape[0], 9 * (N**2)))
+    for i, u_vec in enumerate(disps):
+        u2 = np.kron(u_vec, u_vec).reshape((N, 3, N, 3))
+        disps_2nd[i] = u2.transpose((0, 2, 1, 3)).reshape(-1)
+
+    if sparse:
+        return csr_array(disps_2nd)
+    return disps_2nd
+
+
+def _NNN333_to_NN33N3(row, N):
+    """Reorder row indices in a sparse matrix (NNN333->NN33N3)."""
+    i, rem = np.divmod(row, 27 * (N**2))
+    j, rem = np.divmod(rem, 27 * N)
+    k, rem = np.divmod(rem, 27)
+    a, rem = np.divmod(rem, 9)
+    b, c = np.divmod(rem, 3)
+
+    vec = i * 27 * (N**2)
+    vec += j * 27 * N
+    vec += k * 3
+    vec += a * 9 * N
+    vec += b * 3 * N + c
+    return vec
+
+
+def _get_training_from_full_basis_set(disps, forces, full_basis):
+    """Return training dataset (X, y).
+
+    The dataset is transformed from displacements, forces, and full basis-set.
+
+    disps: (n_samples, N3)
+    forces: (n_samples, N3)
+    full_basis: (n_basis, N, N, N, 3, 3, 3)
+
+    X: features (calculated from Kronecker products of displacements)
+                (n_samples * N3, N_basis)
+    y: observations (forces), (n_samples * N3)
+
+    X = displacements @ compress_mat @ compress_eigvecs
+
+    displacements: (n_samples, NN33)
+    compress_mat: (NNN333, n_compr)
+    compress_eigvecs: (n_compr, n_basis)
+    Matrix reshapings are appropriately applied to compress_mat
+    and its products.
+
+    """
+    N3 = full_basis.shape[1] * 3
+    NN33 = N3 * N3
+    n_basis = full_basis.shape[0]
+
+    disps_reshape = set_2nd_disps(disps, sparse=False)
+
+    full_basis = full_basis.transpose((1, 2, 4, 5, 3, 6, 0)).reshape((NN33, -1))
+    X = -0.5 * np.dot(disps_reshape, full_basis)
+    X = X.reshape((-1, n_basis))
+    y = forces.reshape(-1)
+    return X, y
+
+
+def _get_training_exact(
+    disps, forces, compress_mat, compress_eigvecs, batch_size=200, use_mkl=False
+):
+    r"""Calculate X.T @ X and X.T @ y.
+
+    X = displacements @ compress_mat @ compress_eigvecs
+
+    displacements: (n_samples, NN33)
+    compress_mat: (NNN333, n_compr)
+    compress_eigvecs: (n_compr, n_basis)
+    Matrix reshapings are appropriately applied to compress_mat
+    and its products.
+
+    X.T @ X and X.T @ y are sequentially calculated using divided dataset.
+        X.T @ X = \sum_i X_i.T @ X_i
+        X.T @ y = \sum_i X_i.T @ y_i (i: batch index)
+
+    """
+    N3 = disps.shape[1]
+    N = N3 // 3
+    NN33 = N3 * N3
+    n_compr = compress_mat.shape[1]
+
+    t1 = time.time()
+    c_perm_fc2 = get_perm_compr_matrix(N)
+    compress_mat = csr_NNN333_to_NN33N3(compress_mat, N).reshape((NN33, -1)).tocsr()
+    compress_mat = -0.5 * (c_perm_fc2.T @ compress_mat)
+    t2 = time.time()
+    print(" precond. compress_mat (for fc3):", t2 - t1)
+
+    sparse_disps = True if use_mkl else False
+    XTX = np.zeros((n_compr, n_compr), dtype=float)
+    XTy = np.zeros(n_compr, dtype=float)
+    begin_batch, end_batch = get_batch_slice(disps.shape[0], batch_size)
+    for begin, end in zip(begin_batch, end_batch):
+        t01 = time.time()
+        disps_batch = set_2nd_disps(disps[begin:end], sparse=sparse_disps)
+        disps_batch = disps_batch @ c_perm_fc2
+        X = dot_product_sparse(
+            disps_batch, compress_mat, use_mkl=use_mkl, dense=True
+        ).reshape((-1, n_compr))
+        y_batch = forces[begin:end].reshape(-1)
+        XTX += X.T @ X
+        XTy += X.T @ y_batch
+        t02 = time.time()
+        print(" solver_block:", end, ":, t =", t02 - t01)
+
+    XTX = compress_eigvecs.T @ XTX @ compress_eigvecs
+    XTy = compress_eigvecs.T @ XTy
+    t3 = time.time()
+    print(" (disp @ compr @ eigvecs).T @ (disp @ compr @ eigvecs):", t3 - t2)
+    return XTX, XTy
diff --git a/src/symfc/solvers/solver_base.py b/src/symfc/solvers/solver_base.py
new file mode 100644
index 0000000..32e44fa
--- /dev/null
+++ b/src/symfc/solvers/solver_base.py
@@ -0,0 +1,54 @@
+"""Base class of force constants solvers."""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from collections.abc import Sequence
+from typing import Optional, Union
+
+import numpy as np
+
+from symfc.basis_sets import FCBasisSetBase
+
+
+class FCSolverBase(ABC):
+    """Abstract base class of force constants solvers."""
+
+    def __init__(
+        self,
+        basis_set: Union[FCBasisSetBase, Sequence[FCBasisSetBase]],
+        use_mkl: bool = False,
+        log_level: int = 0,
+    ):
+        """Init method."""
+        self._basis_set: FCBasisSetBase = basis_set
+        self._use_mkl: bool = use_mkl
+        self._log_level: int = log_level
+        if isinstance(self._basis_set, Sequence):
+            _basis_set: FCBasisSetBase = self._basis_set[0]
+        else:
+            _basis_set = basis_set
+        self._natom: int = _basis_set.translation_permutations.shape[1]
+        self._coefs: Optional[Union[np.ndarray, Sequence[np.ndarray]]] = None
+
+    @property
+    def coefs(self) -> Optional[np.ndarray]:
+        """Return coefficients of force constants with respect to basis set."""
+        return self._coefs
+
+    @property
+    @abstractmethod
+    def full_fc(self) -> Optional[np.ndarray]:
+        """Return full force constants."""
+        pass
+
+    @property
+    @abstractmethod
+    def compact_fc(self) -> Optional[np.ndarray]:
+        """Return compact force constants."""
+        pass
+
+    @abstractmethod
+    def solve(self):
+        """Solve coefficients of basis set from displacements and forces."""
+        pass
diff --git a/src/symfc/spg_reps/__init__.py b/src/symfc/spg_reps/__init__.py
new file mode 100644
index 0000000..90b8eeb
--- /dev/null
+++ b/src/symfc/spg_reps/__init__.py
@@ -0,0 +1,6 @@
+"""Representations of space group for N-th order force constants."""
+
+from .spg_reps_base import SpgRepsBase  # noqa 401
+from .spg_reps_O1 import SpgRepsO1  # noqa 401
+from .spg_reps_O2 import SpgRepsO2  # noqa 401
+from .spg_reps_O3 import SpgRepsO3  # noqa 401
diff --git a/src/symfc/spg_reps/spg_reps_O1.py b/src/symfc/spg_reps/spg_reps_O1.py
new file mode 100644
index 0000000..42c763e
--- /dev/null
+++ b/src/symfc/spg_reps/spg_reps_O1.py
@@ -0,0 +1,74 @@
+"""O1 reps of space group ops with respect to atomic coordinate basis."""
+
+from __future__ import annotations
+
+import numpy as np
+from scipy.sparse import csr_array
+
+from symfc.utils.utils import SymfcAtoms
+
+from .spg_reps_base import SpgRepsBase
+
+
+class SpgRepsO1(SpgRepsBase):
+    """Class of reps of space group operations for fc1."""
+
+    def __init__(self, supercell: SymfcAtoms):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+
+        """
+        self._r1_reps: list[csr_array]
+        self._col: np.ndarray
+        self._data: np.ndarray
+        super().__init__(supercell)
+
+    @property
+    def r_reps(self) -> list[csr_array]:
+        """Return 1st rank tensor rotation matricies."""
+        return self._r1_reps
+
+    def get_sigma1_rep(self, i: int) -> csr_array:
+        """Compute and return i-th atomic permutation matrix.
+
+        Parameters
+        ----------
+        i : int
+            Index of coset presentations of space group operations.
+
+        """
+        data, row, col, shape = self._get_sigma1_rep_data(i)
+        return csr_array((data, (row, col)), shape=shape)
+
+    def _prepare(self):
+        super()._prepare()
+        N = len(self._numbers)
+        self._col = np.arange(N, dtype=int)
+        self._data = np.ones(N, dtype=int)
+        self._compute_r1_reps()
+
+    def _compute_r1_reps(self, tol: float = 1e-10):
+        """Compute and return 1st rank tensor rotation matricies.
+
+        This is equivalent to rotation matrix of atomic position in Cartesian
+        coordinates.
+
+        """
+        r1_reps = []
+        for r in self._unique_rotations:
+            r1_rep: np.ndarray = self._lattice.T @ r @ np.linalg.inv(self._lattice.T)
+            row, col = np.nonzero(np.abs(r1_rep) > tol)
+            data = r1_rep[(row, col)]
+            r1_reps.append(csr_array((data, (row, col)), shape=r1_rep.shape))
+        self._r1_reps = r1_reps
+
+    def _get_sigma1_rep_data(self, i: int) -> csr_array:
+        uri = self._unique_rotation_indices
+        permutation = self._permutations[uri[i]]
+        N = len(self._numbers)
+        row = permutation
+        return self._data, row, self._col, (N, N)
diff --git a/src/symfc/spg_reps/spg_reps_O2.py b/src/symfc/spg_reps/spg_reps_O2.py
new file mode 100644
index 0000000..1f9bf2f
--- /dev/null
+++ b/src/symfc/spg_reps/spg_reps_O2.py
@@ -0,0 +1,73 @@
+"""O2 reps of space group ops with respect to atomic coordinate basis."""
+
+from __future__ import annotations
+
+import numpy as np
+from scipy.sparse import csr_array
+
+from symfc.utils.utils import SymfcAtoms
+
+from .spg_reps_base import SpgRepsBase
+
+
+class SpgRepsO2(SpgRepsBase):
+    """Class of reps of space group operations for fc2."""
+
+    def __init__(self, supercell: SymfcAtoms):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+
+        """
+        self._r2_reps: list[csr_array]
+        self._col: np.ndarray
+        self._data: np.ndarray
+        super().__init__(supercell)
+
+    @property
+    def r_reps(self) -> list[csr_array]:
+        """Return 2nd rank tensor rotation matricies."""
+        return self._r2_reps
+
+    def get_sigma2_rep(self, i: int) -> csr_array:
+        """Compute and return i-th atomic pair permutation matrix.
+
+        Parameters
+        ----------
+        i : int
+            Index of coset presentations of space group operations.
+
+        """
+        data, row, col, shape = self._get_sigma2_rep_data(i)
+        return csr_array((data, (row, col)), shape=shape)
+
+    def _prepare(self):
+        super()._prepare()
+        N = len(self._numbers)
+        a = np.arange(N)
+        self._atom_pairs = np.stack(np.meshgrid(a, a), axis=-1).reshape(-1, 2)
+        self._coeff = np.array([1, N], dtype=int)
+        self._col = self._atom_pairs @ self._coeff
+        self._data = np.ones(N * N, dtype=int)
+        self._compute_r2_reps()
+
+    def _compute_r2_reps(self, tol: float = 1e-10):
+        """Compute and return 2nd rank tensor rotation matricies."""
+        r2_reps = []
+        for r in self._unique_rotations:
+            r_c = self._lattice.T @ r @ np.linalg.inv(self._lattice.T)
+            r2_rep = np.kron(r_c, r_c)
+            row, col = np.nonzero(np.abs(r2_rep) > tol)
+            data = r2_rep[(row, col)]
+            r2_reps.append(csr_array((data, (row, col)), shape=r2_rep.shape))
+        self._r2_reps = r2_reps
+
+    def _get_sigma2_rep_data(self, i: int) -> csr_array:
+        uri = self._unique_rotation_indices
+        permutation = self._permutations[uri[i]]
+        NN = len(self._numbers) ** 2
+        row = permutation[self._atom_pairs] @ self._coeff
+        return self._data, row, self._col, (NN, NN)
diff --git a/src/symfc/spg_reps/spg_reps_O3.py b/src/symfc/spg_reps/spg_reps_O3.py
new file mode 100644
index 0000000..2d4cab8
--- /dev/null
+++ b/src/symfc/spg_reps/spg_reps_O3.py
@@ -0,0 +1,72 @@
+"""O3 reps of space group ops with respect to atomic coordinate basis."""
+
+from __future__ import annotations
+
+import numpy as np
+from scipy.sparse import csr_array
+
+from symfc.spg_reps import SpgRepsBase
+from symfc.utils.utils import SymfcAtoms
+
+
+class SpgRepsO3(SpgRepsBase):
+    """Class of reps of space group operations for fc3."""
+
+    def __init__(self, supercell: SymfcAtoms):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+
+        """
+        self._r3_reps: list[csr_array]
+        self._col: np.ndarray
+        self._data: np.ndarray
+        super().__init__(supercell)
+
+    @property
+    def r_reps(self) -> list[csr_array]:
+        """Return 3rd rank tensor rotation matricies."""
+        return self._r3_reps
+
+    def get_sigma3_rep(self, i: int) -> csr_array:
+        """Compute and return i-th atomic pair permutation matrix.
+
+        Parameters
+        ----------
+        i : int
+            Index of coset presentations of space group operations.
+
+        """
+        data, row, col, shape = self._get_sigma3_rep_data(i)
+        return csr_array((data, (row, col)), shape=shape)
+
+    def _prepare(self):
+        super()._prepare()
+        N = len(self._numbers)
+        a = np.arange(N)
+        self._atom_triplets = np.stack(np.meshgrid(a, a, a), axis=-1).reshape(-1, 3)
+        self._coeff = np.array([1, N, N**2], dtype=int)
+        self._col = self._atom_triplets @ self._coeff
+        self._data = np.ones(N * N * N, dtype=int)
+        self._compute_r3_reps()
+
+    def _compute_r3_reps(self, tol: float = 1e-10):
+        """Compute and return 3rd rank tensor rotation matricies."""
+        r3_reps = []
+        for r in self._unique_rotations:
+            r_c = self._lattice.T @ r @ np.linalg.inv(self._lattice.T)
+            r3_rep = np.kron(r_c, np.kron(r_c, r_c))
+            row, col = np.nonzero(np.abs(r3_rep) > tol)
+            data = r3_rep[(row, col)]
+            r3_reps.append(csr_array((data, (row, col)), shape=r3_rep.shape))
+        self._r3_reps = r3_reps
+
+    def _get_sigma3_rep_data(self, i: int) -> csr_array:
+        uri = self._unique_rotation_indices
+        permutation = self._permutations[uri[i]]
+        NNN = len(self._numbers) ** 3
+        row = permutation[self._atom_triplets] @ self._coeff
+        return self._data, row, self._col, (NNN, NNN)
diff --git a/src/symfc/spg_reps/spg_reps_base.py b/src/symfc/spg_reps/spg_reps_base.py
new file mode 100644
index 0000000..dd5b1df
--- /dev/null
+++ b/src/symfc/spg_reps/spg_reps_base.py
@@ -0,0 +1,102 @@
+"""Reps base of space group ops with respect to atomic coordinate basis."""
+
+from __future__ import annotations
+
+from typing import Optional
+
+import numpy as np
+import spglib
+
+from symfc.utils.utils import SymfcAtoms, compute_sg_permutations
+
+
+class SpgRepsBase:
+    """Base class of reps of space group operations."""
+
+    def __init__(self, supercell: SymfcAtoms):
+        """Init method.
+
+        Parameters
+        ----------
+        supercell : SymfcAtoms
+            Supercell.
+
+        """
+        self._lattice = np.array(supercell.cell, dtype="double", order="C")
+        self._positions = np.array(
+            supercell.scaled_positions, dtype="double", order="C"
+        )
+        self._numbers = supercell.numbers
+        self._unique_rotations: Optional[np.ndarray] = None
+        self._unique_rotation_indices: Optional[np.ndarray] = None
+        self._translation_permutations: Optional[np.ndarray] = None
+        self._prepare()
+
+    @property
+    def translation_permutations(self) -> np.ndarray:
+        """Return permutations by lattice translation.
+
+        Returns
+        --------
+        Atom indices after lattice translations.
+        shape=(lattice_translations, supercell_atoms), dtype=int
+
+        """
+        return self._translation_permutations
+
+    @property
+    def unique_rotation_indices(self) -> np.ndarray:
+        """Return indices of coset representatives of space group operations."""
+        return self._unique_rotation_indices
+
+    def _prepare(self) -> np.ndarray:
+        rotations, translations = self._get_symops()
+        (
+            self._unique_rotation_indices,
+            self._unique_rotations,
+        ) = self._get_unique_rotation_indices(rotations)
+        self._permutations = compute_sg_permutations(
+            self._positions, rotations, translations, self._lattice.T, 1e-5
+        )
+        self._translation_permutations = self._get_translation_permutations(rotations)
+
+    def _get_translation_permutations(self, rotations) -> np.ndarray:
+        eye3 = np.eye(3, dtype=int)
+        trans_perms = []
+        for r, perm in zip(rotations, self._permutations):
+            if np.array_equal(r, eye3):
+                trans_perms.append(perm)
+        return np.array(trans_perms, dtype="intc", order="C")
+
+    def _get_unique_rotation_indices(self, rotations: np.ndarray) -> list[int]:
+        unique_rotations: list[np.ndarray] = []
+        indices = []
+        for i, r in enumerate(rotations):
+            is_found = False
+            for ur in unique_rotations:
+                if np.array_equal(r, ur):
+                    is_found = True
+                    break
+            if not is_found:
+                unique_rotations.append(r)
+                indices.append(i)
+        return indices, unique_rotations
+
+    def _get_symops(self) -> tuple[np.ndarray, np.ndarray]:
+        """Return symmetry operations.
+
+        The set of inverse operations is the same as the set of the operations.
+
+        Returns
+        -------
+        rotations : array_like
+            A set of rotation matrices of inverse space group operations.
+            (n_symops, 3, 3), dtype='intc', order='C'
+        translations : array_like
+            A set of translation vectors. It is assumed that inverse matrices are
+            included in this set.
+            (n_symops, 3), dtype='double'.
+
+        """
+        symops = spglib.get_symmetry((self._lattice, self._positions, self._numbers))
+        return symops["rotations"], symops["translations"]
diff --git a/src/symfc/utils/__init__.py b/src/symfc/utils/__init__.py
new file mode 100644
index 0000000..dc51750
--- /dev/null
+++ b/src/symfc/utils/__init__.py
@@ -0,0 +1,1 @@
+"""Utilities."""
diff --git a/src/symfc/utils/eig_tools.py b/src/symfc/utils/eig_tools.py
new file mode 100644
index 0000000..b662fe6
--- /dev/null
+++ b/src/symfc/utils/eig_tools.py
@@ -0,0 +1,241 @@
+"""Utility functions for eigenvalue solutions."""
+
+import itertools
+from collections import defaultdict
+
+import numpy as np
+import scipy
+from scipy.linalg import get_lapack_funcs
+from scipy.sparse import csr_array
+
+try:
+    from sparse_dot_mkl import dot_product_mkl
+except ImportError:
+    pass
+
+
+def dot_product_sparse(
+    A: csr_array,
+    B: csr_array,
+    use_mkl: bool = False,
+    dense: bool = False,
+) -> csr_array:
+    """Compute dot-product of sparse matrices."""
+    if use_mkl:
+        return dot_product_mkl(A, B, dense=dense)
+    return A @ B
+
+
+def eigsh_projector(
+    p: csr_array, verbose: bool = True, log_interval: int = 10000
+) -> csr_array:
+    """Solve eigenvalue problem for matrix p.
+
+    Return sparse matrix for eigenvectors of matrix p.
+
+    This algorithm begins with finding block diagonal structures in matrix p.
+    For each block submatrix, eigenvectors are solved. When p = diag(A,B), Av =
+    v, and Bw = w, p[v,0] = [v,0] and p[0,w] = [0,w] are solutions.
+
+    This function avoids solving numpy.eigh for duplicate block matrices. This
+    function is efficient for matrix p composed of many duplicate block
+    matrices.
+
+    """
+    n_components, labels = scipy.sparse.csgraph.connected_components(p)
+    group = defaultdict(list)
+    for i, l in enumerate(labels):
+        group[l].append(i)
+
+    r = np.array([i for ids in group.values() for i in ids for j in ids])
+    c = np.array([j for ids in group.values() for i in ids for j in ids])
+    s_end = np.cumsum([len(ids) * len(ids) for ids in group.values()])
+    s_begin = np.zeros(len(s_end), dtype=int)
+    s_begin[1:] = s_end[:-1]
+    p_data = np.ravel(p[r, c])
+
+    if verbose:
+        print(" N (blocks) =", n_components)
+        rank = int(round(sum(p.diagonal())))
+        print(" rank (P) =", rank)
+
+    uniq_eigvecs = dict()
+    row, col, data = [], [], []
+    col_id = 0
+    for i, (s1, s2, ids) in enumerate(zip(s_begin, s_end, group.values())):
+        if verbose and (i + 1) % log_interval == 0:
+            print(" eigsh_block:", i + 1)
+
+        p_block1 = p_data[s1:s2]
+        if len(ids) > 1:
+            key = tuple(p_block1)
+            try:
+                eigvecs = uniq_eigvecs[key]
+            except KeyError:
+                size = len(ids)
+                p_numpy = p_block1.reshape((size, size))
+                rank = int(round(np.trace(p_numpy)))
+                if rank == 0:
+                    eigvals, eigvecs = None, None
+                else:
+                    eigvals, eigvecs = np.linalg.eigh(p_numpy)
+                    nonzero = np.isclose(eigvals, 1.0)
+                    eigvecs = eigvecs[:, nonzero]
+                uniq_eigvecs[key] = eigvecs
+            if eigvecs is not None:
+                n_row, n_col = eigvecs.shape
+                row.extend([ids[i] for i in range(n_row) for j in range(n_col)])
+                col.extend(
+                    [j for i in range(n_row) for j in range(col_id, col_id + n_col)]
+                )
+                data.extend(eigvecs.reshape(-1))
+                col_id += n_col
+        else:
+            if not np.isclose(p_block1[0], 0.0):
+                row.append(ids[0])
+                col.append(col_id)
+                data.append(1.0)
+                col_id += 1
+
+    n_col = col_id
+    return csr_array((data, (row, col)), shape=(p.shape[0], n_col))
+
+
+def eigsh_projector_sumrule(p: csr_array, verbose: bool = True) -> np.ndarray:
+    """Solve eigenvalue problem for matrix p.
+
+    Return dense matrix for eigenvectors of matrix p.
+
+    This algorithm begins with finding block diagonal structures in matrix p.
+    For each block submatrix, eigenvectors are solved.
+    When p = diag(A,B), Av = v, and Bw = w, p[v,0] = [v,0] and p[0,w] = [0,w]
+    are solutions.
+
+    This function solves numpy.eigh for all block matrices.
+    This function is efficient for matrix p composed of nonequivalent
+    block matrices.
+
+    """
+    n_components, labels = scipy.sparse.csgraph.connected_components(p)
+    group = defaultdict(list)
+    for i, l in enumerate(labels):
+        group[l].append(i)
+
+    if verbose:
+        print(" n_blocks in P =", n_components)
+
+    eigvecs_full = np.zeros(p.shape, dtype="double")
+    col_id = 0
+    for i, ids in enumerate(group.values()):
+        if verbose:
+            print(" eigsh_block:", i, ": block_size =", len(ids))
+
+        p_block = p[np.ix_(ids, ids)].toarray()
+        rank = int(round(np.trace(p_block)))
+        if rank > 0:
+            if rank < 30000:
+                eigvals, eigvecs = np.linalg.eigh(p_block)
+            else:
+                print("eigh: solver dsyevr is used.")
+                (syevr,) = get_lapack_funcs(("syevr",), ilp64=False)
+                eigvals, eigvecs, _, _, _ = syevr(p_block, compute_v=True)
+
+                """
+                print('eigh: solver dsyevd is used.')
+                (syevd,) = get_lapack_funcs(('syevd',), ilp64=False)
+                eigvals, eigvecs, _ = syevd(p_block, compute_v=True)
+                """
+                # eigvals, eigvecs, _, _, _ = scipy.linalg.lapack.dsyevr(
+                #       p_block, compute_v=True)
+                # eigvals, eigvecs, _ = scipy.linalg.lapack.dsyevd(
+                #       p_block, compute_v=True)
+
+            nonzero = np.isclose(eigvals, 1.0)
+            eigvecs = eigvecs[:, nonzero]
+            col_ids = np.arange(col_id, col_id + eigvecs.shape[1])
+            eigvecs_full[np.ix_(ids, col_ids)] = eigvecs
+            col_id += eigvecs.shape[1]
+
+    return eigvecs_full[:, :col_id]
+
+
+def connected_components(p: csr_array) -> np.ndarray:
+    """Find connected matrix elements.
+
+    This algorithm is simple but inefficient.
+
+    todo: A breadth-first search or depth-first search should be implemented.
+
+    """
+    labels = np.arange(p.shape[0], dtype=np.int64)
+    for i in range(p.shape[0]):
+        if labels[i] == i:
+            adj = p.getcol(i).indices
+            labels[adj] = i
+    return labels
+
+
+def eigsh_projector_memory_efficient(
+    p: csr_array, verbose: bool = True, log_interval: int = 10000
+):
+    """Eigenvalue solver for connected matrix elements.
+
+    Sparse matrix p must be projector.
+
+    """
+    if len(p.data) < 2147483647:
+        _, labels = scipy.sparse.csgraph.connected_components(p)
+    else:
+        labels = connected_components(p)
+
+    group = defaultdict(list)
+    for i, l in enumerate(labels):
+        group[l].append(i)
+
+    if verbose:
+        print(" N (blocks) =", len(group.keys()))
+        rank = int(round(sum(p.diagonal())))
+        print(" rank (P) =", rank)
+
+    uniq_eigvecs = dict()
+    row, col, data = [], [], []
+    col_id = 0
+    for i, (key, ids) in enumerate(group.items()):
+        if verbose and (i + 1) % log_interval == 0:
+            print(" eigsh_block:", i + 1)
+
+        if len(ids) > 1:
+            r, c = np.array(list(zip(*itertools.product(ids, ids))))
+            p_block1 = np.ravel(p[r, c])
+            key = tuple(p_block1)
+            if key in uniq_eigvecs:
+                eigvecs = uniq_eigvecs[key]
+            else:
+                size = len(ids)
+                p_numpy = p_block1.reshape((size, size))
+                rank = int(round(np.trace(p_numpy)))
+                if rank == 0:
+                    eigvals, eigvecs = None, None
+                else:
+                    eigvals, eigvecs = np.linalg.eigh(p_numpy)
+                    nonzero = np.isclose(eigvals, 1.0)
+                    eigvecs = eigvecs[:, nonzero]
+                uniq_eigvecs[key] = eigvecs
+
+            if eigvecs is not None:
+                n_row, n_col = eigvecs.shape
+                row.extend([ids[i] for i in range(n_row) for j in range(n_col)])
+                col.extend(
+                    [j for i in range(n_row) for j in range(col_id, col_id + n_col)]
+                )
+                data.extend(eigvecs.reshape(-1))
+                col_id = n_col
+        else:
+            if not np.isclose(p[ids[0], ids[0]], 0.0):
+                row.append(ids[0])
+                col.append(col_id)
+                data.append(1.0)
+                col_id += 1
+
+    n_col = col_id
+    return csr_array((data, (row, col)), shape=(p.shape[0], n_col))
diff --git a/src/symfc/utils/matrix_tools_O1.py b/src/symfc/utils/matrix_tools_O1.py
new file mode 100644
index 0000000..f280da2
--- /dev/null
+++ b/src/symfc/utils/matrix_tools_O1.py
@@ -0,0 +1,50 @@
+"""Matrix utility functions for 1st order force constants."""
+
+import numpy as np
+import scipy
+from scipy.sparse import csr_array
+
+
+def compressed_projector_sum_rules(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+):
+    """Return projection matrix for sum rule compressed by C."""
+    proj_cplmt = _compressed_complement_projector_sum_rules(
+        compress_mat, N, use_mkl=use_mkl
+    )
+    return scipy.sparse.identity(proj_cplmt.shape[0]) - proj_cplmt
+
+
+def _compressed_complement_projector_sum_rules_algo1(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+):
+    r"""Return complementary projection matrix for sum rule compressed by C.
+
+    proj_sum_cplmt = [C.T @ Csum(c)] @ [Csum(c).T @ C]
+                   = c_sum_cplmt_compr.T @ c_sum_cplmt_compr
+    Matrix shape of proj_sum_cplmt is (C.shape[1], C.shape[1]).
+    C.shape[0] must be equal to NN33.
+
+    Sum rules are given as sums over i: \sum_i \phi_{ia,jb} = 0
+
+    """
+    N3 = 3 * N
+
+    row = np.arange(N3)
+    col = np.tile(range(3), N)
+    data = np.zeros(N3)
+    data[:] = 1 / np.sqrt(N)
+    c_sum_cplmt = csr_array((data, (row, col)), shape=(N3, 3))
+
+    c_sum_cplmt_compr = c_sum_cplmt.T @ compress_mat
+    proj_sum_cplmt = c_sum_cplmt_compr.T @ c_sum_cplmt_compr
+    return proj_sum_cplmt
+
+
+def _compressed_complement_projector_sum_rules(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+):
+    """Return complementary projection matrix for sum rule compressed by C."""
+    return _compressed_complement_projector_sum_rules_algo1(
+        compress_mat, N, use_mkl=use_mkl
+    )
diff --git a/src/symfc/utils/matrix_tools_O2.py b/src/symfc/utils/matrix_tools_O2.py
new file mode 100644
index 0000000..b04bf2d
--- /dev/null
+++ b/src/symfc/utils/matrix_tools_O2.py
@@ -0,0 +1,61 @@
+"""Matrix utility functions for 2nd order force constants."""
+
+import numpy as np
+import scipy
+from scipy.sparse import csr_array
+
+from symfc.utils.eig_tools import dot_product_sparse
+
+
+def compressed_projector_sum_rules(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+):
+    """Return projection matrix for sum rule compressed by C."""
+    proj_cplmt = _compressed_complement_projector_sum_rules(
+        compress_mat, N, use_mkl=use_mkl
+    )
+    return scipy.sparse.identity(proj_cplmt.shape[0]) - proj_cplmt
+
+
+def _compressed_complement_projector_sum_rules_algo1(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+):
+    r"""Return complementary projection matrix for sum rule compressed by C.
+
+    proj_sum_cplmt = [C.T @ Csum(c)] @ [Csum(c).T @ C]
+                   = c_sum_cplmt_compr.T @ c_sum_cplmt_compr
+    Matrix shape of proj_sum_cplmt is (C.shape[1], C.shape[1]).
+    C.shape[0] must be equal to NN33.
+
+    Sum rules are given as sums over i: \sum_i \phi_{ia,jb} = 0
+
+    """
+    NN33 = 9 * N**2
+    N33 = 9 * N
+
+    row = np.arange(NN33)
+    col = np.tile(range(N33), N)
+    data = np.zeros(NN33)
+    data[:] = 1 / np.sqrt(N)
+    c_sum_cplmt = csr_array((data, (row, col)), shape=(NN33, N33))
+
+    if use_mkl:
+        compress_mat = compress_mat.tocsr()
+        c_sum_cplmt = c_sum_cplmt.tocsr()
+
+    # bottleneck part
+    c_sum_cplmt_compr = dot_product_sparse(c_sum_cplmt.T, compress_mat, use_mkl=use_mkl)
+    proj_sum_cplmt = dot_product_sparse(
+        c_sum_cplmt_compr.T, c_sum_cplmt_compr, use_mkl=use_mkl
+    )
+    # bottleneck part: end
+    return proj_sum_cplmt
+
+
+def _compressed_complement_projector_sum_rules(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+):
+    """Return complementary projection matrix for sum rule compressed by C."""
+    return _compressed_complement_projector_sum_rules_algo1(
+        compress_mat, N, use_mkl=use_mkl
+    )
diff --git a/src/symfc/utils/matrix_tools_O3.py b/src/symfc/utils/matrix_tools_O3.py
new file mode 100644
index 0000000..57a5d6b
--- /dev/null
+++ b/src/symfc/utils/matrix_tools_O3.py
@@ -0,0 +1,342 @@
+"""Matrix utility functions for 3rd order force constants."""
+
+import itertools
+import math
+
+import numpy as np
+import scipy
+from scipy.sparse import csr_array
+
+from symfc.utils.eig_tools import dot_product_sparse
+from symfc.utils.solver_funcs import get_batch_slice
+from symfc.utils.utils_O3 import get_lat_trans_decompr_indices_O3
+
+
+def N3N3N3_to_NNN333(combinations_perm: np.ndarray, N: int) -> np.ndarray:
+    """Transform index order."""
+    vec = combinations_perm[:, 0] // 3 * 27 * N**2
+    vec += combinations_perm[:, 1] // 3 * 27 * N
+    vec += combinations_perm[:, 2] // 3 * 27
+    vec += combinations_perm[:, 0] % 3 * 9
+    vec += combinations_perm[:, 1] % 3 * 3
+    vec += combinations_perm[:, 2] % 3
+    return vec
+
+
+def get_perm_compr_matrix_O3(natom: int) -> csr_array:
+    """Return compression matrix by permutation symmetry.
+
+    Matrix shape is (NNN333, (N*3)(N*3+1)(N*3+2)/6).
+
+    """
+    NNN333 = 27 * natom**3
+    combinations3 = np.array(
+        list(itertools.combinations(range(3 * natom), 3)), dtype=int
+    )
+    combinations2 = np.array(
+        list(itertools.combinations(range(3 * natom), 2)), dtype=int
+    )
+    combinations1 = np.array([[i, i, i] for i in range(3 * natom)], dtype=int)
+
+    n_col3 = combinations3.shape[0]
+    n_col2 = combinations2.shape[0] * 2
+    n_col1 = combinations1.shape[0]
+    n_col = n_col3 + n_col2 + n_col1
+    n_data3 = combinations3.shape[0] * 6
+    n_data2 = combinations2.shape[0] * 6
+    n_data1 = combinations1.shape[0]
+    n_data = n_data3 + n_data2 + n_data1
+
+    row = np.zeros(n_data, dtype="int_")
+    col = np.zeros(n_data, dtype="int_")
+    data = np.zeros(n_data, dtype="double")
+
+    # (3) for FC3 with three distinguished indices (ia,jb,kc)
+    begin_id, end_id = 0, n_data3
+    perms = [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]]
+    combinations_perm = combinations3[:, perms].reshape((-1, 3))
+    row[begin_id:end_id] = N3N3N3_to_NNN333(combinations_perm, natom)
+    col[begin_id:end_id] = np.repeat(range(n_col3), 6)
+    data[begin_id:end_id] = 1 / math.sqrt(6)
+
+    # (2) for FC3 with two distinguished indices (ia,ia,jb)
+    begin_id = end_id
+    end_id = begin_id + n_data2
+    perms = [[0, 0, 1], [0, 1, 0], [1, 0, 0], [0, 1, 1], [1, 0, 1], [1, 1, 0]]
+    combinations_perm = combinations2[:, perms].reshape((-1, 3))
+    row[begin_id:end_id] = N3N3N3_to_NNN333(combinations_perm, natom)
+    col[begin_id:end_id] = np.repeat(range(n_col3, n_col3 + n_col2), 3)
+    data[begin_id:end_id] = 1 / math.sqrt(3)
+
+    # (1) for FC3 with single index ia
+    begin_id = end_id
+    row[begin_id:] = N3N3N3_to_NNN333(combinations1, natom)
+    col[begin_id:] = np.array(range(n_col3 + n_col2, n_col))
+    data[begin_id:] = 1.0
+
+    return csr_array((data, (row, col)), shape=(NNN333, n_col))
+
+
+def compressed_complement_projector_sum_rules_algo1(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+) -> csr_array:
+    r"""Return complementary projection matrix for sum rule compressed by C.
+
+    C = compress_mat.
+
+    proj_sum_cplmt = [C.T @ Csum(c)] @ [Csum(c).T @ C]
+                   = c_sum_cplmt_compr.T @ c_sum_cplmt_compr
+    Matrix shape of proj_sum_cplmt is (C.shape[1], C.shape[1]).
+    C.shape[0] must be equal to NNN333.
+
+    Sum rules are given as sums over i: \sum_i \phi_{ia,jb,kc} = 0
+
+    """
+    NNN333 = 27 * N**3
+    NN333 = 27 * N**2
+
+    row = np.arange(NNN333)
+    col = np.tile(range(NN333), N)
+    data = np.zeros(NNN333)
+    data[:] = 1 / math.sqrt(N)
+    c_sum_cplmt = csr_array((data, (row, col)), shape=(NNN333, NN333))
+
+    if use_mkl:
+        compress_mat = compress_mat.tocsr()
+        c_sum_cplmt = c_sum_cplmt.tocsr()
+
+    # bottleneck part
+    c_sum_cplmt_compr = dot_product_sparse(c_sum_cplmt.T, compress_mat, use_mkl=use_mkl)
+    proj_sum_cplmt = dot_product_sparse(
+        c_sum_cplmt_compr.T, c_sum_cplmt_compr, use_mkl=use_mkl
+    )
+    # bottleneck part: end
+    return proj_sum_cplmt
+
+
+def compressed_complement_projector_sum_rules(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+) -> csr_array:
+    """Return complementary projection matrix for sum rule compressed by C."""
+    return compressed_complement_projector_sum_rules_algo1(
+        compress_mat, N, use_mkl=use_mkl
+    )
+
+
+def compressed_projector_sum_rules(
+    compress_mat: csr_array, N: int, use_mkl: bool = False
+) -> csr_array:
+    """Return projection matrix for sum rule compressed by C."""
+    proj_cplmt = compressed_complement_projector_sum_rules(
+        compress_mat, N, use_mkl=use_mkl
+    )
+    return scipy.sparse.identity(proj_cplmt.shape[0]) - proj_cplmt
+
+
+"""
+Memory-efficient updated functions for running FCBasisSetO3.
+But they have not been yet implemented in FCBasisSetO3.
+"""
+
+
+def set_complement_sum_rules(trans_perms) -> csr_array:
+    """Calculate a decomposition of complementary projector for sum rules.
+
+    It is compressed by C_trans without allocating C_trans.
+    P_sum^(c) = C_trans.T @ C_sum^(c) @ C_sum^(c).T @ C_trans
+              = [C_sum^(c).T @ C_trans].T @ [C_sum^(c).T @ C_trans]
+
+    Return
+    ------
+    Product of C_sum^(c).T and C_trans.
+    """
+    n_lp, natom = trans_perms.shape
+    NNN27 = natom**3 * 27
+    NN27 = natom**2 * 27
+
+    decompr_idx = get_lat_trans_decompr_indices_O3(trans_perms)
+    decompr_idx = decompr_idx.reshape((natom, NN27)).T.reshape(-1)
+
+    row = np.repeat(np.arange(NN27), natom)
+    c_sum_cplmt = csr_array(
+        (
+            np.ones(NNN27, dtype="double"),
+            (row, decompr_idx),
+        ),
+        shape=(NN27, NNN27 // n_lp),
+        dtype="double",
+    )
+    c_sum_cplmt /= n_lp * natom
+    return c_sum_cplmt
+
+
+def compressed_complement_projector_sum_rules_from_compact_compr_mat(
+    trans_perms, n_a_compress_mat: csr_array, use_mkl: bool = False
+) -> csr_array:
+    """Calculate a complementary projector for sum rules.
+
+    This is compressed by C_trans and n_a_compress_mat without
+    allocating C_trans. Batch calculations are used to reduce
+    the amount of memory allocation.
+
+    Return
+    ------
+    Compressed projector
+    P^(c) = n_a_compress_mat.T @ C_trans.T @ C_sum^(c)
+            @ C_sum^(c).T @ C_trans @ n_a_compress_mat
+    """
+    n_lp, natom = trans_perms.shape
+    NNN27 = natom**3 * 27
+    NN27 = natom**2 * 27
+
+    proj_size = n_a_compress_mat.shape[1]
+    proj_sum_cplmt = csr_array(
+        ([], ([], [])),
+        shape=(proj_size, proj_size),
+        dtype="double",
+    )
+
+    decompr_idx = get_lat_trans_decompr_indices_O3(trans_perms)
+    decompr_idx = decompr_idx.reshape((natom, NN27)).T.reshape(-1)
+
+    n_batch = 9
+    if n_batch == 3:
+        batch_size_vector = natom**3 * 9
+        batch_size_matrix = natom**2 * 9
+    elif n_batch == 9:
+        batch_size_vector = natom**3 * 3
+        batch_size_matrix = natom**2 * 3
+    elif n_batch == 27:
+        batch_size_vector = natom**3
+        batch_size_matrix = natom**2
+    elif n_batch == natom:
+        batch_size_vector = natom**2 * 27
+        batch_size_matrix = natom * 27
+    else:
+        raise ValueError("n_batch = 9, 27, or N.")
+
+    for begin, end in zip(*get_batch_slice(NNN27, batch_size_vector)):
+        print("Proj_complement (sum.T @ trans) batch:", end)
+        c_sum_cplmt = csr_array(
+            (
+                np.ones(batch_size_vector, dtype="double"),
+                (
+                    np.repeat(np.arange(batch_size_matrix), natom),
+                    decompr_idx[begin:end],
+                ),
+            ),
+            shape=(batch_size_matrix, NNN27 // n_lp),
+            dtype="double",
+        )
+        c_sum_cplmt = dot_product_sparse(c_sum_cplmt, n_a_compress_mat, use_mkl=use_mkl)
+        proj_sum_cplmt += dot_product_sparse(
+            c_sum_cplmt.T, c_sum_cplmt, use_mkl=use_mkl
+        )
+
+    proj_sum_cplmt /= n_lp * natom
+    return proj_sum_cplmt
+
+
+def compressed_projector_sum_rules_from_compact_compr_mat(
+    trans_perms, n_a_compress_mat: csr_array, use_mkl: bool = False
+) -> csr_array:
+    """Return projection matrix for sum rule.
+
+    This is compressed by C_compr = C_trans @ n_a_compress_mat.
+    """
+    proj_cplmt = compressed_complement_projector_sum_rules_from_compact_compr_mat(
+        trans_perms, n_a_compress_mat, use_mkl=use_mkl
+    )
+    return scipy.sparse.identity(proj_cplmt.shape[0]) - proj_cplmt
+
+
+def get_combinations(n, r):
+    """Return numpy array of combinations.
+
+    combinations = np.array(
+       list(itertools.combinations(range(n), r)), dtype=int
+    )
+    """
+    combs = np.ones((r, n - r + 1), dtype=int)
+    combs[0] = np.arange(n - r + 1)
+    for j in range(1, r):
+        reps = (n - r + j) - combs[j - 1]
+        combs = np.repeat(combs, reps, axis=1)
+        ind = np.add.accumulate(reps)
+        combs[j, ind[:-1]] = 1 - reps[1:]
+        combs[j, 0] = j
+        combs[j] = np.add.accumulate(combs[j])
+    return combs.T
+
+
+def projector_permutation_lat_trans(trans_perms, n_batch=6, use_mkl=False):
+    """Calculate a projector for permutation rules.
+
+    This is compressed by C_trans without allocating C_trans and C_perm.
+    Batch calculations are used to reduce memory allocation.
+
+    Return
+    ------
+    Compressed projector for permutation
+    P_pt = C_trans.T @ C_perm @ C_perm.T @ C_trans
+    """
+    n_lp, natom = trans_perms.shape
+    decompr_idx = get_lat_trans_decompr_indices_O3(trans_perms)
+    NNN27 = natom**3 * 27
+
+    # (1) for FC3 with single index ia
+    combinations = np.array([[i, i, i] for i in range(3 * natom)], dtype=int)
+    n_perm1 = combinations.shape[0]
+
+    combinations = N3N3N3_to_NNN333(combinations, natom)
+    c_pt = csr_array(
+        (
+            np.full(n_perm1, 1.0 / math.sqrt(n_lp)),
+            (np.arange(n_perm1), decompr_idx[combinations]),
+        ),
+        shape=(n_perm1, NNN27 // n_lp),
+        dtype="double",
+    )
+    proj_pt = dot_product_sparse(c_pt.T, c_pt, use_mkl=use_mkl)
+
+    # (2) for FC3 with two distinguished indices (ia,ia,jb)
+    # combinations = np.stack(np.triu_indices(3 * natom, k=1), axis=-1)
+    combinations = get_combinations(3 * natom, 2)
+    n_perm2 = combinations.shape[0] * 2
+
+    perms = [[0, 0, 1], [0, 1, 0], [1, 0, 0], [0, 1, 1], [1, 0, 1], [1, 1, 0]]
+    combinations = combinations[:, perms].reshape((-1, 3))
+    combinations = N3N3N3_to_NNN333(combinations, natom)
+
+    c_pt = csr_array(
+        (
+            np.full(n_perm2 * 3, 1 / math.sqrt(3 * n_lp)),
+            (np.repeat(range(n_perm2), 3), decompr_idx[combinations]),
+        ),
+        shape=(n_perm2, NNN27 // n_lp),
+        dtype="double",
+    )
+    proj_pt += dot_product_sparse(c_pt.T, c_pt, use_mkl=use_mkl)
+
+    # (3) for FC3 with three distinguished indices (ia,jb,kc)
+    combinations = get_combinations(3 * natom, 3)
+    n_perm3 = combinations.shape[0]
+
+    perms = [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]]
+    for begin, end in zip(*get_batch_slice(n_perm3, n_perm3 // n_batch)):
+        print("Proj (perm.T @ trans) batch:", end)
+        batch_size = end - begin
+        combinations_perm = combinations[begin:end][:, perms].reshape((-1, 3))
+        combinations_perm = N3N3N3_to_NNN333(combinations_perm, natom)
+
+        c_pt = csr_array(
+            (
+                np.full(batch_size * 6, 1 / math.sqrt(6 * n_lp)),
+                (np.repeat(np.arange(batch_size), 6), decompr_idx[combinations_perm]),
+            ),
+            shape=(batch_size, NNN27 // n_lp),
+            dtype="double",
+        )
+        proj_pt += dot_product_sparse(c_pt.T, c_pt, use_mkl=use_mkl)
+
+    return proj_pt
diff --git a/src/symfc/utils/solver_funcs.py b/src/symfc/utils/solver_funcs.py
new file mode 100644
index 0000000..d99199f
--- /dev/null
+++ b/src/symfc/utils/solver_funcs.py
@@ -0,0 +1,41 @@
+"""Solver functions."""
+
+import numpy as np
+from scipy.linalg.lapack import get_lapack_funcs
+
+
+def solve_linear_equation(A: np.ndarray, b: np.ndarray):
+    """Solve linear equation using lapack in scipy.
+
+    numpy and scipy implementations
+    x = np.linalg.solve(A, b)
+    x = scipy.linalg.solve(A, b, check_finite=False, assume_a='pos')
+
+    """
+    (posv,) = get_lapack_funcs(("posv",), (A, b))
+    _, x, _ = posv(A, b, lower=False, overwrite_a=False, overwrite_b=False)
+    return x
+
+
+def fit(X: np.ndarray, y: np.ndarray):
+    """Solve a normal equation in least-squares.
+
+    (X.T @ X) @ coefs = X.T @ y
+
+    n_samples, n_features = X.shape
+
+    """
+    A = np.dot(X.T, X)
+    Xy = np.dot(X.T, y)
+    coefs = solve_linear_equation(A, Xy)
+    return coefs
+
+
+def get_batch_slice(n_data, batch_size):
+    """Calculate slice indices for a given batch size."""
+    begin_batch = list(range(0, n_data, batch_size))
+    if len(begin_batch) > 1:
+        end_batch = list(begin_batch[1:]) + [n_data]
+    else:
+        end_batch = [n_data]
+    return begin_batch, end_batch
diff --git a/src/symfc/utils/utils.py b/src/symfc/utils/utils.py
new file mode 100644
index 0000000..729b79b
--- /dev/null
+++ b/src/symfc/utils/utils.py
@@ -0,0 +1,164 @@
+"""Utility functions."""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+def get_indep_atoms_by_lat_trans(trans_perms: np.ndarray) -> np.ndarray:
+    """Return independent atoms by lattice translation symmetry.
+
+    Parameters
+    ----------
+    trans_perms : np.ndarray
+        Atom indices after lattice translations.
+        shape=(lattice_translations, supercell_atoms)
+
+    Returns
+    -------
+    np.ndarray
+        Independent atoms.
+        shape=(n_indep_atoms_by_lattice_translation,), dtype=int
+
+    """
+    unique_atoms: list[int] = []
+    assert np.array_equal(trans_perms[0, :], range(trans_perms.shape[1]))
+    for i, perms in enumerate(trans_perms.T):
+        is_found = False
+        for j in unique_atoms:
+            if j in perms:
+                is_found = True
+                break
+        if not is_found:
+            unique_atoms.append(i)
+    return np.array(unique_atoms, dtype=int)
+
+
+def compute_sg_permutations(
+    positions: np.ndarray,
+    rotations: np.ndarray,
+    translations: np.ndarray,
+    lattice: np.ndarray,
+    symprec: float = 1e-5,
+) -> np.ndarray:
+    """Compute permutations of atoms by space group operations in supercell.
+
+    Permutations of atoms of pure translations and coset representatives are
+    first computed. Then permutations of atoms of all the given space group
+    operations are made as the combitation of these two permutations.
+
+    Parameters
+    ----------
+    positions : ndarray
+        Fractional positions (like SymfcAtoms.scaled_positions) before applying
+        the space group operation.
+    rotations : ndarray
+        Matrix (rotation) parts of space group operations.
+        shape=(len(operations), 3, 3), dtype='intc'
+    translations : ndarray
+        Vector (translation) parts of space group operations.
+        shape=(len(operations), 3), dtype='double'
+    lattice : ndarray
+        Basis vectors in column vectors (like SymfcAtoms.cell.T).
+    symprec : float
+        Symmetry tolerance of the distance unit.
+
+    Returns
+    -------
+    perms : ndarray
+        shape=(len(translations), len(positions)), dtype='intc', order='C'
+
+    """
+    trans_perms = []
+    pure_trans = []
+    for i, (r, t) in enumerate(zip(rotations, translations)):
+        if (r != np.eye(3, dtype=int)).any():
+            continue
+        trans_positions = positions + t
+        diffs = positions[None, :, :] - trans_positions[:, None, :]
+        diffs -= np.rint(diffs)
+        dists = np.linalg.norm(diffs @ lattice.T, axis=2)
+        rows, cols = np.where(dists < symprec)
+        assert len(positions) == len(np.unique(rows)) == len(np.unique(cols))
+        trans_perms.append(cols[np.argsort(rows)])
+        pure_trans.append(t)
+    trans_perms = np.array(trans_perms, dtype=int)
+    pure_trans = np.array(pure_trans)
+
+    unique_r = []
+    unique_t = []
+    r2ur = []
+    unique_rotated_positions = []
+    for i, (r, t) in enumerate(zip(rotations, translations)):
+        is_found = False
+        for j, ur in enumerate(unique_r):
+            if (r == ur).all():
+                is_found = True
+                r2ur.append(j)
+                break
+        if not is_found:
+            r2ur.append(len(unique_r))
+            unique_r.append(r)
+            unique_t.append(t)
+            unique_rotated_positions.append(positions @ r.T + t)
+
+    unique_rotation_perms = []
+    for i, rotated_positions in enumerate(unique_rotated_positions):
+        diffs = positions[None, :, :] - rotated_positions[:, None, :]
+        diffs -= np.rint(diffs)
+        dists = np.linalg.norm(diffs @ lattice.T, axis=2)
+        rows, cols = np.where(dists < symprec)
+        assert len(positions) == len(np.unique(rows)) == len(np.unique(cols))
+        unique_rotation_perms.append(cols[np.argsort(rows)])
+    unique_rotation_perms = np.array(unique_rotation_perms, dtype=int)
+
+    out = []
+    for i, t in enumerate(translations):
+        perms = unique_rotation_perms[r2ur[i]]
+        lattice_trans = t - unique_t[r2ur[i]]
+        diffs = pure_trans - lattice_trans
+        diffs -= np.rint(diffs)
+        dists = np.linalg.norm(diffs @ lattice.T, axis=1)
+        lat_trans_idx = np.where(dists < symprec)
+        assert len(lat_trans_idx) == 1
+        out.append(trans_perms[lat_trans_idx[0], perms])
+    out = np.array(out, dtype="intc", order="C")
+    return out
+
+
+class SymfcAtoms:
+    """Class to represent crystal structure mimicing PhonopyAtoms."""
+
+    def __init__(
+        self,
+        numbers=None,
+        scaled_positions=None,
+        cell=None,
+    ):
+        """Init method."""
+        self._cell = np.array(cell, dtype="double")
+        self._scaled_positions = np.array(scaled_positions, dtype="double")
+        self._numbers = np.array(numbers, dtype="intc")
+
+    def __len__(self):
+        """Return number of atoms."""
+        return len(self.numbers)
+
+    @property
+    def cell(self):
+        """Setter and getter of basis vectors. For getter, copy is returned."""
+        return self._cell.copy()
+
+    @property
+    def scaled_positions(self):
+        """Setter and getter of scaled positions. For getter, copy is returned."""
+        return self._scaled_positions.copy()
+
+    @property
+    def numbers(self):
+        """Setter and getter of atomic numbers. For getter, copy is returned."""
+        return self._numbers.copy()
+
+    def totuple(self):
+        """Return (cell, scaled_position, numbers)."""
+        return (self._cell, self._scaled_positions, self._numbers)
diff --git a/src/symfc/utils/utils_O1.py b/src/symfc/utils/utils_O1.py
new file mode 100644
index 0000000..ad5bfb5
--- /dev/null
+++ b/src/symfc/utils/utils_O1.py
@@ -0,0 +1,173 @@
+"""Utility functions for 1st order force constants."""
+
+import numpy as np
+from scipy.sparse import csr_array, kron
+
+from symfc.spg_reps import SpgRepsO1
+
+from .utils import get_indep_atoms_by_lat_trans
+
+
+def get_lat_trans_decompr_indices(trans_perms: np.ndarray) -> np.ndarray:
+    """Return indices to de-compress compressed matrix by lat-trans-sym.
+
+    Usage
+    -----
+    vec[indices] of shape (n_a*3,) gives an array of shape=(N*3,).
+    1/sqrt(n_lp) must be multiplied manually after decompression to mimic
+    get_lat_trans_compr_matrix.
+
+    Parameters
+    ----------
+    trans_perms : ndarray
+        Permutation of atomic indices by lattice translational symmetry.
+        dtype='intc'.
+        shape=(n_l, N), where n_l and N are the numbers of lattce points and
+        atoms in supercell.
+
+    Returns
+    -------
+    indices : ndarray
+        Indices of n_a * 3 elements.
+        shape=(N*3,), dtype='int_'.
+
+    """
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    n_a = len(indep_atoms)
+    N = trans_perms.shape[1]
+    n_lp = N // n_a
+    size_row = N * 3
+
+    n = 0
+    indices = np.zeros(size_row, dtype="int_")
+    for i_patom in indep_atoms:
+        index_shift = trans_perms[:, i_patom] * 3
+        for a in range(3):
+            indices[index_shift + a] = n
+            n += 1
+    assert n * n_lp == size_row
+    return indices
+
+
+def get_lat_trans_compr_indices(trans_perms: np.ndarray) -> np.ndarray:
+    """Return indices to compress matrix by lat-trans-sym.
+
+    Usage
+    -----
+    vec[indices] of shape (N*3,) vec gives an array of shape=(n_a*3, n_lp).
+    1/sqrt(n_lp) must be multiplied manually after compression to mimic
+    get_lat_trans_compr_matrix.
+
+    Parameters
+    ----------
+    trans_perms : ndarray
+        Permutation of atomic indices by lattice translational symmetry.
+        dtype='intc'. shape=(n_l, N), where n_l and N are the numbers of lattce
+        points and atoms in supercell.
+
+    Returns
+    -------
+    indices : ndarray
+        shape=(n_a*N9, n_lp), dtype='int_'.
+
+    """
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    n_a = len(indep_atoms)
+    N = trans_perms.shape[1]
+    n_lp = N // n_a
+    size_row = N * 3
+
+    n = 0
+    indices = np.zeros((n_a * 3, n_lp), dtype="int_")
+    for i_patom in indep_atoms:
+        for a in range(3):
+            indices[n, :] = trans_perms[:, i_patom] * 3 + a
+            n += 1
+    assert n * n_lp == size_row
+    return indices
+
+
+def get_lat_trans_compr_matrix(decompr_idx: np.ndarray, N: int, n_lp: int) -> csr_array:
+    """Return compression matrix by lattice translation symmetry.
+
+    `decompr_idx` is obtained by `get_lat_trans_decompr_indices`.
+
+    Matrix shape is (N3, n_a*3), where n_a is the number of independent
+    atoms by lattice translation symmetry.
+
+    Data order is (N, 3, n_a, 3) if it is in dense array.
+
+    """
+    N3 = N * 3
+    compression_mat = csr_array(
+        (
+            np.full(N3, 1 / np.sqrt(n_lp), dtype="double"),
+            (np.arange(N3, dtype=int), decompr_idx),
+        ),
+        shape=(N3, N3 // n_lp),
+        dtype="double",
+    )
+    return compression_mat
+
+
+def _get_atomic_lat_trans_decompr_indices(trans_perms: np.ndarray) -> np.ndarray:
+    """Return indices to de-compress compressed matrix by atom-lat-trans-sym.
+
+    This is atomic permutation only version of get_lat_trans_decompr_indices.
+
+    Usage
+    -----
+    vec[indices] of shape (n_a,) gives an array of shape=(N,).
+    1/sqrt(n_lp) must be multiplied manually after decompression.
+
+    Parameters
+    ----------
+    trans_perms : ndarray
+        Permutation of atomic indices by lattice translational symmetry.
+        dtype='intc'.
+        shape=(n_l, N), where n_l and N are the numbers of lattce points and
+        atoms in supercell.
+
+    Returns
+    -------
+    indices : ndarray
+        Indices of n_a * N elements.
+        shape=(N^2*,), dtype='int_'.
+
+    """
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    n_lp, N = trans_perms.shape
+    size_row = N
+
+    n = 0
+    indices = np.zeros(size_row, dtype="int_")
+    for i_patom in indep_atoms:
+        index_shift = trans_perms[:, i_patom]
+        indices[index_shift] = n
+        n += 1
+    assert n * n_lp == size_row
+    return indices
+
+
+def get_compr_coset_reps_sum(spg_reps: SpgRepsO1):
+    """Return compressed projector of coset reps sum."""
+    trans_perms = spg_reps.translation_permutations
+    n_lp, N = trans_perms.shape
+    size = N * 3 // n_lp
+    coset_reps_sum = csr_array(([], ([], [])), shape=(size, size), dtype="double")
+    atomic_decompr_idx = _get_atomic_lat_trans_decompr_indices(trans_perms)
+    C = csr_array(
+        (
+            np.ones(N, dtype=int),
+            (np.arange(N, dtype=int), atomic_decompr_idx),
+        ),
+        shape=(N, N // n_lp),
+    )
+    factor = 1 / n_lp / len(spg_reps.unique_rotation_indices)
+    for i, _ in enumerate(spg_reps.unique_rotation_indices):
+        mat = spg_reps.get_sigma1_rep(i)
+        mat = mat @ C
+        mat = C.T @ mat
+        coset_reps_sum += kron(mat, spg_reps.r_reps[i] * factor)
+
+    return coset_reps_sum
diff --git a/src/symfc/utils/utils_O2.py b/src/symfc/utils/utils_O2.py
new file mode 100644
index 0000000..431839e
--- /dev/null
+++ b/src/symfc/utils/utils_O2.py
@@ -0,0 +1,281 @@
+"""Utility functions for 2nd order force constants."""
+
+import itertools
+
+import numpy as np
+from scipy.sparse import csr_array, kron
+
+from symfc.spg_reps import SpgRepsO2
+
+from .utils import get_indep_atoms_by_lat_trans
+
+
+def get_spg_perm_projector(spg_reps: SpgRepsO2, decompr_idx: np.ndarray) -> csr_array:
+    """Compute compact spg+perm projector matrix using kron.
+
+    This computes C.T @ spg_proj @ perm_proj @ C, where C is
+    ``compression_mat``.
+
+    """
+    trans_perms = spg_reps.translation_permutations
+    n_lp, N = trans_perms.shape
+    compression_mat = get_lat_trans_compr_matrix(decompr_idx, N, n_lp)
+    coset_reps_sum = get_compr_coset_reps_sum(spg_reps)
+    C_perm = get_perm_compr_matrix(N)
+    perm = C_perm.T @ compression_mat
+    perm = perm.T @ perm
+    return coset_reps_sum @ perm
+
+
+def get_lat_trans_decompr_indices(trans_perms: np.ndarray) -> np.ndarray:
+    """Return indices to de-compress compressed matrix by lat-trans-sym.
+
+    Usage
+    -----
+    vec[indices] of shape (n_a*N*9,) gives an array of shape=(N**2*9,).
+    1/sqrt(n_lp) must be multiplied manually after decompression to mimic
+    get_lat_trans_compr_matrix.
+
+    Parameters
+    ----------
+    trans_perms : ndarray
+        Permutation of atomic indices by lattice translational symmetry.
+        dtype='intc'.
+        shape=(n_l, N), where n_l and N are the numbers of lattce points and
+        atoms in supercell.
+
+    Returns
+    -------
+    indices : ndarray
+        Indices of n_a * N9 elements.
+        shape=(N^2*9,), dtype='int_'.
+
+    """
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    n_a = len(indep_atoms)
+    N = trans_perms.shape[1]
+    n_lp = N // n_a
+    size_row = (N * 3) ** 2
+
+    n = 0
+    indices = np.zeros(size_row, dtype="int_")
+    for i_patom in indep_atoms:
+        index_shift_i = trans_perms[:, i_patom] * N * 9
+        for j in range(N):
+            index_shift = index_shift_i + trans_perms[:, j] * 9
+            for ab in range(9):
+                indices[index_shift + ab] = n
+                n += 1
+    assert n * n_lp == size_row
+    return indices
+
+
+def get_lat_trans_compr_indices(trans_perms: np.ndarray) -> np.ndarray:
+    """Return indices to compress matrix by lat-trans-sym.
+
+    Usage
+    -----
+    vec[indices] of shape (N**2*9,) vec gives an array of shape=(n_a*N*9, n_lp).
+    1/sqrt(n_lp) must be multiplied manually after compression to mimic
+    get_lat_trans_compr_matrix.
+
+    Parameters
+    ----------
+    trans_perms : ndarray
+        Permutation of atomic indices by lattice translational symmetry.
+        dtype='intc'. shape=(n_l, N), where n_l and N are the numbers of lattce
+        points and atoms in supercell.
+
+    Returns
+    -------
+    indices : ndarray
+        shape=(n_a*N9, n_lp), dtype='int_'.
+
+    """
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    n_a = len(indep_atoms)
+    N = trans_perms.shape[1]
+    n_lp = N // n_a
+    size_row = (N * 3) ** 2
+
+    n = 0
+    indices = np.zeros((n_a * N * 9, n_lp), dtype="int_")
+    for i_patom in indep_atoms:
+        for j in range(N):
+            for ab in range(9):
+                indices[n, :] = (
+                    trans_perms[:, i_patom] * 9 * N + trans_perms[:, j] * 9 + ab
+                )
+                n += 1
+    assert n * n_lp == size_row
+    return indices
+
+
+def get_lat_trans_compr_matrix(decompr_idx: np.ndarray, N: int, n_lp: int) -> csr_array:
+    """Return compression matrix by lattice translation symmetry.
+
+    `decompr_idx` is obtained by `get_lat_trans_decompr_indices`.
+
+    Matrix shape is (NN33, n_a*N33), where n_a is the number of independent
+    atoms by lattice translation symmetry.
+
+    Data order is (N, N, 3, 3, n_a, N, 3, 3) if it is in dense array.
+
+    """
+    NN9 = N**2 * 9
+    compression_mat = csr_array(
+        (
+            np.full(NN9, 1 / np.sqrt(n_lp), dtype="double"),
+            (np.arange(NN9, dtype=int), decompr_idx),
+        ),
+        shape=(NN9, NN9 // n_lp),
+        dtype="double",
+    )
+    return compression_mat
+
+
+def get_perm_compr_matrix(natom: int) -> csr_array:
+    """Return compression matrix by permutation symmetry.
+
+    Parameters
+    ----------
+    natom : int
+        Number of atoms in supercell.
+
+    Matrix shape is (NN33,(N*3)(N*3+1)/2).
+    Non-zero only ijab and jiba column elements for ijab rows.
+    Rows upper right NN33 matrix elements are selected for rows.
+
+    For the computational performance, get_perm_compr_matrix is implemented in a
+    tricky way effectively using numpy features, and so it is not easy to read.
+    What is expected may be found reading _get_perm_compr_matrix_reference.
+
+    """
+    N = natom
+    A = np.transpose(
+        np.arange(N**2 * 9).reshape(N, N, 3, 3), axes=(0, 2, 1, 3)
+    ).reshape(N * 3, N * 3)
+    ut = np.triu_indices_from(A, k=1)
+    diag = np.diagonal(A)
+    row = np.hstack((np.stack((A[ut], A.T[ut]), axis=1).ravel(), diag))
+    col = np.hstack(
+        (
+            np.repeat(np.arange(len(ut[0]), dtype=int), 2),
+            np.arange(len(ut[0]), len(ut[0]) + len(diag), dtype=int),
+        )
+    )
+    data = np.hstack((np.full(len(ut[0]) * 2, np.sqrt(2) / 2), np.full(len(diag), 1)))
+    return csr_array(
+        (data, (row, col)),
+        shape=(N**2 * 9, (N * 3 * (N * 3 + 1)) // 2),
+        dtype="double",
+    )
+
+
+def _get_atomic_lat_trans_decompr_indices(trans_perms: np.ndarray) -> np.ndarray:
+    """Return indices to de-compress compressed matrix by atom-lat-trans-sym.
+
+    This is atomic permutation only version of get_lat_trans_decompr_indices.
+
+    Usage
+    -----
+    vec[indices] of shape (n_a*N,) gives an array of shape=(N**2,).
+    1/sqrt(n_lp) must be multiplied manually after decompression.
+
+    Parameters
+    ----------
+    trans_perms : ndarray
+        Permutation of atomic indices by lattice translational symmetry.
+        dtype='intc'.
+        shape=(n_l, N), where n_l and N are the numbers of lattce points and
+        atoms in supercell.
+
+    Returns
+    -------
+    indices : ndarray
+        Indices of n_a * N elements.
+        shape=(N^2*,), dtype='int_'.
+
+    """
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    n_lp, N = trans_perms.shape
+    size_row = N**2
+
+    n = 0
+    indices = np.zeros(size_row, dtype="int_")
+    for i_patom in indep_atoms:
+        index_shift_i = trans_perms[:, i_patom] * N
+        for j in range(N):
+            index_shift = index_shift_i + trans_perms[:, j]
+            indices[index_shift] = n
+            n += 1
+    assert n * n_lp == size_row
+    return indices
+
+
+def get_compr_coset_reps_sum(spg_reps: SpgRepsO2):
+    """Return compressed projector of coset reps sum."""
+    trans_perms = spg_reps.translation_permutations
+    n_lp, N = trans_perms.shape
+    size = N**2 * 9 // n_lp
+    coset_reps_sum = csr_array(([], ([], [])), shape=(size, size), dtype="double")
+    atomic_decompr_idx = _get_atomic_lat_trans_decompr_indices(trans_perms)
+    C = csr_array(
+        (
+            np.ones(N**2, dtype=int),
+            (np.arange(N**2, dtype=int), atomic_decompr_idx),
+        ),
+        shape=(N**2, N**2 // n_lp),
+    )
+    factor = 1 / n_lp / len(spg_reps.unique_rotation_indices)
+    for i, _ in enumerate(spg_reps.unique_rotation_indices):
+        mat = spg_reps.get_sigma2_rep(i)
+        mat = mat @ C
+        mat = C.T @ mat
+        coset_reps_sum += kron(mat, spg_reps.r_reps[i] * factor)
+
+    return coset_reps_sum
+
+
+def _get_perm_compr_matrix_reference(natom: int) -> csr_array:
+    """Return compression matrix by permutation symmetry.
+
+    This is a reference implementation of get_perm_compr_matrix. The order of
+    columns is difference from get_perm_compr_matrix, but it is OK if C.T@C
+    is the same.
+
+    Matrix shape is (NN33,(N*3)(N*3+1)/2). Non-zero only ijab and jiba column
+    elements for ijab rows. Rows upper right NN33 matrix elements are selected
+    for rows.
+
+    """
+
+    def to_serial(i: int, a: int, j: int, b: int, natom: int) -> int:
+        """Return NN33-1D index."""
+        return (i * 9 * natom) + (j * 9) + (a * 3) + b
+
+    col, row, data = [], [], []
+    val = np.sqrt(2) / 2
+    size_row = natom**2 * 9
+
+    n = 0
+    for ia, jb in itertools.combinations_with_replacement(range(natom * 3), 2):
+        i_i = ia // 3
+        i_a = ia % 3
+        i_j = jb // 3
+        i_b = jb % 3
+        col.append(n)
+        row.append(to_serial(i_i, i_a, i_j, i_b, natom))
+        if i_i == i_j and i_a == i_b:
+            data.append(1)
+        else:
+            data.append(val)
+            col.append(n)
+            row.append(to_serial(i_j, i_b, i_i, i_a, natom))
+            data.append(val)
+        n += 1
+    if (natom * 3) % 2 == 1:
+        assert (natom * 3) * ((natom * 3 + 1) // 2) == n, f"{natom}, {n}"
+    else:
+        assert ((natom * 3) // 2) * (natom * 3 + 1) == n, f"{natom}, {n}"
+    return csr_array((data, (row, col)), shape=(size_row, n), dtype="double")
diff --git a/src/symfc/utils/utils_O3.py b/src/symfc/utils/utils_O3.py
new file mode 100644
index 0000000..ef6041f
--- /dev/null
+++ b/src/symfc/utils/utils_O3.py
@@ -0,0 +1,154 @@
+"""Utility functions for 3rd order force constants."""
+
+import numpy as np
+from scipy.sparse import csr_array, kron
+
+from symfc.spg_reps import SpgRepsO3
+from symfc.utils.utils import get_indep_atoms_by_lat_trans
+
+
+def get_atomic_lat_trans_decompr_indices_O3(trans_perms: np.ndarray) -> np.ndarray:
+    """Return indices to de-compress compressed matrix by atom-lat-trans-sym.
+
+    This is atomic permutation only version of get_lat_trans_decompr_indices.
+
+    Usage
+    -----
+    vec[indices] of shape (n_a*N*N,) gives an array of shape=(N**3,).
+    1/sqrt(n_lp) must be multiplied manually after decompression.
+
+    Parameters
+    ----------
+    trans_perms : ndarray
+        Permutation of atomic indices by lattice translational symmetry.
+        dtype='intc'.
+        shape=(n_l, N), where n_l and N are the numbers of lattce points and
+        atoms in supercell.
+
+    Returns
+    -------
+    indices : ndarray
+        Indices of n_a * N * N elements.
+        shape=(N**3,), dtype='int_'.
+
+    """
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    n_lp, N = trans_perms.shape
+    size_row = N**3
+
+    n = 0
+    indices = np.zeros(size_row, dtype="int_")
+    for i_patom in indep_atoms:
+        index_shift_i = trans_perms[:, i_patom] * N**2
+        for j in range(N):
+            index_shift_j = index_shift_i + trans_perms[:, j] * N
+            for k in range(N):
+                index_shift = index_shift_j + trans_perms[:, k]
+                indices[index_shift] = n
+                n += 1
+    assert n * n_lp == size_row
+    return indices
+
+
+def get_lat_trans_compr_matrix_O3(trans_perms):
+    """Return lat trans compression matrix."""
+    n_lp, N = trans_perms.shape
+    decompr_idx = get_lat_trans_decompr_indices_O3(trans_perms)
+    c_trans = _get_lat_trans_compr_matrix_O3(decompr_idx, N, n_lp)
+    return c_trans
+
+
+def get_lat_trans_decompr_indices_O3(trans_perms: np.ndarray) -> np.ndarray:
+    """Return indices to de-compress compressed matrix by lat-trans-sym.
+
+    Usage
+    -----
+    vec[indices] of shape (n_a*N*N*27,) gives an array of shape=(N**3*27,).
+    1/sqrt(n_lp) must be multiplied manually after decompression to mimic
+    get_lat_trans_compr_matrix.
+
+    Parameters
+    ----------
+    trans_perms : ndarray
+        Permutation of atomic indices by lattice translational symmetry.
+        dtype='intc'.
+        shape=(n_l, N), where n_l and N are the numbers of lattce points and
+        atoms in supercell.
+
+    Returns
+    -------
+    indices : ndarray
+        Indices of n_a * N * N * 27 elements.
+        shape=(N^3*27,), dtype='int_'.
+
+    """
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    n_a = len(indep_atoms)
+    N = trans_perms.shape[1]
+    n_lp = N // n_a
+    size_row = 27 * N**3
+
+    trans_perms = trans_perms.astype("int_")
+    n = 0
+    indices = np.zeros(size_row, dtype="int_")
+    for i_patom in indep_atoms:
+        index_shift_i = trans_perms[:, i_patom] * N**2 * 27
+        for j in range(N):
+            index_shift_j = index_shift_i + trans_perms[:, j] * N * 27
+            for k in range(N):
+                index_shift = index_shift_j + trans_perms[:, k] * 27
+                for ab in range(27):
+                    indices[index_shift + ab] = n
+                    n += 1
+    assert n * n_lp == size_row
+    return indices
+
+
+def get_compr_coset_reps_sum_O3(spg_reps: SpgRepsO3) -> csr_array:
+    """Return compr matrix of sum of coset reps."""
+    trans_perms = spg_reps.translation_permutations
+    n_lp, N = trans_perms.shape
+    size = N**3 * 27 // n_lp
+    coset_reps_sum = csr_array(([], ([], [])), shape=(size, size), dtype="double")
+    atomic_decompr_idx = get_atomic_lat_trans_decompr_indices_O3(trans_perms)
+    C = csr_array(
+        (
+            np.ones(N**3, dtype=int),
+            (np.arange(N**3, dtype=int), atomic_decompr_idx),
+        ),
+        shape=(N**3, N**3 // n_lp),
+    )
+    factor = 1 / n_lp / len(spg_reps.unique_rotation_indices)
+    for i, _ in enumerate(spg_reps.unique_rotation_indices):
+        mat = spg_reps.get_sigma3_rep(i)
+        mat = mat @ C
+        mat = C.T @ mat
+        coset_reps_sum += kron(mat, spg_reps.r_reps[i] * factor)
+
+    return coset_reps_sum
+
+
+def _get_lat_trans_compr_matrix_O3(
+    decompr_idx: np.ndarray, N: int, n_lp: int
+) -> csr_array:
+    """Return compression matrix by lattice translation symmetry.
+
+    `decompr_idx` is obtained by `get_lat_trans_decompr_indices`.
+
+    Matrix shape is (NNN333, n_a*NN333), where n_a is the number of independent
+    atoms by lattice translation symmetry.
+
+    Data order is (N, N, N, 3, 3, 3, n_a, N, N, 3, 3, 3)
+    if it is in dense array.
+
+    """
+    NNN27 = N**3 * 27
+    compression_mat = csr_array(
+        (
+            np.full(NNN27, 1 / np.sqrt(n_lp), dtype="double"),
+            (np.arange(NNN27, dtype=int), decompr_idx),
+        ),
+        shape=(NNN27, NNN27 // n_lp),
+        dtype="double",
+    )
+    return compression_mat
diff --git a/src/symfc/version.py b/src/symfc/version.py
new file mode 100644
index 0000000..5be14f4
--- /dev/null
+++ b/src/symfc/version.py
@@ -0,0 +1,3 @@
+"""Version number of the package."""
+
+__version__ = "1.0.0"


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..411c660
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,41 @@
+name: symfc test using conda-forge environment
+
+on:
+  pull_request:
+    branches: [ develop ]
+
+jobs:
+  build-linux:
+    runs-on: ubuntu-latest
+    defaults:
+      run:
+        shell: bash -l {0}
+    strategy:
+      matrix:
+        python-version: ["3.11"]
+
+    steps:
+    - uses: actions/checkout@v4
+    - uses: conda-incubator/setup-miniconda@v3
+      with:
+        auto-update-conda: true
+        channels: conda-forge
+        channel-priority: strict
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        conda activate test
+        conda install --yes -c conda-forge python=${{ matrix.python-version }}
+        conda install --yes -c conda-forge spglib scipy numpy pytest codecov pytest-cov
+    - name: Setup symfc
+      run: |
+        conda activate test
+        pip install -e . -vvv
+    - name: Test with pytest
+      run: |
+        conda activate test
+        pytest -v --cov=./ --cov-report=xml
+    - name: Upload coverage reports to Codecov
+      uses: codecov/codecov-action@v3
+      with:
+        verbose: true
diff --git a/tests/NaCl-unitcell.yaml b/tests/NaCl-unitcell.yaml
new file mode 100644
index 0000000..e165aeb
--- /dev/null
+++ b/tests/NaCl-unitcell.yaml
@@ -0,0 +1,29 @@
+lattice:
+- [     5.690301476175671,     0.000000000000000,     0.000000000000000 ] # a
+- [     0.000000000000000,     5.690301476175671,     0.000000000000000 ] # b
+- [     0.000000000000000,     0.000000000000000,     5.690301476175671 ] # c
+points:
+- symbol: Na # 1
+  coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+  mass: 22.989769
+- symbol: Na # 2
+  coordinates: [  0.000000000000000,  0.500000000000000,  0.500000000000000 ]
+  mass: 22.989769
+- symbol: Na # 3
+  coordinates: [  0.500000000000000,  0.000000000000000,  0.500000000000000 ]
+  mass: 22.989769
+- symbol: Na # 4
+  coordinates: [  0.500000000000000,  0.500000000000000,  0.000000000000000 ]
+  mass: 22.989769
+- symbol: Cl # 5
+  coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
+  mass: 35.453000
+- symbol: Cl # 6
+  coordinates: [  0.500000000000000,  0.000000000000000,  0.000000000000000 ]
+  mass: 35.453000
+- symbol: Cl # 7
+  coordinates: [  0.000000000000000,  0.500000000000000,  0.000000000000000 ]
+  mass: 35.453000
+- symbol: Cl # 8
+  coordinates: [  0.000000000000000,  0.000000000000000,  0.500000000000000 ]
+  mass: 35.453000
diff --git a/tests/basis_sets/test_basis_sets_O2.py b/tests/basis_sets/test_basis_sets_O2.py
new file mode 100644
index 0000000..6caa52f
--- /dev/null
+++ b/tests/basis_sets/test_basis_sets_O2.py
@@ -0,0 +1,92 @@
+"""Tests of FCBasisSetO2."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from symfc.basis_sets import FCBasisSetO2
+from symfc.solvers.solver_O2 import FCSolverO2
+from symfc.utils.utils import SymfcAtoms
+
+cwd = Path(__file__).parent
+
+
+def test_fc_basis_set_o2():
+    """Test symmetry adapted basis sets of FCBasisSetO2."""
+    lattice = np.array([[2, 0, 0], [0, 2, 0], [0, 0, 2]])
+    positions = np.array([[0, 0, 0], [0.5, 0.5, 0.5]])
+    numbers = [1, 1]
+    supercell = SymfcAtoms(cell=lattice, scaled_positions=positions, numbers=numbers)
+    sbs = FCBasisSetO2(supercell, log_level=1).run()
+
+    np.testing.assert_allclose(
+        np.sort(sbs.basis_set), [[-np.sqrt(2) / 2], [np.sqrt(2) / 2]]
+    )
+
+    comp_mat = sbs.compression_matrix
+    np.testing.assert_allclose(comp_mat.data, [0.40824829046386313] * comp_mat.size)
+    ref_col = [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0]
+    ref_row = [0, 4, 8, 9, 13, 17, 18, 22, 26, 27, 31, 35]
+    np.testing.assert_array_equal(comp_mat.tocoo().col, ref_col)
+    np.testing.assert_array_equal(comp_mat.tocoo().row, ref_row)
+
+    compact_comp_mat = sbs.compact_compression_matrix
+    np.testing.assert_allclose(
+        compact_comp_mat.data, [0.40824829046386313] * compact_comp_mat.size
+    )
+    ref_col = [0, 0, 0, 1, 1, 1]
+    ref_row = [0, 4, 8, 9, 13, 17]
+    np.testing.assert_array_equal(compact_comp_mat.tocoo().col, ref_col)
+    np.testing.assert_array_equal(compact_comp_mat.tocoo().row, ref_row)
+
+    assert np.linalg.norm(sbs.basis_set) == pytest.approx(1.0)
+
+
+@pytest.mark.parametrize("is_compact_fc", [True, False])
+def test_fc2_NaCl_222(
+    ph_nacl_222: tuple[SymfcAtoms, np.ndarray, np.ndarray], is_compact_fc: bool
+):
+    """Test force constants by NaCl 64 atoms supercell."""
+    _assert_fc(ph_nacl_222, "NaCl_222", is_compact_fc)
+
+
+@pytest.mark.parametrize("is_compact_fc", [True, False])
+def test_fc2_SnO2_223(
+    ph_sno2_223: tuple[SymfcAtoms, np.ndarray, np.ndarray], is_compact_fc: bool
+):
+    """Test force constants by SnO2 72 atoms supercell."""
+    _assert_fc(ph_sno2_223, "SnO2_223", is_compact_fc)
+
+
+def test_fc2_SiO2_221(ph_sio2_221: tuple[SymfcAtoms, np.ndarray, np.ndarray]):
+    """Test force constants by SiO2 36 atoms supercell."""
+    _assert_fc(ph_sio2_221, "SiO2_221")
+
+
+def test_fc2_GaN_442(ph_gan_442: tuple[SymfcAtoms, np.ndarray, np.ndarray]):
+    """Test force constants by GaN 128 atoms supercell."""
+    _assert_fc(ph_gan_442, "GaN_442")
+
+
+def test_fc2_GaN_222(ph_gan_222: tuple[SymfcAtoms, np.ndarray, np.ndarray]):
+    """Test force constants by GaN 32 atoms supercell."""
+    _assert_fc(ph_gan_222, "GaN_222")
+
+
+def _assert_fc(
+    ph: tuple[SymfcAtoms, np.ndarray, np.ndarray], name: str, is_compact_fc: bool = True
+):
+    supercell, displacements, forces = ph
+    basis_set = FCBasisSetO2(supercell, log_level=1).run()
+    fc_solver = FCSolverO2(basis_set, log_level=1).solve(displacements, forces)
+    if is_compact_fc:
+        fc = fc_solver.compact_fc
+        # np.savetxt(f"compact_fc_{name}.xz", fc.ravel())
+        fc_ref = np.loadtxt(cwd / ".." / f"compact_fc_{name}.xz").reshape(fc.shape)
+    else:
+        fc = fc_solver.full_fc
+        fc_ref = np.loadtxt(cwd / ".." / f"full_fc_{name}.xz").reshape(fc.shape)
+    np.testing.assert_allclose(fc, fc_ref, atol=1e-6)
diff --git a/tests/basis_sets/test_basis_sets_O3.py b/tests/basis_sets/test_basis_sets_O3.py
new file mode 100644
index 0000000..d44e3d1
--- /dev/null
+++ b/tests/basis_sets/test_basis_sets_O3.py
@@ -0,0 +1,100 @@
+"""Tests of FCBasisSetO3."""
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from symfc.basis_sets import FCBasisSetO2, FCBasisSetO3
+from symfc.solvers import FCSolverO2O3
+from symfc.utils.utils import SymfcAtoms
+from symfc.utils.utils_O3 import get_lat_trans_compr_matrix_O3
+
+cwd = Path(__file__).parent
+
+
+def test_fc_basis_set_o3():
+    """Test symmetry adapted basis sets of FCBasisSetO3."""
+    a = 5.4335600299999998
+    lattice = np.array([[a, 0, 0], [0, a, 0], [0, 0, a]])
+    positions = np.array(
+        [
+            [0.875, 0.875, 0.875],
+            [0.875, 0.375, 0.375],
+            [0.375, 0.875, 0.375],
+            [0.375, 0.375, 0.875],
+            [0.125, 0.125, 0.125],
+            [0.125, 0.625, 0.625],
+            [0.625, 0.125, 0.625],
+            [0.625, 0.625, 0.125],
+        ]
+    )
+    numbers = [1, 1, 1, 1, 1, 1, 1, 1]
+    supercell = SymfcAtoms(cell=lattice, scaled_positions=positions, numbers=numbers)
+    sbs = FCBasisSetO3(supercell, log_level=1).run()
+
+    basis_ref = [
+        0.671875,
+        0.656250,
+        0.765625,
+        1.000000,
+        0.750000,
+        0.875000,
+        0.562500,
+        0.765625,
+        0.718750,
+        0.671875,
+        1.000000,
+        0.875000,
+        0.875000,
+        0.625000,
+        0.750000,
+        0.562500,
+        0.875000,
+    ]
+    np.testing.assert_allclose(
+        np.sort([v @ v for v in sbs.basis_set]), np.sort(basis_ref)
+    )
+
+    assert np.linalg.norm(sbs.basis_set) ** 2 == pytest.approx(13.0)
+
+    np.testing.assert_array_equal(
+        sbs.compact_compression_matrix.tocoo().row[[0, 6, 9]], [5, 32, 42]
+    )
+    np.testing.assert_allclose(
+        sbs.compression_matrix.data[[0, 6, 9]],
+        [-0.1443375672974064, 0.058925565098878946, 0.0833333333333333],
+    )
+
+    np.testing.assert_array_equal(
+        sbs.compression_matrix.tocoo().row[[0, 6, 9]], [5, 32, 42]
+    )
+    np.testing.assert_allclose(
+        sbs.compression_matrix.data[[0, 6, 9]],
+        [-0.1443375672974064, 0.058925565098878946, 0.0833333333333333],
+    )
+
+    lat_trans_compr_matrix_O3 = get_lat_trans_compr_matrix_O3(
+        sbs.translation_permutations
+    )
+    np.testing.assert_allclose(lat_trans_compr_matrix_O3.data, [0.5] * 13824)
+    assert lat_trans_compr_matrix_O3.indices[-1] == 2726
+
+
+def test_si_111_fc3(ph3_si_111_fc3: tuple[SymfcAtoms, np.ndarray, np.ndarray]):
+    """Test fc2 and fc3 by Si-111-222 supercells and compared with ALM.
+
+    This test with ALM is skipped when ALM is not installed.
+
+    """
+    supercell, displacements, forces = ph3_si_111_fc3
+    basis_set_o2 = FCBasisSetO2(supercell, log_level=1).run()
+    basis_set_o3 = FCBasisSetO3(supercell, log_level=1).run()
+    fc_solver = FCSolverO2O3([basis_set_o2, basis_set_o3], log_level=1).solve(
+        displacements, forces
+    )
+    fc2, fc3 = fc_solver.compact_fc
+    fc2_ref = np.loadtxt(cwd / ".." / "compact_fc_Si_111_fc3_2.xz").reshape(fc2.shape)
+    fc3_ref = np.loadtxt(cwd / ".." / "compact_fc_Si_111_fc3_3.xz").reshape(fc3.shape)
+    np.testing.assert_allclose(fc2_ref, fc2, atol=1e-6)
+    np.testing.assert_allclose(fc3_ref, fc3, atol=1e-6)
diff --git a/tests/basis_sets/test_basis_sets_base.py b/tests/basis_sets/test_basis_sets_base.py
new file mode 100644
index 0000000..51a0e1d
--- /dev/null
+++ b/tests/basis_sets/test_basis_sets_base.py
@@ -0,0 +1,15 @@
+"""Tests of FCBasisSetBase."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from symfc.basis_sets import FCBasisSetBase
+from symfc.utils.utils import SymfcAtoms
+
+
+def test_base_fc_basis_set(ph_nacl_222: tuple[SymfcAtoms, np.ndarray, np.ndarray]):
+    """Test that FCBasisSet can not be instantiate."""
+    with pytest.raises(TypeError):
+        _ = FCBasisSetBase(ph_nacl_222[0])
diff --git a/tests/compact_fc_GaN_222.xz b/tests/compact_fc_GaN_222.xz
new file mode 100644
index 0000000..dd5f7ca
Binary files /dev/null and b/tests/compact_fc_GaN_222.xz differ
diff --git a/tests/compact_fc_GaN_442.xz b/tests/compact_fc_GaN_442.xz
new file mode 100644
index 0000000..1b02bf9
Binary files /dev/null and b/tests/compact_fc_GaN_442.xz differ
diff --git a/tests/compact_fc_NaCl_222.xz b/tests/compact_fc_NaCl_222.xz
new file mode 100644
index 0000000..5c3e3c4
Binary files /dev/null and b/tests/compact_fc_NaCl_222.xz differ
diff --git a/tests/compact_fc_SiO2_221.xz b/tests/compact_fc_SiO2_221.xz
new file mode 100644
index 0000000..cc0efe7
Binary files /dev/null and b/tests/compact_fc_SiO2_221.xz differ
diff --git a/tests/compact_fc_Si_111_fc3_2.xz b/tests/compact_fc_Si_111_fc3_2.xz
new file mode 100644
index 0000000..2420a87
Binary files /dev/null and b/tests/compact_fc_Si_111_fc3_2.xz differ
diff --git a/tests/compact_fc_Si_111_fc3_3.xz b/tests/compact_fc_Si_111_fc3_3.xz
new file mode 100644
index 0000000..1e14cef
Binary files /dev/null and b/tests/compact_fc_Si_111_fc3_3.xz differ
diff --git a/tests/compact_fc_SnO2_223.xz b/tests/compact_fc_SnO2_223.xz
new file mode 100644
index 0000000..9119f29
Binary files /dev/null and b/tests/compact_fc_SnO2_223.xz differ
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..cd49528
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,877 @@
+"""Pytest conftest.py."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from symfc.utils.utils import SymfcAtoms
+
+cwd = Path(__file__).parent
+
+
+def pytest_addoption(parser):
+    """Add command option to pytest."""
+    parser.addoption(
+        "--runbig", action="store_true", default=False, help="run big tests"
+    )
+
+
+def pytest_configure(config):
+    """Set up marker big."""
+    config.addinivalue_line("markers", "big: mark test as big to run")
+
+
+def pytest_collection_modifyitems(config, items):
+    """Add mechanism to run with --runbig."""
+    if config.getoption("--runbig"):
+        # --runbig given in cli: do not skip slow tests
+        return
+    skip_big = pytest.mark.skip(reason="need --runbig option to run")
+    for item in items:
+        if "big" in item.keywords:
+            item.add_marker(skip_big)
+
+
+@pytest.fixture(scope="session")
+def cell_nacl_111() -> SymfcAtoms:
+    """Return unitcell of NaCl."""
+    lattice = [
+        [5.690301476175671, 0.000000000000000, 0.000000000000000],
+        [0.000000000000000, 5.690301476175671, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 5.690301476175671],
+    ]
+    points = [
+        [0.000000000000000, 0.000000000000000, 0.000000000000000],
+        [0.000000000000000, 0.500000000000000, 0.500000000000000],
+        [0.500000000000000, 0.000000000000000, 0.500000000000000],
+        [0.500000000000000, 0.500000000000000, 0.000000000000000],
+        [0.500000000000000, 0.500000000000000, 0.500000000000000],
+        [0.500000000000000, 0.000000000000000, 0.000000000000000],
+        [0.000000000000000, 0.500000000000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 0.500000000000000],
+    ]
+    numbers = [11, 11, 11, 11, 17, 17, 17, 17]
+    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)
+    return cell
+
+
+@pytest.fixture(scope="session")
+def cell_gan_111() -> SymfcAtoms:
+    """Return unitcell of GaN."""
+    lattice = [
+        [3.180765520000000, 0.000000000000000, 0.000000000000000],
+        [-1.590382760000000, 2.754623740000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 5.182605840000000],
+    ]
+    points = [
+        [0.333333329999992, 0.666666669999983, 0.124191919999999],
+        [0.666666669999991, 0.333333329999998, 0.624191919999998],
+        [0.333333329999992, 0.666666669999983, 0.500808080000002],
+        [0.666666669999991, 0.333333329999998, 0.000808080000000],
+    ]
+    numbers = [7, 7, 31, 31]
+    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)
+    return cell
+
+
+@pytest.fixture(scope="session")
+def ph_nacl_222() -> tuple[SymfcAtoms, np.ndarray, np.ndarray]:
+    """Return NaCl-222 data."""
+    lattice = [
+        [11.281120000000000, 0.000000000000000, 0.000000000000000],
+        [0.000000000000000, 11.281120000000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 11.281120000000000],
+    ]
+    points = [
+        [0.000000000000000, 0.000000000000000, 0.000000000000000],
+        [0.500000000000000, 0.000000000000000, 0.000000000000000],
+        [0.000000000000000, 0.500000000000000, 0.000000000000000],
+        [0.500000000000000, 0.500000000000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 0.500000000000000],
+        [0.500000000000000, 0.000000000000000, 0.500000000000000],
+        [0.000000000000000, 0.500000000000000, 0.500000000000000],
+        [0.500000000000000, 0.500000000000000, 0.500000000000000],
+        [0.000000000000000, 0.250000000000000, 0.250000000000000],
+        [0.500000000000000, 0.250000000000000, 0.250000000000000],
+        [0.000000000000000, 0.750000000000000, 0.250000000000000],
+        [0.500000000000000, 0.750000000000000, 0.250000000000000],
+        [0.000000000000000, 0.250000000000000, 0.750000000000000],
+        [0.500000000000000, 0.250000000000000, 0.750000000000000],
+        [0.000000000000000, 0.750000000000000, 0.750000000000000],
+        [0.500000000000000, 0.750000000000000, 0.750000000000000],
+        [0.250000000000000, 0.000000000000000, 0.250000000000000],
+        [0.750000000000000, 0.000000000000000, 0.250000000000000],
+        [0.250000000000000, 0.500000000000000, 0.250000000000000],
+        [0.750000000000000, 0.500000000000000, 0.250000000000000],
+        [0.250000000000000, 0.000000000000000, 0.750000000000000],
+        [0.750000000000000, 0.000000000000000, 0.750000000000000],
+        [0.250000000000000, 0.500000000000000, 0.750000000000000],
+        [0.750000000000000, 0.500000000000000, 0.750000000000000],
+        [0.250000000000000, 0.250000000000000, 0.000000000000000],
+        [0.750000000000000, 0.250000000000000, 0.000000000000000],
+        [0.250000000000000, 0.750000000000000, 0.000000000000000],
+        [0.750000000000000, 0.750000000000000, 0.000000000000000],
+        [0.250000000000000, 0.250000000000000, 0.500000000000000],
+        [0.750000000000000, 0.250000000000000, 0.500000000000000],
+        [0.250000000000000, 0.750000000000000, 0.500000000000000],
+        [0.750000000000000, 0.750000000000000, 0.500000000000000],
+        [0.250000000000000, 0.250000000000000, 0.250000000000000],
+        [0.750000000000000, 0.250000000000000, 0.250000000000000],
+        [0.250000000000000, 0.750000000000000, 0.250000000000000],
+        [0.750000000000000, 0.750000000000000, 0.250000000000000],
+        [0.250000000000000, 0.250000000000000, 0.750000000000000],
+        [0.750000000000000, 0.250000000000000, 0.750000000000000],
+        [0.250000000000000, 0.750000000000000, 0.750000000000000],
+        [0.750000000000000, 0.750000000000000, 0.750000000000000],
+        [0.250000000000000, 0.000000000000000, 0.000000000000000],
+        [0.750000000000000, 0.000000000000000, 0.000000000000000],
+        [0.250000000000000, 0.500000000000000, 0.000000000000000],
+        [0.750000000000000, 0.500000000000000, 0.000000000000000],
+        [0.250000000000000, 0.000000000000000, 0.500000000000000],
+        [0.750000000000000, 0.000000000000000, 0.500000000000000],
+        [0.250000000000000, 0.500000000000000, 0.500000000000000],
+        [0.750000000000000, 0.500000000000000, 0.500000000000000],
+        [0.000000000000000, 0.250000000000000, 0.000000000000000],
+        [0.500000000000000, 0.250000000000000, 0.000000000000000],
+        [0.000000000000000, 0.750000000000000, 0.000000000000000],
+        [0.500000000000000, 0.750000000000000, 0.000000000000000],
+        [0.000000000000000, 0.250000000000000, 0.500000000000000],
+        [0.500000000000000, 0.250000000000000, 0.500000000000000],
+        [0.000000000000000, 0.750000000000000, 0.500000000000000],
+        [0.500000000000000, 0.750000000000000, 0.500000000000000],
+        [0.000000000000000, 0.000000000000000, 0.250000000000000],
+        [0.500000000000000, 0.000000000000000, 0.250000000000000],
+        [0.000000000000000, 0.500000000000000, 0.250000000000000],
+        [0.500000000000000, 0.500000000000000, 0.250000000000000],
+        [0.000000000000000, 0.000000000000000, 0.750000000000000],
+        [0.500000000000000, 0.000000000000000, 0.750000000000000],
+        [0.000000000000000, 0.500000000000000, 0.750000000000000],
+        [0.500000000000000, 0.500000000000000, 0.750000000000000],
+    ]
+    numbers = [
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        11,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+        17,
+    ]
+    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)
+    N = 20
+    dfset = np.loadtxt(cwd / "dfset_NaCl_222_rd.xz")
+    d = dfset[:, :3].reshape(N, -1, 3)
+    f = dfset[:, 3:].reshape(N, -1, 3)
+    return cell, d, f
+
+
+@pytest.fixture(scope="session")
+def ph_sno2_223() -> tuple[SymfcAtoms, np.ndarray, np.ndarray]:
+    """Return SnO2-223 data."""
+    lattice = [
+        [9.545162500000000, 0.000000000000000, 0.000000000000000],
+        [0.000000000000000, 9.545162500000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 9.649138170000001],
+    ]
+    points = [
+        [0.096865350000001, 0.403134649999997, 0.166666666666667],
+        [0.596865350000001, 0.403134649999997, 0.166666666666667],
+        [0.096865350000001, 0.903134649999997, 0.166666666666667],
+        [0.596865350000001, 0.903134649999997, 0.166666666666667],
+        [0.096865350000001, 0.403134649999997, 0.500000000000000],
+        [0.596865350000001, 0.403134649999997, 0.500000000000000],
+        [0.096865350000001, 0.903134649999997, 0.500000000000000],
+        [0.596865350000001, 0.903134649999997, 0.500000000000000],
+        [0.096865350000001, 0.403134649999997, 0.833333333333333],
+        [0.596865350000001, 0.403134649999997, 0.833333333333333],
+        [0.096865350000001, 0.903134649999997, 0.833333333333333],
+        [0.596865350000001, 0.903134649999997, 0.833333333333333],
+        [0.403134649999997, 0.096865350000001, 0.166666666666667],
+        [0.903134649999997, 0.096865350000001, 0.166666666666667],
+        [0.403134649999997, 0.596865350000001, 0.166666666666667],
+        [0.903134649999997, 0.596865350000001, 0.166666666666667],
+        [0.403134649999997, 0.096865350000001, 0.500000000000000],
+        [0.903134649999997, 0.096865350000001, 0.500000000000000],
+        [0.403134649999997, 0.596865350000001, 0.500000000000000],
+        [0.903134649999997, 0.596865350000001, 0.500000000000000],
+        [0.403134649999997, 0.096865350000001, 0.833333333333333],
+        [0.903134649999997, 0.096865350000001, 0.833333333333333],
+        [0.403134649999997, 0.596865350000001, 0.833333333333333],
+        [0.903134649999997, 0.596865350000001, 0.833333333333333],
+        [0.153134649999998, 0.153134649999998, 0.000000000000000],
+        [0.653134649999997, 0.153134649999998, 0.000000000000000],
+        [0.153134649999998, 0.653134649999997, 0.000000000000000],
+        [0.653134649999997, 0.653134649999997, 0.000000000000000],
+        [0.153134649999998, 0.153134649999998, 0.333333333333333],
+        [0.653134649999997, 0.153134649999998, 0.333333333333333],
+        [0.153134649999998, 0.653134649999997, 0.333333333333333],
+        [0.653134649999997, 0.653134649999997, 0.333333333333333],
+        [0.153134649999998, 0.153134649999998, 0.666666666666667],
+        [0.653134649999997, 0.153134649999998, 0.666666666666667],
+        [0.153134649999998, 0.653134649999997, 0.666666666666667],
+        [0.653134649999997, 0.653134649999997, 0.666666666666667],
+        [0.346865350000003, 0.346865350000003, 0.000000000000000],
+        [0.846865350000003, 0.346865350000003, 0.000000000000000],
+        [0.346865350000003, 0.846865350000003, 0.000000000000000],
+        [0.846865350000003, 0.846865350000003, 0.000000000000000],
+        [0.346865350000003, 0.346865350000003, 0.333333333333333],
+        [0.846865350000003, 0.346865350000003, 0.333333333333333],
+        [0.346865350000003, 0.846865350000003, 0.333333333333333],
+        [0.846865350000003, 0.846865350000003, 0.333333333333333],
+        [0.346865350000003, 0.346865350000003, 0.666666666666667],
+        [0.846865350000003, 0.346865350000003, 0.666666666666667],
+        [0.346865350000003, 0.846865350000003, 0.666666666666667],
+        [0.846865350000003, 0.846865350000003, 0.666666666666667],
+        [0.250000000000000, 0.250000000000000, 0.166666666666667],
+        [0.750000000000000, 0.250000000000000, 0.166666666666667],
+        [0.250000000000000, 0.750000000000000, 0.166666666666667],
+        [0.750000000000000, 0.750000000000000, 0.166666666666667],
+        [0.250000000000000, 0.250000000000000, 0.500000000000000],
+        [0.750000000000000, 0.250000000000000, 0.500000000000000],
+        [0.250000000000000, 0.750000000000000, 0.500000000000000],
+        [0.750000000000000, 0.750000000000000, 0.500000000000000],
+        [0.250000000000000, 0.250000000000000, 0.833333333333333],
+        [0.750000000000000, 0.250000000000000, 0.833333333333333],
+        [0.250000000000000, 0.750000000000000, 0.833333333333333],
+        [0.750000000000000, 0.750000000000000, 0.833333333333333],
+        [0.000000000000000, 0.000000000000000, 0.000000000000000],
+        [0.500000000000000, 0.000000000000000, 0.000000000000000],
+        [0.000000000000000, 0.500000000000000, 0.000000000000000],
+        [0.500000000000000, 0.500000000000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 0.333333333333333],
+        [0.500000000000000, 0.000000000000000, 0.333333333333333],
+        [0.000000000000000, 0.500000000000000, 0.333333333333333],
+        [0.500000000000000, 0.500000000000000, 0.333333333333333],
+        [0.000000000000000, 0.000000000000000, 0.666666666666667],
+        [0.500000000000000, 0.000000000000000, 0.666666666666667],
+        [0.000000000000000, 0.500000000000000, 0.666666666666667],
+        [0.500000000000000, 0.500000000000000, 0.666666666666667],
+    ]
+    numbers = [
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+        50,
+    ]
+    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)
+    N = 40
+    dfset = np.loadtxt(cwd / "dfset_SnO2_223_rd.xz")
+    d = dfset[:, :3].reshape(N, -1, 3)
+    f = dfset[:, 3:].reshape(N, -1, 3)
+    return cell, d, f
+
+
+@pytest.fixture(scope="session")
+def ph_sio2_221() -> tuple[SymfcAtoms, np.ndarray, np.ndarray]:
+    """Return SiO2-221 data."""
+    lattice = [
+        [9.919749440000000, 0.000000000000000, 0.000000000000000],
+        [-4.959874720000000, 8.590755020000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 5.453559970000000],
+    ]
+    points = [
+        [0.293442579999996, 0.367112664999997, 0.213280300000002],
+        [0.793442579999996, 0.367112664999997, 0.213280300000002],
+        [0.293442579999996, 0.867112664999997, 0.213280300000002],
+        [0.793442579999996, 0.867112664999997, 0.213280300000002],
+        [0.073670084999997, 0.206557419999994, 0.879946960000002],
+        [0.573670084999997, 0.206557419999994, 0.879946960000002],
+        [0.073670084999997, 0.706557419999994, 0.879946960000002],
+        [0.573670084999997, 0.706557419999994, 0.879946960000002],
+        [0.132887334999998, 0.426329914999998, 0.546613630000002],
+        [0.632887334999998, 0.426329914999998, 0.546613630000002],
+        [0.132887334999998, 0.926329914999998, 0.546613630000002],
+        [0.632887334999998, 0.926329914999998, 0.546613630000002],
+        [0.367112665000005, 0.293442580000006, 0.786719699999998],
+        [0.867112665000005, 0.293442580000006, 0.786719699999998],
+        [0.367112665000005, 0.793442580000006, 0.786719699999998],
+        [0.867112665000005, 0.793442580000006, 0.786719699999998],
+        [0.426329914999999, 0.132887335000003, 0.453386369999998],
+        [0.926329914999999, 0.132887335000003, 0.453386369999998],
+        [0.426329914999999, 0.632887335000003, 0.453386369999998],
+        [0.926329914999999, 0.632887335000003, 0.453386369999998],
+        [0.206557420000001, 0.073670085000000, 0.120053040000000],
+        [0.706557420000001, 0.073670085000000, 0.120053040000000],
+        [0.206557420000001, 0.573670085000000, 0.120053040000000],
+        [0.706557420000001, 0.573670085000000, 0.120053040000000],
+        [0.235454159999997, 0.235454160000002, 0.000000000000000],
+        [0.735454159999998, 0.235454160000002, 0.000000000000000],
+        [0.235454159999997, 0.735454160000002, 0.000000000000000],
+        [0.735454159999998, 0.735454160000002, 0.000000000000000],
+        [0.264545839999997, 0.000000000000000, 0.333333330000000],
+        [0.764545839999997, 0.000000000000000, 0.333333330000000],
+        [0.264545839999997, 0.500000000000000, 0.333333330000000],
+        [0.764545839999997, 0.500000000000000, 0.333333330000000],
+        [0.999999999999995, 0.264545839999998, 0.666666670000000],
+        [0.499999999999996, 0.264545839999998, 0.666666670000000],
+        [0.999999999999995, 0.764545839999998, 0.666666670000000],
+        [0.499999999999996, 0.764545839999998, 0.666666670000000],
+    ]
+    numbers = [
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        8,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+    ]
+    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)
+    N = 40
+    dfset = np.loadtxt(cwd / "dfset_SiO2_221_rd.xz")
+    d = dfset[:, :3].reshape(N, -1, 3)
+    f = dfset[:, 3:].reshape(N, -1, 3)
+    return cell, d, f
+
+
+@pytest.fixture(scope="session")
+def ph_gan_222() -> tuple[SymfcAtoms, np.ndarray, np.ndarray]:
+    """Return GaN-222 data."""
+    lattice = [
+        [6.361531040000000, 0.000000000000000, 0.000000000000000],
+        [-3.180765520000000, 5.509247480000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 10.365211680000000],
+    ]
+    points = [
+        [0.166666664999996, 0.333333334999992, 0.062095959999999],
+        [0.666666664999996, 0.333333334999992, 0.062095959999999],
+        [0.166666664999996, 0.833333334999992, 0.062095959999999],
+        [0.666666664999996, 0.833333334999992, 0.062095959999999],
+        [0.166666664999996, 0.333333334999992, 0.562095960000000],
+        [0.666666664999996, 0.333333334999992, 0.562095960000000],
+        [0.166666664999996, 0.833333334999992, 0.562095960000000],
+        [0.666666664999996, 0.833333334999992, 0.562095960000000],
+        [0.333333334999996, 0.166666664999999, 0.312095959999999],
+        [0.833333334999996, 0.166666664999999, 0.312095959999999],
+        [0.333333334999996, 0.666666664999999, 0.312095959999999],
+        [0.833333334999996, 0.666666664999999, 0.312095959999999],
+        [0.333333334999996, 0.166666664999999, 0.812095959999999],
+        [0.833333334999996, 0.166666664999999, 0.812095959999999],
+        [0.333333334999996, 0.666666664999999, 0.812095959999999],
+        [0.833333334999996, 0.666666664999999, 0.812095959999999],
+        [0.166666664999996, 0.333333334999992, 0.250404040000001],
+        [0.666666664999996, 0.333333334999992, 0.250404040000001],
+        [0.166666664999996, 0.833333334999992, 0.250404040000001],
+        [0.666666664999996, 0.833333334999992, 0.250404040000001],
+        [0.166666664999996, 0.333333334999992, 0.750404040000001],
+        [0.666666664999996, 0.333333334999992, 0.750404040000001],
+        [0.166666664999996, 0.833333334999992, 0.750404040000001],
+        [0.666666664999996, 0.833333334999992, 0.750404040000001],
+        [0.333333334999996, 0.166666664999999, 0.000404040000000],
+        [0.833333334999996, 0.166666664999999, 0.000404040000000],
+        [0.333333334999996, 0.666666664999999, 0.000404040000000],
+        [0.833333334999996, 0.666666664999999, 0.000404040000000],
+        [0.333333334999996, 0.166666664999999, 0.500404040000000],
+        [0.833333334999996, 0.166666664999999, 0.500404040000000],
+        [0.333333334999996, 0.666666664999999, 0.500404040000000],
+        [0.833333334999996, 0.666666664999999, 0.500404040000000],
+    ]
+    numbers = [
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+    ]
+    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)
+    N = 40
+    dfset = np.loadtxt(cwd / "dfset_GaN_222_rd.xz")
+    d = dfset[:, :3].reshape(N, -1, 3)
+    f = dfset[:, 3:].reshape(N, -1, 3)
+    return cell, d, f
+
+
+@pytest.fixture(scope="session")
+def ph_gan_442() -> tuple[SymfcAtoms, np.ndarray, np.ndarray]:
+    """Return GaN-442 data."""
+    lattice = [
+        [12.723062080000000, 0.000000000000000, 0.000000000000000],
+        [-6.361531040000000, 11.018494960000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 10.365211680000000],
+    ]
+    points = [
+        [0.083333332499998, 0.166666667499996, 0.062095959999999],
+        [0.333333332499998, 0.166666667499996, 0.062095959999999],
+        [0.583333332499998, 0.166666667499996, 0.062095959999999],
+        [0.833333332499998, 0.166666667499996, 0.062095959999999],
+        [0.083333332499998, 0.416666667499996, 0.062095959999999],
+        [0.333333332499998, 0.416666667499996, 0.062095959999999],
+        [0.583333332499998, 0.416666667499996, 0.062095959999999],
+        [0.833333332499998, 0.416666667499996, 0.062095959999999],
+        [0.083333332499998, 0.666666667499996, 0.062095959999999],
+        [0.333333332499998, 0.666666667499996, 0.062095959999999],
+        [0.583333332499998, 0.666666667499996, 0.062095959999999],
+        [0.833333332499998, 0.666666667499996, 0.062095959999999],
+        [0.083333332499998, 0.916666667499996, 0.062095959999999],
+        [0.333333332499998, 0.916666667499996, 0.062095959999999],
+        [0.583333332499998, 0.916666667499996, 0.062095959999999],
+        [0.833333332499998, 0.916666667499996, 0.062095959999999],
+        [0.083333332499998, 0.166666667499996, 0.562095960000000],
+        [0.333333332499998, 0.166666667499996, 0.562095960000000],
+        [0.583333332499998, 0.166666667499996, 0.562095960000000],
+        [0.833333332499998, 0.166666667499996, 0.562095960000000],
+        [0.083333332499998, 0.416666667499996, 0.562095960000000],
+        [0.333333332499998, 0.416666667499996, 0.562095960000000],
+        [0.583333332499998, 0.416666667499996, 0.562095960000000],
+        [0.833333332499998, 0.416666667499996, 0.562095960000000],
+        [0.083333332499998, 0.666666667499996, 0.562095960000000],
+        [0.333333332499998, 0.666666667499996, 0.562095960000000],
+        [0.583333332499998, 0.666666667499996, 0.562095960000000],
+        [0.833333332499998, 0.666666667499996, 0.562095960000000],
+        [0.083333332499998, 0.916666667499996, 0.562095960000000],
+        [0.333333332499998, 0.916666667499996, 0.562095960000000],
+        [0.583333332499998, 0.916666667499996, 0.562095960000000],
+        [0.833333332499998, 0.916666667499996, 0.562095960000000],
+        [0.166666667499998, 0.083333332499999, 0.312095959999999],
+        [0.416666667499998, 0.083333332499999, 0.312095959999999],
+        [0.666666667499998, 0.083333332499999, 0.312095959999999],
+        [0.916666667499998, 0.083333332499999, 0.312095959999999],
+        [0.166666667499998, 0.333333332500000, 0.312095959999999],
+        [0.416666667499998, 0.333333332500000, 0.312095959999999],
+        [0.666666667499998, 0.333333332500000, 0.312095959999999],
+        [0.916666667499998, 0.333333332500000, 0.312095959999999],
+        [0.166666667499998, 0.583333332500000, 0.312095959999999],
+        [0.416666667499998, 0.583333332500000, 0.312095959999999],
+        [0.666666667499998, 0.583333332500000, 0.312095959999999],
+        [0.916666667499998, 0.583333332500000, 0.312095959999999],
+        [0.166666667499998, 0.833333332500000, 0.312095959999999],
+        [0.416666667499998, 0.833333332500000, 0.312095959999999],
+        [0.666666667499998, 0.833333332500000, 0.312095959999999],
+        [0.916666667499998, 0.833333332500000, 0.312095959999999],
+        [0.166666667499998, 0.083333332499999, 0.812095959999999],
+        [0.416666667499998, 0.083333332499999, 0.812095959999999],
+        [0.666666667499998, 0.083333332499999, 0.812095959999999],
+        [0.916666667499998, 0.083333332499999, 0.812095959999999],
+        [0.166666667499998, 0.333333332500000, 0.812095959999999],
+        [0.416666667499998, 0.333333332500000, 0.812095959999999],
+        [0.666666667499998, 0.333333332500000, 0.812095959999999],
+        [0.916666667499998, 0.333333332500000, 0.812095959999999],
+        [0.166666667499998, 0.583333332500000, 0.812095959999999],
+        [0.416666667499998, 0.583333332500000, 0.812095959999999],
+        [0.666666667499998, 0.583333332500000, 0.812095959999999],
+        [0.916666667499998, 0.583333332500000, 0.812095959999999],
+        [0.166666667499998, 0.833333332500000, 0.812095959999999],
+        [0.416666667499998, 0.833333332500000, 0.812095959999999],
+        [0.666666667499998, 0.833333332500000, 0.812095959999999],
+        [0.916666667499998, 0.833333332500000, 0.812095959999999],
+        [0.083333332499998, 0.166666667499996, 0.250404040000001],
+        [0.333333332499998, 0.166666667499996, 0.250404040000001],
+        [0.583333332499998, 0.166666667499996, 0.250404040000001],
+        [0.833333332499998, 0.166666667499996, 0.250404040000001],
+        [0.083333332499998, 0.416666667499996, 0.250404040000001],
+        [0.333333332499998, 0.416666667499996, 0.250404040000001],
+        [0.583333332499998, 0.416666667499996, 0.250404040000001],
+        [0.833333332499998, 0.416666667499996, 0.250404040000001],
+        [0.083333332499998, 0.666666667499996, 0.250404040000001],
+        [0.333333332499998, 0.666666667499996, 0.250404040000001],
+        [0.583333332499998, 0.666666667499996, 0.250404040000001],
+        [0.833333332499998, 0.666666667499996, 0.250404040000001],
+        [0.083333332499998, 0.916666667499996, 0.250404040000001],
+        [0.333333332499998, 0.916666667499996, 0.250404040000001],
+        [0.583333332499998, 0.916666667499996, 0.250404040000001],
+        [0.833333332499998, 0.916666667499996, 0.250404040000001],
+        [0.083333332499998, 0.166666667499996, 0.750404040000001],
+        [0.333333332499998, 0.166666667499996, 0.750404040000001],
+        [0.583333332499998, 0.166666667499996, 0.750404040000001],
+        [0.833333332499998, 0.166666667499996, 0.750404040000001],
+        [0.083333332499998, 0.416666667499996, 0.750404040000001],
+        [0.333333332499998, 0.416666667499996, 0.750404040000001],
+        [0.583333332499998, 0.416666667499996, 0.750404040000001],
+        [0.833333332499998, 0.416666667499996, 0.750404040000001],
+        [0.083333332499998, 0.666666667499996, 0.750404040000001],
+        [0.333333332499998, 0.666666667499996, 0.750404040000001],
+        [0.583333332499998, 0.666666667499996, 0.750404040000001],
+        [0.833333332499998, 0.666666667499996, 0.750404040000001],
+        [0.083333332499998, 0.916666667499996, 0.750404040000001],
+        [0.333333332499998, 0.916666667499996, 0.750404040000001],
+        [0.583333332499998, 0.916666667499996, 0.750404040000001],
+        [0.833333332499998, 0.916666667499996, 0.750404040000001],
+        [0.166666667499998, 0.083333332499999, 0.000404040000000],
+        [0.416666667499998, 0.083333332499999, 0.000404040000000],
+        [0.666666667499998, 0.083333332499999, 0.000404040000000],
+        [0.916666667499998, 0.083333332499999, 0.000404040000000],
+        [0.166666667499998, 0.333333332500000, 0.000404040000000],
+        [0.416666667499998, 0.333333332500000, 0.000404040000000],
+        [0.666666667499998, 0.333333332500000, 0.000404040000000],
+        [0.916666667499998, 0.333333332500000, 0.000404040000000],
+        [0.166666667499998, 0.583333332500000, 0.000404040000000],
+        [0.416666667499998, 0.583333332500000, 0.000404040000000],
+        [0.666666667499998, 0.583333332500000, 0.000404040000000],
+        [0.916666667499998, 0.583333332500000, 0.000404040000000],
+        [0.166666667499998, 0.833333332500000, 0.000404040000000],
+        [0.416666667499998, 0.833333332500000, 0.000404040000000],
+        [0.666666667499998, 0.833333332500000, 0.000404040000000],
+        [0.916666667499998, 0.833333332500000, 0.000404040000000],
+        [0.166666667499998, 0.083333332499999, 0.500404040000000],
+        [0.416666667499998, 0.083333332499999, 0.500404040000000],
+        [0.666666667499998, 0.083333332499999, 0.500404040000000],
+        [0.916666667499998, 0.083333332499999, 0.500404040000000],
+        [0.166666667499998, 0.333333332500000, 0.500404040000000],
+        [0.416666667499998, 0.333333332500000, 0.500404040000000],
+        [0.666666667499998, 0.333333332500000, 0.500404040000000],
+        [0.916666667499998, 0.333333332500000, 0.500404040000000],
+        [0.166666667499998, 0.583333332500000, 0.500404040000000],
+        [0.416666667499998, 0.583333332500000, 0.500404040000000],
+        [0.666666667499998, 0.583333332500000, 0.500404040000000],
+        [0.916666667499998, 0.583333332500000, 0.500404040000000],
+        [0.166666667499998, 0.833333332500000, 0.500404040000000],
+        [0.416666667499998, 0.833333332500000, 0.500404040000000],
+        [0.666666667499998, 0.833333332500000, 0.500404040000000],
+        [0.916666667499998, 0.833333332500000, 0.500404040000000],
+    ]
+    numbers = [
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        7,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+        31,
+    ]
+    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)
+    N = 40
+    dfset = np.loadtxt(cwd / "dfset_GaN_442_rd.xz")
+    d = dfset[:, :3].reshape(N, -1, 3)
+    f = dfset[:, 3:].reshape(N, -1, 3)
+    return cell, d, f
+
+
+@pytest.fixture(scope="session")
+def ph3_si_111_fc3() -> tuple[SymfcAtoms, np.ndarray, np.ndarray]:
+    """Return Si-111 data for fc3."""
+    lattice = [
+        [5.433560030000000, 0.000000000000000, 0.000000000000000],
+        [0.000000000000000, 5.433560030000000, 0.000000000000000],
+        [0.000000000000000, 0.000000000000000, 5.433560030000000],
+    ]
+    points = [
+        [0.875000000000000, 0.875000000000000, 0.875000000000000],
+        [0.875000000000000, 0.375000000000000, 0.375000000000000],
+        [0.375000000000000, 0.875000000000000, 0.375000000000000],
+        [0.375000000000000, 0.375000000000000, 0.875000000000000],
+        [0.125000000000000, 0.125000000000000, 0.125000000000000],
+        [0.125000000000000, 0.625000000000000, 0.625000000000000],
+        [0.625000000000000, 0.125000000000000, 0.625000000000000],
+        [0.625000000000000, 0.625000000000000, 0.125000000000000],
+    ]
+    numbers = [
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+        14,
+    ]
+    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)
+    N = 20
+    dfset = np.loadtxt(cwd / "dfset_Si_111_fc3_rd.xz")
+    d = dfset[:, :3].reshape(N, -1, 3)
+    f = dfset[:, 3:].reshape(N, -1, 3)
+    return cell, d, f
diff --git a/tests/dfset_GaN_222_rd.xz b/tests/dfset_GaN_222_rd.xz
new file mode 100644
index 0000000..17179c1
Binary files /dev/null and b/tests/dfset_GaN_222_rd.xz differ
diff --git a/tests/dfset_GaN_442_rd.xz b/tests/dfset_GaN_442_rd.xz
new file mode 100644
index 0000000..5b2aa2c
Binary files /dev/null and b/tests/dfset_GaN_442_rd.xz differ
diff --git a/tests/dfset_NaCl_222_rd.xz b/tests/dfset_NaCl_222_rd.xz
new file mode 100644
index 0000000..3b9c1e4
Binary files /dev/null and b/tests/dfset_NaCl_222_rd.xz differ
diff --git a/tests/dfset_SiO2_221_rd.xz b/tests/dfset_SiO2_221_rd.xz
new file mode 100644
index 0000000..f3a7719
Binary files /dev/null and b/tests/dfset_SiO2_221_rd.xz differ
diff --git a/tests/dfset_Si_111_fc3_rd.xz b/tests/dfset_Si_111_fc3_rd.xz
new file mode 100644
index 0000000..c81d279
Binary files /dev/null and b/tests/dfset_Si_111_fc3_rd.xz differ
diff --git a/tests/dfset_SnO2_223_rd.xz b/tests/dfset_SnO2_223_rd.xz
new file mode 100644
index 0000000..14bb203
Binary files /dev/null and b/tests/dfset_SnO2_223_rd.xz differ
diff --git a/tests/full_fc_NaCl_222.xz b/tests/full_fc_NaCl_222.xz
new file mode 100644
index 0000000..47083f1
Binary files /dev/null and b/tests/full_fc_NaCl_222.xz differ
diff --git a/tests/full_fc_Si_111_fc3_2.xz b/tests/full_fc_Si_111_fc3_2.xz
new file mode 100644
index 0000000..b38ed2b
Binary files /dev/null and b/tests/full_fc_Si_111_fc3_2.xz differ
diff --git a/tests/full_fc_Si_111_fc3_3.xz b/tests/full_fc_Si_111_fc3_3.xz
new file mode 100644
index 0000000..ee0b8d1
Binary files /dev/null and b/tests/full_fc_Si_111_fc3_3.xz differ
diff --git a/tests/full_fc_SnO2_223.xz b/tests/full_fc_SnO2_223.xz
new file mode 100644
index 0000000..fcc8a9d
Binary files /dev/null and b/tests/full_fc_SnO2_223.xz differ
diff --git a/tests/perms_super.dat b/tests/perms_super.dat
new file mode 100644
index 0000000..d45992a
--- /dev/null
+++ b/tests/perms_super.dat
@@ -0,0 +1,96 @@
+0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
+12 14 15 13 8 10 11 9 2 0 1 3 6 4 5 7 24 26 27 25 28 30 31 29 22 20 21 23 18 16 17 19
+2 1 3 0 6 5 7 4 11 8 10 9 15 12 14 13 18 17 19 16 22 21 23 20 27 24 26 25 31 28 30 29
+15 14 13 12 11 10 9 8 3 2 1 0 7 6 5 4 27 26 25 24 31 30 29 28 23 22 21 20 19 18 17 16
+3 1 0 2 7 5 4 6 9 11 10 8 13 15 14 12 19 17 16 18 23 21 20 22 25 27 26 24 29 31 30 28
+13 14 12 15 9 10 8 11 0 3 1 2 4 7 5 6 25 26 24 27 29 30 28 31 20 23 21 22 16 19 17 18
+15 14 12 13 11 10 8 9 2 3 1 0 6 7 5 4 27 26 24 25 31 30 28 29 22 23 21 20 18 19 17 16
+2 1 0 3 6 5 4 7 8 11 10 9 12 15 14 13 18 17 16 19 22 21 20 23 24 27 26 25 28 31 30 29
+12 14 13 15 8 10 9 11 0 2 1 3 4 6 5 7 24 26 25 27 28 30 29 31 20 22 21 23 16 18 17 19
+0 1 3 2 4 5 7 6 9 8 10 11 13 12 14 15 16 17 19 18 20 21 23 22 25 24 26 27 29 28 30 31
+13 14 15 12 9 10 11 8 3 0 1 2 7 4 5 6 25 26 27 24 29 30 31 28 23 20 21 22 19 16 17 18
+3 1 2 0 7 5 6 4 11 9 10 8 15 13 14 12 19 17 18 16 23 21 22 20 27 25 26 24 31 29 30 28
+4 5 6 7 0 1 2 3 12 13 14 15 8 9 10 11 20 21 22 23 16 17 18 19 28 29 30 31 24 25 26 27
+8 10 11 9 12 14 15 13 6 4 5 7 2 0 1 3 28 30 31 29 24 26 27 25 18 16 17 19 22 20 21 23
+6 5 7 4 2 1 3 0 15 12 14 13 11 8 10 9 22 21 23 20 18 17 19 16 31 28 30 29 27 24 26 25
+11 10 9 8 15 14 13 12 7 6 5 4 3 2 1 0 31 30 29 28 27 26 25 24 19 18 17 16 23 22 21 20
+7 5 4 6 3 1 0 2 13 15 14 12 9 11 10 8 23 21 20 22 19 17 16 18 29 31 30 28 25 27 26 24
+9 10 8 11 13 14 12 15 4 7 5 6 0 3 1 2 29 30 28 31 25 26 24 27 16 19 17 18 20 23 21 22
+11 10 8 9 15 14 12 13 6 7 5 4 2 3 1 0 31 30 28 29 27 26 24 25 18 19 17 16 22 23 21 20
+6 5 4 7 2 1 0 3 12 15 14 13 8 11 10 9 22 21 20 23 18 17 16 19 28 31 30 29 24 27 26 25
+8 10 9 11 12 14 13 15 4 6 5 7 0 2 1 3 28 30 29 31 24 26 25 27 16 18 17 19 20 22 21 23
+4 5 7 6 0 1 3 2 13 12 14 15 9 8 10 11 20 21 23 22 16 17 19 18 29 28 30 31 25 24 26 27
+9 10 11 8 13 14 15 12 7 4 5 6 3 0 1 2 29 30 31 28 25 26 27 24 19 16 17 18 23 20 21 22
+7 5 6 4 3 1 2 0 15 13 14 12 11 9 10 8 23 21 22 20 19 17 18 16 31 29 30 28 27 25 26 24
+3 2 1 0 7 6 5 4 11 10 9 8 15 14 13 12 19 18 17 16 23 22 21 20 27 26 25 24 31 30 29 28
+15 13 12 14 11 9 8 10 1 3 2 0 5 7 6 4 27 25 24 26 31 29 28 30 21 23 22 20 17 19 18 16
+1 2 0 3 5 6 4 7 8 11 9 10 12 15 13 14 17 18 16 19 21 22 20 23 24 27 25 26 28 31 29 30
+12 13 14 15 8 9 10 11 0 1 2 3 4 5 6 7 24 25 26 27 28 29 30 31 20 21 22 23 16 17 18 19
+0 2 3 1 4 6 7 5 10 8 9 11 14 12 13 15 16 18 19 17 20 22 23 21 26 24 25 27 30 28 29 31
+14 13 15 12 10 9 11 8 3 0 2 1 7 4 6 5 26 25 27 24 30 29 31 28 23 20 22 21 19 16 18 17
+12 13 15 14 8 9 11 10 1 0 2 3 5 4 6 7 24 25 27 26 28 29 31 30 21 20 22 23 17 16 18 19
+1 2 3 0 5 6 7 4 11 8 9 10 15 12 13 14 17 18 19 16 21 22 23 20 27 24 25 26 31 28 29 30
+15 13 14 12 11 9 10 8 3 1 2 0 7 5 6 4 27 25 26 24 31 29 30 28 23 21 22 20 19 17 18 16
+3 2 0 1 7 6 4 5 10 11 9 8 14 15 13 12 19 18 16 17 23 22 20 21 26 27 25 24 30 31 29 28
+14 13 12 15 10 9 8 11 0 3 2 1 4 7 6 5 26 25 24 27 30 29 28 31 20 23 22 21 16 19 18 17
+0 2 1 3 4 6 5 7 8 10 9 11 12 14 13 15 16 18 17 19 20 22 21 23 24 26 25 27 28 30 29 31
+7 6 5 4 3 2 1 0 15 14 13 12 11 10 9 8 23 22 21 20 19 18 17 16 31 30 29 28 27 26 25 24
+11 9 8 10 15 13 12 14 5 7 6 4 1 3 2 0 31 29 28 30 27 25 24 26 17 19 18 16 21 23 22 20
+5 6 4 7 1 2 0 3 12 15 13 14 8 11 9 10 21 22 20 23 17 18 16 19 28 31 29 30 24 27 25 26
+8 9 10 11 12 13 14 15 4 5 6 7 0 1 2 3 28 29 30 31 24 25 26 27 16 17 18 19 20 21 22 23
+4 6 7 5 0 2 3 1 14 12 13 15 10 8 9 11 20 22 23 21 16 18 19 17 30 28 29 31 26 24 25 27
+10 9 11 8 14 13 15 12 7 4 6 5 3 0 2 1 30 29 31 28 26 25 27 24 19 16 18 17 23 20 22 21
+8 9 11 10 12 13 15 14 5 4 6 7 1 0 2 3 28 29 31 30 24 25 27 26 17 16 18 19 21 20 22 23
+5 6 7 4 1 2 3 0 15 12 13 14 11 8 9 10 21 22 23 20 17 18 19 16 31 28 29 30 27 24 25 26
+11 9 10 8 15 13 14 12 7 5 6 4 3 1 2 0 31 29 30 28 27 25 26 24 19 17 18 16 23 21 22 20
+7 6 4 5 3 2 0 1 14 15 13 12 10 11 9 8 23 22 20 21 19 18 16 17 30 31 29 28 26 27 25 24
+10 9 8 11 14 13 12 15 4 7 6 5 0 3 2 1 30 29 28 31 26 25 24 27 16 19 18 17 20 23 22 21
+4 6 5 7 0 2 1 3 12 14 13 15 8 10 9 11 20 22 21 23 16 18 17 19 28 30 29 31 24 26 25 27
+2 3 0 1 6 7 4 5 10 11 8 9 14 15 12 13 18 19 16 17 22 23 20 21 26 27 24 25 30 31 28 29
+14 12 13 15 10 8 9 11 0 2 3 1 4 6 7 5 26 24 25 27 30 28 29 31 20 22 23 21 16 18 19 17
+0 3 1 2 4 7 5 6 9 10 8 11 13 14 12 15 16 19 17 18 20 23 21 22 25 26 24 27 29 30 28 31
+13 12 15 14 9 8 11 10 1 0 3 2 5 4 7 6 25 24 27 26 29 28 31 30 21 20 23 22 17 16 19 18
+1 3 2 0 5 7 6 4 11 9 8 10 15 13 12 14 17 19 18 16 21 23 22 20 27 25 24 26 31 29 28 30
+15 12 14 13 11 8 10 9 2 1 3 0 6 5 7 4 27 24 26 25 31 28 30 29 22 21 23 20 18 17 19 16
+13 12 14 15 9 8 10 11 0 1 3 2 4 5 7 6 25 24 26 27 29 28 30 31 20 21 23 22 16 17 19 18
+0 3 2 1 4 7 6 5 10 9 8 11 14 13 12 15 16 19 18 17 20 23 22 21 26 25 24 27 30 29 28 31
+14 12 15 13 10 8 11 9 2 0 3 1 6 4 7 5 26 24 27 25 30 28 31 29 22 20 23 21 18 16 19 17
+2 3 1 0 6 7 5 4 11 10 8 9 15 14 12 13 18 19 17 16 22 23 21 20 27 26 24 25 31 30 28 29
+15 12 13 14 11 8 9 10 1 2 3 0 5 6 7 4 27 24 25 26 31 28 29 30 21 22 23 20 17 18 19 16
+1 3 0 2 5 7 4 6 9 11 8 10 13 15 12 14 17 19 16 18 21 23 20 22 25 27 24 26 29 31 28 30
+6 7 4 5 2 3 0 1 14 15 12 13 10 11 8 9 22 23 20 21 18 19 16 17 30 31 28 29 26 27 24 25
+10 8 9 11 14 12 13 15 4 6 7 5 0 2 3 1 30 28 29 31 26 24 25 27 16 18 19 17 20 22 23 21
+4 7 5 6 0 3 1 2 13 14 12 15 9 10 8 11 20 23 21 22 16 19 17 18 29 30 28 31 25 26 24 27
+9 8 11 10 13 12 15 14 5 4 7 6 1 0 3 2 29 28 31 30 25 24 27 26 17 16 19 18 21 20 23 22
+5 7 6 4 1 3 2 0 15 13 12 14 11 9 8 10 21 23 22 20 17 19 18 16 31 29 28 30 27 25 24 26
+11 8 10 9 15 12 14 13 6 5 7 4 2 1 3 0 31 28 30 29 27 24 26 25 18 17 19 16 22 21 23 20
+9 8 10 11 13 12 14 15 4 5 7 6 0 1 3 2 29 28 30 31 25 24 26 27 16 17 19 18 20 21 23 22
+4 7 6 5 0 3 2 1 14 13 12 15 10 9 8 11 20 23 22 21 16 19 18 17 30 29 28 31 26 25 24 27
+10 8 11 9 14 12 15 13 6 4 7 5 2 0 3 1 30 28 31 29 26 24 27 25 18 16 19 17 22 20 23 21
+6 7 5 4 2 3 1 0 15 14 12 13 11 10 8 9 22 23 21 20 18 19 17 16 31 30 28 29 27 26 24 25
+11 8 9 10 15 12 13 14 5 6 7 4 1 2 3 0 31 28 29 30 27 24 25 26 17 18 19 16 21 22 23 20
+5 7 4 6 1 3 0 2 13 15 12 14 9 11 8 10 21 23 20 22 17 19 16 18 29 31 28 30 25 27 24 26
+1 0 3 2 5 4 7 6 9 8 11 10 13 12 15 14 17 16 19 18 21 20 23 22 25 24 27 26 29 28 31 30
+13 15 14 12 9 11 10 8 3 1 0 2 7 5 4 6 25 27 26 24 29 31 30 28 23 21 20 22 19 17 16 18
+3 0 2 1 7 4 6 5 10 9 11 8 14 13 15 12 19 16 18 17 23 20 22 21 26 25 27 24 30 29 31 28
+14 15 12 13 10 11 8 9 2 3 0 1 6 7 4 5 26 27 24 25 30 31 28 29 22 23 20 21 18 19 16 17
+2 0 1 3 6 4 5 7 8 10 11 9 12 14 15 13 18 16 17 19 22 20 21 23 24 26 27 25 28 30 31 29
+12 15 13 14 8 11 9 10 1 2 0 3 5 6 4 7 24 27 25 26 28 31 29 30 21 22 20 23 17 18 16 19
+14 15 13 12 10 11 9 8 3 2 0 1 7 6 4 5 26 27 25 24 30 31 29 28 23 22 20 21 19 18 16 17
+3 0 1 2 7 4 5 6 9 10 11 8 13 14 15 12 19 16 17 18 23 20 21 22 25 26 27 24 29 30 31 28
+13 15 12 14 9 11 8 10 1 3 0 2 5 7 4 6 25 27 24 26 29 31 28 30 21 23 20 22 17 19 16 18
+1 0 2 3 5 4 6 7 8 9 11 10 12 13 15 14 17 16 18 19 21 20 22 23 24 25 27 26 28 29 31 30
+12 15 14 13 8 11 10 9 2 1 0 3 6 5 4 7 24 27 26 25 28 31 30 29 22 21 20 23 18 17 16 19
+2 0 3 1 6 4 7 5 10 8 11 9 14 12 15 13 18 16 19 17 22 20 23 21 26 24 27 25 30 28 31 29
+5 4 7 6 1 0 3 2 13 12 15 14 9 8 11 10 21 20 23 22 17 16 19 18 29 28 31 30 25 24 27 26
+9 11 10 8 13 15 14 12 7 5 4 6 3 1 0 2 29 31 30 28 25 27 26 24 19 17 16 18 23 21 20 22
+7 4 6 5 3 0 2 1 14 13 15 12 10 9 11 8 23 20 22 21 19 16 18 17 30 29 31 28 26 25 27 24
+10 11 8 9 14 15 12 13 6 7 4 5 2 3 0 1 30 31 28 29 26 27 24 25 18 19 16 17 22 23 20 21
+6 4 5 7 2 0 1 3 12 14 15 13 8 10 11 9 22 20 21 23 18 16 17 19 28 30 31 29 24 26 27 25
+8 11 9 10 12 15 13 14 5 6 4 7 1 2 0 3 28 31 29 30 24 27 25 26 17 18 16 19 21 22 20 23
+10 11 9 8 14 15 13 12 7 6 4 5 3 2 0 1 30 31 29 28 26 27 25 24 19 18 16 17 23 22 20 21
+7 4 5 6 3 0 1 2 13 14 15 12 9 10 11 8 23 20 21 22 19 16 17 18 29 30 31 28 25 26 27 24
+9 11 8 10 13 15 12 14 5 7 4 6 1 3 0 2 29 31 28 30 25 27 24 26 17 19 16 18 21 23 20 22
+5 4 6 7 1 0 2 3 12 13 15 14 8 9 11 10 21 20 22 23 17 16 18 19 28 29 31 30 24 25 27 26
+8 11 10 9 12 15 14 13 6 5 4 7 2 1 0 3 28 31 30 29 24 27 26 25 18 17 16 19 22 21 20 23
+6 4 7 5 2 0 3 1 14 12 15 13 10 8 11 9 22 20 23 21 18 16 19 17 30 28 31 29 26 24 27 25
diff --git a/tests/phonopy_to_python.py b/tests/phonopy_to_python.py
new file mode 100644
index 0000000..7b572d8
--- /dev/null
+++ b/tests/phonopy_to_python.py
@@ -0,0 +1,39 @@
+"""This is a script to convert the phonopy.yaml to raw data in python."""
+
+import numpy as np
+import phono3py
+import phonopy
+
+ph3 = phono3py.load("phono3py_params_Si-111-222-rd.yaml.xz", produce_fc=False)
+ph = phonopy.load("phonopy_GaN_222_rd.yaml.xz", produce_fc=False)
+
+print("    lattice = [")
+for v in ph.supercell.cell:
+    print("[", ", ".join([f"{x:.15f}" for x in v]), "],")
+print("    ]")
+
+print("    points = [")
+for v in ph.supercell.scaled_positions:
+    print("[", ", ".join([f"{x:.15f}" for x in v]), "],")
+print("    ]")
+
+print("    numbers = [")
+for v in ph.supercell.numbers:
+    print(f"{v},")
+print("    ]")
+
+print("    cell = SymfcAtoms(cell=lattice, scaled_positions=points, numbers=numbers)")
+
+
+print(f"    N = {ph.displacements.shape[0]}")
+print('    dfset = np.loadtxt(cwd / "dfset_Si_111_fc3_rd.xz")')
+print("    d = dfset[:, :3].reshape(N, -1, 3)")
+print("    f = dfset[:, 3:].reshape(N, -1, 3)")
+print("    return cell, d, f")
+
+d = ph.displacements.reshape(-1, 3)
+f = ph.forces.reshape(-1, 3)
+
+with open("dfset_Si_111_fc3_rd", "w") as w:
+    for v in np.c_[(d, f)]:
+        print(" ".join([f"{x:18.15f}" for x in v]), file=w)
diff --git a/tests/spg_reps/test_spg_reps_O1.py b/tests/spg_reps/test_spg_reps_O1.py
new file mode 100644
index 0000000..d9895f3
--- /dev/null
+++ b/tests/spg_reps/test_spg_reps_O1.py
@@ -0,0 +1,86 @@
+"""Tests of SpgRepsO1 class."""
+
+import numpy as np
+
+from symfc.spg_reps import SpgRepsO1
+from symfc.utils.utils import SymfcAtoms
+
+
+def test_spg_reps_o1_get_sigma1_rep_NaCl111(cell_nacl_111: SymfcAtoms):
+    """Test of SpgRepsO1.get_sigma1_rep."""
+    spg_reps_o1 = SpgRepsO1(cell_nacl_111)
+    N = len(cell_nacl_111)
+    for i, _ in enumerate(spg_reps_o1.unique_rotation_indices):
+        s_r = spg_reps_o1.get_sigma1_rep(i).toarray()
+        r_atoms = s_r @ np.arange(N, dtype=int)
+        r_c = spg_reps_o1.r_reps[i]
+        L = spg_reps_o1._lattice.T
+        r = np.linalg.inv(L) @ r_c @ L
+        r_pos = spg_reps_o1._positions @ r.T
+        for j, pos in enumerate(spg_reps_o1._positions):
+            diff = r_pos - pos
+            diff -= np.rint(diff)
+            idx = np.where(np.linalg.norm(diff @ L.T, axis=1) < 1e-10)[0][0]
+            assert r_atoms[j] == idx
+
+
+def test_spg_reps_o1_r_reps_NaCl111(cell_nacl_111: SymfcAtoms):
+    """Test of SpgRepsO1.r_reps."""
+    ref_r_reps = np.array(
+        [
+            [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
+            [-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0],
+            [0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0],
+            [0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0],
+            [-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0],
+            [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0],
+            [0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0],
+            [0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0],
+            [1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0],
+            [-1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
+            [0.0, -1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0],
+            [0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0],
+            [-1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0],
+            [1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0],
+            [0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, -1.0],
+            [0.0, -1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0],
+            [0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, -1.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0],
+            [0.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0],
+            [0.0, 0.0, -1.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0],
+            [0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0, -1.0, 0.0],
+            [0.0, 0.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, 1.0, 0.0, 1.0, 0.0, -1.0, 0.0, 0.0],
+            [0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 1.0, 0.0, 0.0],
+            [0.0, 0.0, -1.0, 1.0, 0.0, 0.0, 0.0, -1.0, 0.0],
+            [0.0, 0.0, 1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, -1.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0],
+            [0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0],
+            [0.0, 0.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, -1.0, 0.0],
+            [0.0, 0.0, -1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0],
+            [0.0, 0.0, 1.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0],
+            [0.0, -1.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0],
+            [-1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 1.0, 0.0],
+            [0.0, -1.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 1.0, 0.0, 0.0],
+            [-1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0],
+            [1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0],
+            [0.0, -1.0, 0.0, 0.0, 0.0, -1.0, 1.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0],
+            [-1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0],
+            [1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0],
+            [0.0, 1.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, 0.0],
+            [0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 1.0, 0.0],
+            [-1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0],
+        ]
+    )
+
+    spg_reps_o1 = SpgRepsO1(cell_nacl_111)
+    for r in spg_reps_o1.r_reps:
+        r_ravel = r.toarray().ravel()
+        ids = np.where(np.linalg.norm(ref_r_reps - r_ravel, axis=1) < 1e-5)[0]
+        assert len(ids) == 1
diff --git a/tests/spg_reps/test_spg_reps_O2.py b/tests/spg_reps/test_spg_reps_O2.py
new file mode 100644
index 0000000..7713d06
--- /dev/null
+++ b/tests/spg_reps/test_spg_reps_O2.py
@@ -0,0 +1,17 @@
+"""Tests of SpgRepsO2 class."""
+
+import numpy as np
+
+from symfc.spg_reps import SpgRepsO2
+from symfc.utils.utils import SymfcAtoms
+
+
+def test_spg_reps_o2(cell_nacl_111: SymfcAtoms):
+    """Test of SpgRepsO2."""
+    spg_reps_o2 = SpgRepsO2(cell_nacl_111)
+    trace_sum = 0
+    for i, _ in enumerate(spg_reps_o2.unique_rotation_indices):
+        s_r = spg_reps_o2.get_sigma2_rep(i).toarray()
+        trace_sum += np.trace(s_r)
+
+    assert trace_sum == 960
diff --git a/tests/spg_reps/test_spg_reps_base.py b/tests/spg_reps/test_spg_reps_base.py
new file mode 100644
index 0000000..351ed75
--- /dev/null
+++ b/tests/spg_reps/test_spg_reps_base.py
@@ -0,0 +1,40 @@
+"""Tests of SpgRepsBase class."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from symfc.spg_reps import SpgRepsBase
+from symfc.utils.utils import SymfcAtoms
+
+
+def test_translation_permutations_NaCl_111(cell_nacl_111: SymfcAtoms):
+    """Test SpgReps.translation_permutations."""
+    cell = cell_nacl_111
+    sym_reps = SpgRepsBase(cell)
+    trans_perms = sym_reps.translation_permutations
+    # for v in trans_perms:
+    #     print("[", ", ".join([f"{x}" for x in v]), "],")
+    perms_ref = [
+        [0, 1, 2, 3, 4, 5, 6, 7],
+        [3, 2, 1, 0, 7, 6, 5, 4],
+        [2, 3, 0, 1, 6, 7, 4, 5],
+        [1, 0, 3, 2, 5, 4, 7, 6],
+    ]
+    done = []
+    for tperm in trans_perms:
+        for i, perm in enumerate(perms_ref):
+            if np.array_equal(tperm, perm):
+                done.append(i)
+                break
+    np.testing.assert_array_equal(np.sort(done), [0, 1, 2, 3])
+
+
+def test_translation_permutations_shape_GaN_222(
+    ph_gan_222: tuple[SymfcAtoms, np.ndarray, np.ndarray]
+):
+    """Test SpgReps.translation_permutations."""
+    supercell, _, _ = ph_gan_222
+    sym_reps = SpgRepsBase(supercell)
+    trans_perms = sym_reps.translation_permutations
+    assert trans_perms.shape == (8, 32)
diff --git a/tests/test_api.py b/tests/test_api.py
new file mode 100644
index 0000000..1cac7cb
--- /dev/null
+++ b/tests/test_api.py
@@ -0,0 +1,86 @@
+"""Tests of Symfc."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from symfc import Symfc
+from symfc.utils.utils import SymfcAtoms
+
+cwd = Path(__file__).parent
+
+
+def test_api_NaCl_222(ph_nacl_222: tuple[SymfcAtoms, np.ndarray, np.ndarray]):
+    """Test Symfc class."""
+    supercell, displacements, forces = ph_nacl_222
+    symfc = Symfc(supercell)
+    symfc._compute_basis_set(2)
+    symfc.displacements = displacements
+    np.testing.assert_array_almost_equal(symfc.displacements, displacements)
+    symfc.forces = forces
+    np.testing.assert_array_almost_equal(symfc.forces, forces)
+    symfc.solve(
+        [
+            2,
+        ]
+    )
+    fc = symfc.force_constants[2]
+    fc_ref = np.loadtxt(cwd / "compact_fc_NaCl_222.xz").reshape(fc.shape)
+    np.testing.assert_allclose(fc, fc_ref)
+
+
+def test_api_NaCl_222_with_dataset(
+    ph_nacl_222: tuple[SymfcAtoms, np.ndarray, np.ndarray]
+):
+    """Test Symfc class with displacements and forces as input."""
+    supercell, displacements, forces = ph_nacl_222
+    symfc = Symfc(
+        supercell,
+        displacements=displacements,
+        forces=forces,
+    ).run(orders=[2])
+    fc = symfc.force_constants[2]
+    fc_ref = np.loadtxt(cwd / "compact_fc_NaCl_222.xz").reshape(fc.shape)
+    np.testing.assert_allclose(fc, fc_ref)
+
+
+def test_api_NaCl_222_exception(ph_nacl_222: tuple[SymfcAtoms, np.ndarray, np.ndarray]):
+    """Test Symfc class with displacements and forces as input."""
+    supercell, _, _ = ph_nacl_222
+    symfc = Symfc(supercell)
+    symfc._compute_basis_set(2)
+    with pytest.raises(RuntimeError):
+        symfc.solve(
+            orders=[
+                2,
+            ]
+        )
+
+
+@pytest.mark.parametrize("is_compact_fc", [True, False])
+def test_api_si_111_fc3(
+    ph3_si_111_fc3: tuple[SymfcAtoms, np.ndarray, np.ndarray], is_compact_fc: bool
+):
+    """Test Symfc class with displacements and forces as input."""
+    supercell, displacements, forces = ph3_si_111_fc3
+    symfc = Symfc(supercell, displacements=displacements, forces=forces).run(
+        orders=[2, 3], is_compact_fc=is_compact_fc
+    )
+    fc2 = symfc.force_constants[2]
+    fc3 = symfc.force_constants[3]
+
+    if is_compact_fc:
+        # np.savetxt(cwd / "compact_fc_Si_111_fc3_2.xz", fc2.ravel())
+        # np.savetxt(cwd / "compact_fc_Si_111_fc3_3.xz", fc3.ravel())
+        fc2_ref = np.loadtxt(cwd / "compact_fc_Si_111_fc3_2.xz").reshape(fc2.shape)
+        fc3_ref = np.loadtxt(cwd / "compact_fc_Si_111_fc3_3.xz").reshape(fc3.shape)
+    else:
+        # np.savetxt(cwd / "full_fc_Si_111_fc3_2.xz", fc2.ravel())
+        # np.savetxt(cwd / "full_fc_Si_111_fc3_3.xz", fc3.ravel())
+        fc2_ref = np.loadtxt(cwd / "full_fc_Si_111_fc3_2.xz").reshape(fc2.shape)
+        fc3_ref = np.loadtxt(cwd / "full_fc_Si_111_fc3_3.xz").reshape(fc3.shape)
+    np.testing.assert_allclose(fc2_ref, fc2, atol=1e-6)
+    np.testing.assert_allclose(fc3_ref, fc3, atol=1e-6)
diff --git a/tests/utils/test_utils.py b/tests/utils/test_utils.py
new file mode 100644
index 0000000..d336668
--- /dev/null
+++ b/tests/utils/test_utils.py
@@ -0,0 +1,73 @@
+"""Tests of matrix manipulating functions."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from symfc.spg_reps import SpgRepsBase
+from symfc.utils.utils import (
+    SymfcAtoms,
+    compute_sg_permutations,
+    get_indep_atoms_by_lat_trans,
+)
+
+cwd = Path(__file__).parent
+
+
+def test_get_indep_atoms_by_lattice_translation(
+    ph_nacl_222: tuple[SymfcAtoms, np.ndarray, np.ndarray]
+):
+    """Test of get_indep_atoms_by_lattice_translation."""
+    supercell, _, _ = ph_nacl_222
+    sym_op_reps = SpgRepsBase(supercell)
+    trans_perms = sym_op_reps.translation_permutations
+    assert trans_perms.shape == (32, 64)
+    indep_atoms = get_indep_atoms_by_lat_trans(trans_perms)
+    np.testing.assert_array_equal(indep_atoms, [0, 32])
+
+
+def test_compute_sg_permutations(
+    ph_gan_222: tuple[SymfcAtoms, np.ndarray, np.ndarray], cell_gan_111: SymfcAtoms
+):
+    """Test compute_sg_permutations."""
+    pytest.importorskip("spglib")
+    from spglib import spglib
+
+    supercell, _, _ = ph_gan_222
+    primitive = cell_gan_111
+    dataset = spglib.get_symmetry_dataset(supercell.totuple())
+    primitive_dataset = spglib.get_symmetry_dataset(primitive.totuple())
+    perms = compute_sg_permutations(
+        primitive.scaled_positions,
+        primitive_dataset["rotations"],
+        primitive_dataset["translations"],
+        primitive.cell,
+    )
+    ref_perms = [
+        [0, 1, 2, 3],
+        [1, 0, 3, 2],
+        [0, 1, 2, 3],
+        [1, 0, 3, 2],
+        [0, 1, 2, 3],
+        [1, 0, 3, 2],
+        [1, 0, 3, 2],
+        [0, 1, 2, 3],
+        [1, 0, 3, 2],
+        [0, 1, 2, 3],
+        [1, 0, 3, 2],
+        [0, 1, 2, 3],
+    ]
+    np.testing.assert_array_equal(ref_perms, perms)
+
+    perms_super = compute_sg_permutations(
+        supercell.scaled_positions,
+        dataset["rotations"],
+        dataset["translations"],
+        supercell.cell,
+    )
+    # np.savetxt("perms_super.dat", perms_super, fmt="%d")
+    perms_super_ref = np.loadtxt(cwd / ".." / "perms_super.dat", dtype=int)
+    np.testing.assert_array_equal(perms_super_ref, perms_super)
diff --git a/tests/utils/test_utils_O2.py b/tests/utils/test_utils_O2.py
new file mode 100644
index 0000000..b5df1ee
--- /dev/null
+++ b/tests/utils/test_utils_O2.py
@@ -0,0 +1,761 @@
+"""Tests of matrix manipulating functions for 2nd order force constants."""
+
+import numpy as np
+
+from symfc.spg_reps import SpgRepsBase
+from symfc.utils.utils import SymfcAtoms
+from symfc.utils.utils_O2 import (
+    _get_atomic_lat_trans_decompr_indices,
+    _get_perm_compr_matrix_reference,
+    get_lat_trans_compr_indices,
+    get_lat_trans_compr_matrix,
+    get_lat_trans_decompr_indices,
+    get_perm_compr_matrix,
+)
+
+
+def test_get_perm_compr_matrix():
+    """Test of get_perm_compr_matrix."""
+    C1 = get_perm_compr_matrix(8)
+    C2 = _get_perm_compr_matrix_reference(8)
+    np.testing.assert_array_almost_equal((C1 @ C1.T).toarray(), (C2 @ C2.T).toarray())
+    np.testing.assert_array_almost_equal((C1.T @ C1).toarray(), (C2.T @ C2).toarray())
+
+
+def test_get_lat_trans_decompr_indices(cell_nacl_111: SymfcAtoms):
+    """Test of get_lat_trans_decompr_indices.
+
+    The one dimensional array with row-size of compr-mat.
+    Every element indicates column position.
+
+    """
+    ref = [
+        0,
+        1,
+        2,
+        3,
+        4,
+        5,
+        6,
+        7,
+        8,
+        9,
+        10,
+        11,
+        12,
+        13,
+        14,
+        15,
+        16,
+        17,
+        18,
+        19,
+        20,
+        21,
+        22,
+        23,
+        24,
+        25,
+        26,
+        27,
+        28,
+        29,
+        30,
+        31,
+        32,
+        33,
+        34,
+        35,
+        36,
+        37,
+        38,
+        39,
+        40,
+        41,
+        42,
+        43,
+        44,
+        45,
+        46,
+        47,
+        48,
+        49,
+        50,
+        51,
+        52,
+        53,
+        54,
+        55,
+        56,
+        57,
+        58,
+        59,
+        60,
+        61,
+        62,
+        63,
+        64,
+        65,
+        66,
+        67,
+        68,
+        69,
+        70,
+        71,
+        9,
+        10,
+        11,
+        12,
+        13,
+        14,
+        15,
+        16,
+        17,
+        0,
+        1,
+        2,
+        3,
+        4,
+        5,
+        6,
+        7,
+        8,
+        27,
+        28,
+        29,
+        30,
+        31,
+        32,
+        33,
+        34,
+        35,
+        18,
+        19,
+        20,
+        21,
+        22,
+        23,
+        24,
+        25,
+        26,
+        45,
+        46,
+        47,
+        48,
+        49,
+        50,
+        51,
+        52,
+        53,
+        36,
+        37,
+        38,
+        39,
+        40,
+        41,
+        42,
+        43,
+        44,
+        63,
+        64,
+        65,
+        66,
+        67,
+        68,
+        69,
+        70,
+        71,
+        54,
+        55,
+        56,
+        57,
+        58,
+        59,
+        60,
+        61,
+        62,
+        18,
+        19,
+        20,
+        21,
+        22,
+        23,
+        24,
+        25,
+        26,
+        27,
+        28,
+        29,
+        30,
+        31,
+        32,
+        33,
+        34,
+        35,
+        0,
+        1,
+        2,
+        3,
+        4,
+        5,
+        6,
+        7,
+        8,
+        9,
+        10,
+        11,
+        12,
+        13,
+        14,
+        15,
+        16,
+        17,
+        54,
+        55,
+        56,
+        57,
+        58,
+        59,
+        60,
+        61,
+        62,
+        63,
+        64,
+        65,
+        66,
+        67,
+        68,
+        69,
+        70,
+        71,
+        36,
+        37,
+        38,
+        39,
+        40,
+        41,
+        42,
+        43,
+        44,
+        45,
+        46,
+        47,
+        48,
+        49,
+        50,
+        51,
+        52,
+        53,
+        27,
+        28,
+        29,
+        30,
+        31,
+        32,
+        33,
+        34,
+        35,
+        18,
+        19,
+        20,
+        21,
+        22,
+        23,
+        24,
+        25,
+        26,
+        9,
+        10,
+        11,
+        12,
+        13,
+        14,
+        15,
+        16,
+        17,
+        0,
+        1,
+        2,
+        3,
+        4,
+        5,
+        6,
+        7,
+        8,
+        63,
+        64,
+        65,
+        66,
+        67,
+        68,
+        69,
+        70,
+        71,
+        54,
+        55,
+        56,
+        57,
+        58,
+        59,
+        60,
+        61,
+        62,
+        45,
+        46,
+        47,
+        48,
+        49,
+        50,
+        51,
+        52,
+        53,
+        36,
+        37,
+        38,
+        39,
+        40,
+        41,
+        42,
+        43,
+        44,
+        72,
+        73,
+        74,
+        75,
+        76,
+        77,
+        78,
+        79,
+        80,
+        81,
+        82,
+        83,
+        84,
+        85,
+        86,
+        87,
+        88,
+        89,
+        90,
+        91,
+        92,
+        93,
+        94,
+        95,
+        96,
+        97,
+        98,
+        99,
+        100,
+        101,
+        102,
+        103,
+        104,
+        105,
+        106,
+        107,
+        108,
+        109,
+        110,
+        111,
+        112,
+        113,
+        114,
+        115,
+        116,
+        117,
+        118,
+        119,
+        120,
+        121,
+        122,
+        123,
+        124,
+        125,
+        126,
+        127,
+        128,
+        129,
+        130,
+        131,
+        132,
+        133,
+        134,
+        135,
+        136,
+        137,
+        138,
+        139,
+        140,
+        141,
+        142,
+        143,
+        81,
+        82,
+        83,
+        84,
+        85,
+        86,
+        87,
+        88,
+        89,
+        72,
+        73,
+        74,
+        75,
+        76,
+        77,
+        78,
+        79,
+        80,
+        99,
+        100,
+        101,
+        102,
+        103,
+        104,
+        105,
+        106,
+        107,
+        90,
+        91,
+        92,
+        93,
+        94,
+        95,
+        96,
+        97,
+        98,
+        117,
+        118,
+        119,
+        120,
+        121,
+        122,
+        123,
+        124,
+        125,
+        108,
+        109,
+        110,
+        111,
+        112,
+        113,
+        114,
+        115,
+        116,
+        135,
+        136,
+        137,
+        138,
+        139,
+        140,
+        141,
+        142,
+        143,
+        126,
+        127,
+        128,
+        129,
+        130,
+        131,
+        132,
+        133,
+        134,
+        90,
+        91,
+        92,
+        93,
+        94,
+        95,
+        96,
+        97,
+        98,
+        99,
+        100,
+        101,
+        102,
+        103,
+        104,
+        105,
+        106,
+        107,
+        72,
+        73,
+        74,
+        75,
+        76,
+        77,
+        78,
+        79,
+        80,
+        81,
+        82,
+        83,
+        84,
+        85,
+        86,
+        87,
+        88,
+        89,
+        126,
+        127,
+        128,
+        129,
+        130,
+        131,
+        132,
+        133,
+        134,
+        135,
+        136,
+        137,
+        138,
+        139,
+        140,
+        141,
+        142,
+        143,
+        108,
+        109,
+        110,
+        111,
+        112,
+        113,
+        114,
+        115,
+        116,
+        117,
+        118,
+        119,
+        120,
+        121,
+        122,
+        123,
+        124,
+        125,
+        99,
+        100,
+        101,
+        102,
+        103,
+        104,
+        105,
+        106,
+        107,
+        90,
+        91,
+        92,
+        93,
+        94,
+        95,
+        96,
+        97,
+        98,
+        81,
+        82,
+        83,
+        84,
+        85,
+        86,
+        87,
+        88,
+        89,
+        72,
+        73,
+        74,
+        75,
+        76,
+        77,
+        78,
+        79,
+        80,
+        135,
+        136,
+        137,
+        138,
+        139,
+        140,
+        141,
+        142,
+        143,
+        126,
+        127,
+        128,
+        129,
+        130,
+        131,
+        132,
+        133,
+        134,
+        117,
+        118,
+        119,
+        120,
+        121,
+        122,
+        123,
+        124,
+        125,
+        108,
+        109,
+        110,
+        111,
+        112,
+        113,
+        114,
+        115,
+        116,
+    ]
+    trans_perms_ref = np.array(
+        [
+            [0, 1, 2, 3, 4, 5, 6, 7],
+            [3, 2, 1, 0, 7, 6, 5, 4],
+            [2, 3, 0, 1, 6, 7, 4, 5],
+            [1, 0, 3, 2, 5, 4, 7, 6],
+        ],
+        dtype=int,
+    )
+
+    spg_reps = SpgRepsBase(cell_nacl_111)
+    trans_perms = spg_reps.translation_permutations
+    assert trans_perms.shape == (4, 8)
+    check_trans_perms(trans_perms, trans_perms_ref)
+    decompr_idx = get_lat_trans_decompr_indices(trans_perms_ref)
+    np.testing.assert_array_equal(ref, decompr_idx)
+
+
+def test_get_lat_trans_compr_indices(cell_nacl_111: SymfcAtoms):
+    """Test get_lat_trans_compr_indices.
+
+    The two dimensional array (n_a * N * 9, n_lp) stores NN33 indices where
+    compression matrix elements are non-zero.
+
+    """
+    trans_perms_ref = np.array(
+        [
+            [0, 1, 2, 3, 4, 5, 6, 7],
+            [3, 2, 1, 0, 7, 6, 5, 4],
+            [2, 3, 0, 1, 6, 7, 4, 5],
+            [1, 0, 3, 2, 5, 4, 7, 6],
+        ],
+        dtype=int,
+    )
+
+    spg_reps = SpgRepsBase(cell_nacl_111)
+    trans_perms = spg_reps.translation_permutations
+    n_lp, N = trans_perms.shape
+    assert trans_perms.shape == (4, 8)
+    check_trans_perms(trans_perms, trans_perms_ref)
+    decompr_mat = get_lat_trans_decompr_indices(trans_perms_ref)
+    compr_mat = get_lat_trans_compr_matrix(decompr_mat, N, n_lp).toarray()
+    compr_idx = get_lat_trans_compr_indices(trans_perms_ref)
+    for c, elem_idx in enumerate(compr_idx):
+        for r in elem_idx:
+            np.testing.assert_almost_equal(compr_mat[r, c], 0.5)
+
+
+def test_get_atomic_lat_trans_decompr_indices(cell_nacl_111: SymfcAtoms):
+    """Test of get_atomic_lat_trans_decompr_indices.
+
+    This function is an atomic version of get_lat_trans_decompr_indices.
+    The one dimensional array with row-size of compr-mat.
+    Every element indicates column position.
+
+    """
+    ref = [
+        0,
+        1,
+        2,
+        3,
+        4,
+        5,
+        6,
+        7,
+        1,
+        0,
+        3,
+        2,
+        5,
+        4,
+        7,
+        6,
+        2,
+        3,
+        0,
+        1,
+        6,
+        7,
+        4,
+        5,
+        3,
+        2,
+        1,
+        0,
+        7,
+        6,
+        5,
+        4,
+        8,
+        9,
+        10,
+        11,
+        12,
+        13,
+        14,
+        15,
+        9,
+        8,
+        11,
+        10,
+        13,
+        12,
+        15,
+        14,
+        10,
+        11,
+        8,
+        9,
+        14,
+        15,
+        12,
+        13,
+        11,
+        10,
+        9,
+        8,
+        15,
+        14,
+        13,
+        12,
+    ]
+    trans_perms_ref = np.array(
+        [
+            [0, 1, 2, 3, 4, 5, 6, 7],
+            [3, 2, 1, 0, 7, 6, 5, 4],
+            [2, 3, 0, 1, 6, 7, 4, 5],
+            [1, 0, 3, 2, 5, 4, 7, 6],
+        ],
+        dtype=int,
+    )
+
+    spg_reps = SpgRepsBase(cell_nacl_111)
+    trans_perms = spg_reps.translation_permutations
+    assert trans_perms.shape == (4, 8)
+    check_trans_perms(trans_perms, trans_perms_ref)
+    atomic_decompr_idx = _get_atomic_lat_trans_decompr_indices(trans_perms_ref)
+    np.testing.assert_array_equal(ref, atomic_decompr_idx)
+
+
+def check_trans_perms(trans_perms, trans_perms_ref):
+    """Check the order of rows of trans_perms.
+
+    The row order of trans_perms may depend on symmetry finding condition.
+
+    """
+    idxs = []
+    for row in trans_perms:
+        match = np.where(((trans_perms_ref - row) == 0).all(axis=1))[0]
+        assert len(match) == 1
+        idxs.append(match[0])
+    np.testing.assert_array_equal(np.sort(idxs), range(len(idxs)))