Merge branch 'rc'

.pre-commit-config.yaml

@@ -11,7 +11,7 @@ repos:
   - id: check-added-large-files
 
 - repo: https://github.com/pycqa/flake8
-  rev: 6.1.0
+  rev: 7.0.0
   hooks:
   - id: flake8
     exclude: |
@@ -35,7 +35,7 @@ repos:
     - "--ignore=E203,W503,E501"
 
 - repo: https://github.com/psf/black
-  rev: 23.10.1
+  rev: 24.2.0
   hooks:
   - id: black
     args:
@@ -47,7 +47,7 @@ repos:
   - id: pydocstyle
 
 - repo: https://github.com/pycqa/isort
-  rev: 5.12.0
+  rev: 5.13.2
   hooks:
     - id: isort
       name: isort (python)

doc/changelog.md

@@ -2,6 +2,12 @@
 
 # Change Log
 
+## Feb-14-2024: Version 2.21.1
+
+- `np.random.randn` is replaced by `np.random.Generator.standard_normal`. Note
+  that by this chang, random seed became incompatible.
+- Some bug fixes.
+
 ## Dec-4-2023: Version 2.21.0
 
 - Maintenance release.

doc/conf.py

@@ -1,4 +1,5 @@
 """Sphinx configuration of phonopy documentation."""
+
 # -*- coding: utf-8 -*-
 #
 # phonopy documentation build configuration file, created by
@@ -52,9 +53,9 @@ copyright = "2009, Atsushi Togo"
 # built documents.
 #
 # The short X.Y version.
-version = "2.20"
+version = "2.21"
 # The full version, including alpha/beta/rc tags.
-release = "2.20.0"
+release = "2.21.1"
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

doc/lammps.md

@@ -14,12 +14,14 @@
 ## Supported LAMMPS structure input format
 
 Two important limitations are:
+
 - Crystal structure has to be described in the similar format for
   [read_data](https://docs.lammps.org/read_data.html).
 - Basis vectors are rotated to match the LAMMPS structure file format of
   [triclinic simulation box](https://docs.lammps.org/Howto_triclinic.html).
 
 ### Supported `read_data` keywords in the header
+
 ```
 atoms
 atom types
@@ -102,7 +104,7 @@ run 0
 
 where `supercell-001` is that generated by phonopy. The `pair_style` of
 `polymlp` is a LAMMPS module of the polynomial machine learning potentials
-provided at https://sekocha.github.io/lammps/index-e.html. For the HCP Ti
+provided at <https://sekocha.github.io/lammps/index-e.html>. For the HCP Ti
 calculation found in the [example
 directory](https://github.com/phonopy/phonopy/tree/develop/example), `mlp.lammp`
 of gtinv-294 was obtained from [Polynomial Machine Learning Potential Repository
@@ -113,7 +115,7 @@ University](http://cms.mtl.kyoto-u.ac.jp/seko/mlp-repository/index.html).
 
 1. Read a lammps input structure file and create supercells with
 
-   ```bash
+   ```
    % phonopy --lammps -c lammps_structure_Ti -d --dim 4 4 3
            _
      _ __ | |__   ___  _ __   ___   _ __  _   _
@@ -157,7 +159,7 @@ University](http://cms.mtl.kyoto-u.ac.jp/seko/mlp-repository/index.html).
 
 3. Make `FORCE_SETS`
 
-   ```bash
+   ```
    % phonopy -f lammps_forces_Ti.0
            _
      _ __ | |__   ___  _ __   ___   _ __  _   _
@@ -252,4 +254,4 @@ write_data dump.unitcell
 ```
 
 More instruction is found at
-https://gist.github.com/lan496/e9dff8449cd7489f6722b276282e66a0.
+<https://gist.github.com/lan496/e9dff8449cd7489f6722b276282e66a0>.

example/Al-QHA/Al-QHA.py

@@ -1,4 +1,5 @@
 """Example of QHA calculation by Al."""
+
 import numpy as np
 import yaml
 from yaml import CLoader as Loader

example/Al2O3/Al2O3.py

@@ -1,4 +1,5 @@
 """Example by corundum Al2O3."""
+
 import numpy as np
 
 import phonopy

example/MgB2/MgB2.py

@@ -1,4 +1,5 @@
 """Example of calculation of irreps of MgB2."""
+
 import numpy as np
 
 import phonopy

example/NaCl/NaCl-band-gv.py

@@ -1,4 +1,5 @@
 """Example of band structure and group velocity calculation by NaCl."""
+
 from typing import List
 
 import numpy as np

example/NaCl/NaCl-band.py

@@ -1,4 +1,5 @@
 """NaCl band structure calculation example."""
+
 import phonopy
 from phonopy.phonon.band_structure import get_band_qpoints
 

example/NaCl/NaCl-gv.py

@@ -1,4 +1,5 @@
 """Group velocity example by NaCl."""
+
 import matplotlib.pyplot as plt
 import numpy as np
 

example/NaCl/NaCl-read_write_fc.py

@@ -1,4 +1,5 @@
 """Example to read and write FORCE_CONSTANTS file."""
+
 import phonopy
 from phonopy.file_IO import parse_FORCE_CONSTANTS, write_FORCE_CONSTANTS
 

example/NaCl/NaCl-yaml.py

@@ -1,4 +1,5 @@
 """Example to obtain PhonopyYaml instance."""
+
 import phonopy
 from phonopy.interface.phonopy_yaml import PhonopyYaml
 

example/NaCl/NaCl.py

@@ -1,4 +1,5 @@
 """Example of NaCl calculation."""
+
 from typing import List
 
 import numpy as np

example/Si-QHA/Si.py

@@ -1,4 +1,5 @@
 """QHA example of Si."""
+
 import tarfile
 
 import matplotlib.pyplot as plt

example/Si-gruneisen/Si-gruneisen.py

@@ -1,4 +1,5 @@
 """Example to calculate mode Grueneisen parameters."""
+
 from typing import List
 
 import numpy as np

example/Si-lammps/generate_displacements.py

@@ -1,4 +1,5 @@
 """A script to generate supercells with displacements for LAMMPS."""
+
 import phonopy
 from phonopy.interface.calculator import write_supercells_with_displacements
 from phonopy.interface.phonopy_yaml import read_cell_yaml

phonopy/__init__.py

@@ -1,4 +1,5 @@
 """Phonon calculation code: Phonopy."""
+
 # Copyright (C) 2015 Atsushi Togo
 # All rights reserved.
 #

phonopy/api_gruneisen.py

@@ -1,4 +1,5 @@
 """API of mode Grueneisen parameter calculation."""
+
 # Copyright (C) 2015 Atsushi Togo
 # All rights reserved.
 #

phonopy/api_phonopy.py

@@ -1,4 +1,5 @@
 """Phonopy class."""
+
 # Copyright (C) 2015 Atsushi Togo
 # All rights reserved.
 #
@@ -38,6 +39,7 @@ from __future__ import annotations
 import sys
 import textwrap
 import warnings
+from collections.abc import Sequence
 from typing import Optional, Union
 
 import numpy as np
@@ -48,7 +50,11 @@ from phonopy.harmonic.displacement import (
     get_least_displacements,
     get_random_displacements_dataset,
 )
-from phonopy.harmonic.dynamical_matrix import DynamicalMatrixGL, get_dynamical_matrix
+from phonopy.harmonic.dynamical_matrix import (
+    DynamicalMatrix,
+    DynamicalMatrixGL,
+    get_dynamical_matrix,
+)
 from phonopy.harmonic.dynmat_to_fc import DynmatToForceConstants
 from phonopy.harmonic.force_constants import cutoff_force_constants
 from phonopy.harmonic.force_constants import get_fc2 as get_phonopy_fc2
@@ -80,6 +86,7 @@ from phonopy.spectrum.dynamic_structure_factor import DynamicStructureFactor
 from phonopy.structure.atoms import PhonopyAtoms
 from phonopy.structure.cells import (
     Primitive,
+    Supercell,
     get_primitive,
     get_primitive_matrix,
     get_supercell,
@@ -139,8 +146,8 @@ class Phonopy:
     def __init__(
         self,
         unitcell,
-        supercell_matrix=None,
-        primitive_matrix=None,
+        supercell_matrix: Optional[Union[Sequence, np.ndarray]] = None,
+        primitive_matrix: Optional[Union[str, Sequence, np.ndarray]] = None,
         nac_params=None,
         factor=VaspToTHz,
         frequency_scale_factor=None,
@@ -245,7 +252,7 @@ class Phonopy:
         self._gv_delta_q = group_velocity_delta_q
 
     @property
-    def version(self):
+    def version(self) -> str:
         """Return phonopy release version number.
 
         str
@@ -281,7 +288,7 @@ class Phonopy:
         return self.primitive
 
     @property
-    def unitcell(self):
+    def unitcell(self) -> PhonopyAtoms:
         """Return input unit cell.
 
         PhonopyAtoms
@@ -313,7 +320,7 @@ class Phonopy:
         self._displacement_dataset = None
 
     @property
-    def supercell(self):
+    def supercell(self) -> Supercell:
         """Return supercell.
 
         Supercell
@@ -331,7 +338,7 @@ class Phonopy:
         return self.supercell
 
     @property
-    def symmetry(self):
+    def symmetry(self) -> Symmetry:
         """Return symmetry of supercell.
 
         Symmetry
@@ -349,7 +356,7 @@ class Phonopy:
         return self.symmetry
 
     @property
-    def primitive_symmetry(self):
+    def primitive_symmetry(self) -> Symmetry:
         """Return symmetry of primitive cell.
 
         Symmetry
@@ -368,7 +375,7 @@ class Phonopy:
         return self.primitive_symmetry
 
     @property
-    def supercell_matrix(self):
+    def supercell_matrix(self) -> np.ndarray:
         """Return transformation matrix to supercell cell from unit cell.
 
         ndarray
@@ -388,7 +395,7 @@ class Phonopy:
         return self.supercell_matrix
 
     @property
-    def primitive_matrix(self):
+    def primitive_matrix(self) -> np.ndarray:
         """Return transformation matrix to primitive cell from unit cell.
 
         ndarray
@@ -408,7 +415,7 @@ class Phonopy:
         return self.primitive_matrix
 
     @property
-    def unit_conversion_factor(self):
+    def unit_conversion_factor(self) -> float:
         """Return phonon frequency unit conversion factor.
 
         float
@@ -432,7 +439,7 @@ class Phonopy:
         return self.unit_conversion_factor
 
     @property
-    def calculator(self):
+    def calculator(self) -> str:
         """Return calculator name.
 
         str
@@ -442,7 +449,7 @@ class Phonopy:
         return self._calculator
 
     @property
-    def dataset(self):
+    def dataset(self) -> dict:
         """Return dataset to store displacements and forces.
 
         Dataset containing information of displacements in supercells.
@@ -521,7 +528,7 @@ class Phonopy:
         self.dataset = displacement_dataset
 
     @property
-    def displacements(self):
+    def displacements(self) -> Union[np.ndarray, list]:
         """Getter and setter of displacements in supercells.
 
         There are two types of displacement dataset. See the docstring
@@ -581,7 +588,7 @@ class Phonopy:
         return self.displacements
 
     @property
-    def force_constants(self):
+    def force_constants(self) -> np.ndarray:
         """Getter and setter of supercell force constants.
 
         Force constants matrix.
@@ -654,7 +661,7 @@ class Phonopy:
             self._set_dynamical_matrix()
 
     @property
-    def forces(self):
+    def forces(self) -> np.ndarray:
         """Return forces of supercells.
 
         ndarray to get and array_like to set
@@ -708,7 +715,7 @@ class Phonopy:
         self.forces = sets_of_forces
 
     @property
-    def dynamical_matrix(self):
+    def dynamical_matrix(self) -> DynamicalMatrix:
         """Return DynamicalMatrix instance.
 
         This is not dynamical matrices but the instance of DynamicalMatrix
@@ -727,7 +734,7 @@ class Phonopy:
         return self.dynamical_matrix
 
     @property
-    def nac_params(self):
+    def nac_params(self) -> dict:
         """Getter and setter of parameters for non-analytical term correction.
 
         dict
@@ -769,7 +776,7 @@ class Phonopy:
         self.nac_params = nac_params
 
     @property
-    def supercells_with_displacements(self):
+    def supercells_with_displacements(self) -> list[PhonopyAtoms]:
         """Return supercells with displacements.
 
         list of PhonopyAtoms
@@ -794,7 +801,7 @@ class Phonopy:
         return self.supercells_with_displacements
 
     @property
-    def mesh_numbers(self):
+    def mesh_numbers(self) -> np.ndarray:
         """Return sampling mesh numbers in reciprocal space."""
         if self._mesh is None:
             return None
@@ -802,62 +809,62 @@ class Phonopy:
             return self._mesh.mesh_numbers
 
     @property
-    def qpoints(self):
+    def qpoints(self) -> QpointsPhonon:
         """Return QpointsPhonon instance."""
         return self._qpoints
 
     @property
-    def band_structure(self):
+    def band_structure(self) -> BandStructure:
         """Return BandStructure instance."""
         return self._band_structure
 
     @property
-    def group_velocity(self):
+    def group_velocity(self) -> GroupVelocity:
         """Return GroupVelocity instance."""
         return self._group_velocity
 
     @property
-    def mesh(self):
+    def mesh(self) -> Union[Mesh, IterMesh]:
         """Return Mesh or IterMesh instance."""
         return self._mesh
 
     @property
-    def random_displacements(self):
+    def random_displacements(self) -> RandomDisplacements:
         """Return RandomDisplacements instance."""
         return self._random_displacements
 
     @property
-    def dynamic_structure_factor(self):
+    def dynamic_structure_factor(self) -> DynamicStructureFactor:
         """Return DynamicStructureFactor instance."""
         return self._dynamic_structure_factor
 
     @property
-    def thermal_properties(self):
+    def thermal_properties(self) -> ThermalProperties:
         """Return ThermalProperties instance."""
         return self._thermal_properties
 
     @property
-    def thermal_displacements(self):
+    def thermal_displacements(self) -> ThermalDisplacements:
         """Return ThermalDisplacements instance."""
         return self._thermal_displacements
 
     @property
-    def thermal_displacement_matrices(self):
+    def thermal_displacement_matrices(self) -> ThermalDisplacementMatrices:
         """Return ThermalDisplacementMatrices instance."""
         return self._thermal_displacement_matrices
 
     @property
-    def irreps(self):
+    def irreps(self) -> IrReps:
         """Return IrReps instance."""
         return self._irreps
 
     @property
-    def moment(self):
+    def moment(self) -> PhononMoment:
         """Return PhononMoment instance."""
         return self._moment
 
     @property
-    def total_dos(self):
+    def total_dos(self) -> TotalDos:
         """Return TotalDos instance."""
         return self._total_dos
 
@@ -871,12 +878,12 @@ class Phonopy:
         return self.projected_dos
 
     @property
-    def projected_dos(self):
-        """Return PartialDos instance."""
+    def projected_dos(self) -> ProjectedDos:
+        """Return ProjectedDOS instance."""
         return self._pdos
 
     @property
-    def masses(self):
+    def masses(self) -> np.ndarray:
         """Getter and setter of masses of primitive cell atoms."""
         return self._primitive.masses
 
@@ -899,7 +906,7 @@ class Phonopy:
 
     def generate_displacements(
         self,
-        distance: float = 0.01,
+        distance: Optional[float] = None,
         is_plusminus: Union[str, bool] = "auto",
         is_diagonal: bool = True,
         is_trigonal: bool = False,
@@ -908,7 +915,9 @@ class Phonopy:
         temperature: Optional[float] = None,
         cutoff_frequency: Optional[float] = None,
         max_distance: Optional[float] = None,
-    ):
+        is_random_distance: bool = False,
+        min_distance: Optional[float] = None,
+    ) -> None:
         """Generate displacement dataset.
 
         There are two modes, finite difference method with systematic
@@ -926,7 +935,9 @@ class Phonopy:
         ----------
         distance : float, optional
             Displacement distance. Unit is the same as that used for crystal
-            structure. Default is 0.01.
+            structure. Default is 0.01. For random direction and distance
+            displacements generation, this value is used when `max_distance` is
+            unspecified.
         is_plusminus : 'auto', True, or False, optional
             For each atom, displacement of one direction (False), both
             direction, i.e., one directiona and its opposite direction (True),
@@ -944,7 +955,7 @@ class Phonopy:
             'distance' parameter, i.e., all atoms in supercell are displaced
             with the same displacement distance in direct space. Default is
             None.
-        random_seed : 32bit unsigned int or None, optional
+        random_seed : int or None, optional
             Random seed for random displacements generation. Default is None.
         temperature : float or None, optional
             With given temperature, random displacements at temperature is
@@ -961,7 +972,17 @@ class Phonopy:
             In random displacements generation from canonical ensemble of
             harmonic phonons, displacements larger than max distance are
             renormalized to the max distance, i.e., a disptalcement d is shorten
-            by d -> d / |d| * max_distance if |d| > max_distance.
+            by d -> d / |d| * max_distance if |d| > max_distance. In random
+            direction and distance displacements generation, this value is is
+            specified.
+        is_random_distance : bool, optional
+            Random direction displacements are generated also with random
+            amplitudes. The maximum value is defined by `distance` and minimum
+            value is given by `min_distance`. Default is False.
+        min_distance : float or None, optional
+            In random direction displacements generation with random distance
+            (`is_random_distance=True`), the minimum distance is given by this
+            value.
 
         """
         if number_of_snapshots is not None and number_of_snapshots > 0:
@@ -972,12 +993,21 @@ class Phonopy:
                 _random_seed = None
                 displacement_dataset = {}
             if temperature is None:
+                if max_distance is None:
+                    if distance is None:
+                        _distance = 0.01
+                    else:
+                        _distance = distance
+                else:
+                    _distance = max_distance
                 d = get_random_displacements_dataset(
                     number_of_snapshots,
-                    distance,
                     len(self._supercell),
+                    _distance,
                     random_seed=_random_seed,
                     is_plusminus=(is_plusminus is True),
+                    is_random_distance=is_random_distance,
+                    min_distance=min_distance,
                 )
                 displacement_dataset["displacements"] = d
             else:
@@ -992,6 +1022,10 @@ class Phonopy:
                 )
                 displacement_dataset["displacements"] = d
         else:
+            if distance is None:
+                _distance = 0.01
+            else:
+                _distance = distance
             displacement_directions = get_least_displacements(
                 self._symmetry,
                 is_plusminus=is_plusminus,
@@ -1000,7 +1034,7 @@ class Phonopy:
                 log_level=self._log_level,
             )
             displacement_dataset = directions_to_displacement_dataset(
-                displacement_directions, distance, self._supercell
+                displacement_directions, _distance, self._supercell
             )
         self.dataset = displacement_dataset
 
@@ -1011,7 +1045,7 @@ class Phonopy:
         fc_calculator=None,
         fc_calculator_options=None,
         show_drift=True,
-    ):
+    ) -> None:
         """Compute supercell force constants from forces-displacements dataset.
 
         Supercell force constants are computed from forces and displacements.
@@ -1068,7 +1102,7 @@ class Phonopy:
         if self._primitive.masses is not None:
             self._set_dynamical_matrix()
 
-    def symmetrize_force_constants(self, level=1, show_drift=True):
+    def symmetrize_force_constants(self, level=1, show_drift=True) -> None:
         """Symmetrize force constants.
 
         This applies translational and permutation symmetries successfully,
@@ -1101,7 +1135,7 @@ class Phonopy:
         if self._primitive.masses is not None:
             self._set_dynamical_matrix()
 
-    def symmetrize_force_constants_by_space_group(self, show_drift=True):
+    def symmetrize_force_constants_by_space_group(self, show_drift=True) -> None:
         """Symmetrize force constants using space group operations.
 
         Space group operations except for pure translations are applied
@@ -1134,7 +1168,7 @@ class Phonopy:
     #####################
 
     # Single q-point
-    def get_dynamical_matrix_at_q(self, q):
+    def get_dynamical_matrix_at_q(self, q) -> np.ndarray:
         """Calculate dynamical matrix at a given q-point.
 
         Parameters
@@ -1158,9 +1192,9 @@ class Phonopy:
             raise RuntimeError(msg)
 
         self._dynamical_matrix.run(q)
-        return self._dynamical_matrix.get_dynamical_matrix()
+        return self._dynamical_matrix.dynamical_matrix
 
-    def get_frequencies(self, q):
+    def get_frequencies(self, q) -> np.ndarray:
         """Calculate phonon frequencies at a given q-point.
 
         Parameters
@@ -1193,7 +1227,7 @@ class Phonopy:
 
         return np.array(frequencies) * self._factor
 
-    def get_frequencies_with_eigenvectors(self, q):
+    def get_frequencies_with_eigenvectors(self, q) -> np.ndarray:
         """Calculate phonon frequencies and eigenvectors at a given q-point.
 
         Parameters
@@ -1245,7 +1279,7 @@ class Phonopy:
         path_connections=None,
         labels=None,
         is_legacy_plot=False,
-    ):
+    ) -> None:
         """Run phonon band structure calculation.
 
         Parameters
@@ -1329,7 +1363,7 @@ class Phonopy:
             is_legacy_plot=is_legacy_plot,
         )
 
-    def get_band_structure_dict(self):
+    def get_band_structure_dict(self) -> dict:
         """Return calculated band structures.
 
         Returns
@@ -1509,7 +1543,7 @@ class Phonopy:
         self._band_structure.plot(axs)
         return plt
 
-    def write_hdf5_band_structure(self, comment=None, filename="band.hdf5"):
+    def write_hdf5_band_structure(self, comment=None, filename="band.hdf5") -> None:
         """Write band structure in hdf5 format.
 
         Parameters
@@ -1522,7 +1556,9 @@ class Phonopy:
         """
         self._band_structure.write_hdf5(comment=comment, filename=filename)
 
-    def write_yaml_band_structure(self, comment=None, filename=None, compression=None):
+    def write_yaml_band_structure(
+        self, comment=None, filename=None, compression=None
+    ) -> None:
         """Write band structure in yaml.
 
         Parameters
@@ -1553,7 +1589,7 @@ class Phonopy:
         with_group_velocities=False,
         is_gamma_center=False,
         use_iter_mesh=False,
-    ):
+    ) -> None:
         """Initialize mesh sampling phonon calculation without starting to run.
 
         Phonon calculation starts explicitly with calling Mesh.run() or
@@ -1662,7 +1698,7 @@ class Phonopy:
         with_eigenvectors=False,
         with_group_velocities=False,
         is_gamma_center=False,
-    ):
+    ) -> None:
         """Run mesh sampling phonon calculation.
 
         See the parameter details in Phonopy.init_mesh.
@@ -1749,7 +1785,7 @@ class Phonopy:
                 is_gamma_center=is_gamma_center,
             )
 
-    def get_mesh_dict(self):
+    def get_mesh_dict(self) -> dict:
         """Return phonon properties calculated by mesh sampling.
 
         Returns
@@ -1837,11 +1873,11 @@ class Phonopy:
             self._mesh.grid_mapping_table,
         )
 
-    def write_hdf5_mesh(self):
+    def write_hdf5_mesh(self) -> None:
         """Write mesh calculation results in hdf5 format."""
         self._mesh.write_hdf5()
 
-    def write_yaml_mesh(self):
+    def write_yaml_mesh(self) -> None:
         """Write mesh calculation results in yaml format."""
         self._mesh.write_yaml()
 
@@ -1956,7 +1992,7 @@ class Phonopy:
         with_group_velocities=False,
         with_dynamical_matrices=False,
         nac_q_direction=None,
-    ):
+    ) -> None:
         """Run phonon calculation at specified q-points.
 
         Parameters
@@ -2071,11 +2107,11 @@ class Phonopy:
         qpt = self.get_qpoints_dict()
         return (qpt["frequencies"], qpt["eigenvectors"])
 
-    def write_hdf5_qpoints_phonon(self):
+    def write_hdf5_qpoints_phonon(self) -> None:
         """Write phonon properties calculated at q-points in hdf5 format."""
         self._qpoints.write_hdf5()
 
-    def write_yaml_qpoints_phonon(self):
+    def write_yaml_qpoints_phonon(self) -> None:
         """Write phonon properties calculated at q-points in yaml format."""
         self._qpoints.write_yaml()
 
@@ -2087,7 +2123,7 @@ class Phonopy:
         freq_max=None,
         freq_pitch=None,
         use_tetrahedron_method=True,
-    ):
+    ) -> None:
         """Run total DOS calculation.
 
         Parameters
@@ -2148,7 +2184,7 @@ class Phonopy:
         with_tight_frequency_range=False,
         write_dat=False,
         filename="total_dos.dat",
-    ):
+    ) -> None:
         """Conveniently calculate and draw total DOS."""
         self.run_mesh(
             mesh=mesh,
@@ -2166,7 +2202,7 @@ class Phonopy:
                 with_tight_frequency_range=with_tight_frequency_range,
             )
 
-    def get_total_dos_dict(self):
+    def get_total_dos_dict(self) -> dict:
         """Return total DOS.
 
         Returns
@@ -2208,13 +2244,13 @@ class Phonopy:
 
         return dos["frequency_points"], dos["total_dos"]
 
-    def set_Debye_frequency(self, freq_max_fit=None):
+    def set_Debye_frequency(self, freq_max_fit=None) -> None:
         """Calculate Debye frequency on top of total DOS."""
         self._total_dos.set_Debye_frequency(
             len(self._primitive), freq_max_fit=freq_max_fit
         )
 
-    def get_Debye_frequency(self):
+    def get_Debye_frequency(self) -> float:
         """Return Debye frequency."""
         return self._total_dos.get_Debye_frequency()
 
@@ -2266,7 +2302,7 @@ class Phonopy:
         )
         self.write_total_dos(filename=filename)
 
-    def write_total_dos(self, filename="total_dos.dat"):
+    def write_total_dos(self, filename="total_dos.dat") -> None:
         """Write total DOS to text file."""
         self._total_dos.write(filename=filename)
 
@@ -2280,7 +2316,7 @@ class Phonopy:
         use_tetrahedron_method=True,
         direction=None,
         xyz_projection=False,
-    ):
+    ) -> None:
         """Run projected DOS calculation.
 
         Parameters
@@ -2373,7 +2409,7 @@ class Phonopy:
         with_tight_frequency_range=False,
         write_dat=False,
         filename="projected_dos.dat",
-    ):
+    ) -> None:
         """Conveniently calculate and draw projected DOS.
 
         Parameters
@@ -2418,7 +2454,7 @@ class Phonopy:
                 with_tight_frequency_range=with_tight_frequency_range,
             )
 
-    def get_projected_dos_dict(self):
+    def get_projected_dos_dict(self) -> dict:
         """Return projected DOS.
 
         Projection is done to atoms and may be also done along directions
@@ -2549,7 +2585,7 @@ class Phonopy:
         )
         self.write_projected_dos(filename=filename)
 
-    def write_projected_dos(self, filename="projected_dos.dat"):
+    def write_projected_dos(self, filename="projected_dos.dat") -> None:
         """Write projected DOS to text file."""
         self._pdos.write(filename=filename)
 
@@ -2564,7 +2600,7 @@ class Phonopy:
         pretend_real=False,
         band_indices=None,
         is_projection=False,
-    ):
+    ) -> None:
         """Run calculation of thermal properties at constant volume.
 
         In phonopy, imaginary frequencies are represented as negative real
@@ -2643,7 +2679,7 @@ class Phonopy:
             pretend_real=pretend_real,
         )
 
-    def get_thermal_properties_dict(self):
+    def get_thermal_properties_dict(self) -> dict:
         """Return thermal properties.
 
         Returns
@@ -2741,7 +2777,7 @@ class Phonopy:
 
         return plt
 
-    def write_yaml_thermal_properties(self, filename="thermal_properties.yaml"):
+    def write_yaml_thermal_properties(self, filename="thermal_properties.yaml") -> None:
         """Write thermal properties in yaml format."""
         self._thermal_properties.write_yaml(filename=filename)
 
@@ -2755,7 +2791,7 @@ class Phonopy:
         direction=None,
         freq_min=None,
         freq_max=None,
-    ):
+    ) -> None:
         """Run thermal displacements calculation.
 
         Parameters
@@ -2835,7 +2871,7 @@ class Phonopy:
             freq_max=freq_max,
         )
 
-    def get_thermal_displacements_dict(self):
+    def get_thermal_displacements_dict(self) -> dict:
         """Return thermal displacements."""
         if self._thermal_displacements is None:
             msg = "run_thermal_displacements has to be done."
@@ -2874,7 +2910,7 @@ class Phonopy:
 
         return plt
 
-    def write_yaml_thermal_displacements(self):
+    def write_yaml_thermal_displacements(self) -> None:
         """Write thermal displacements in yaml format."""
         self._thermal_displacements.write_yaml()
 
@@ -2887,7 +2923,7 @@ class Phonopy:
         temperatures=None,
         freq_min=None,
         freq_max=None,
-    ):
+    ) -> None:
         """Run thermal displacement matrices calculation.
 
         Parameters
@@ -2958,7 +2994,7 @@ class Phonopy:
             freq_max=freq_max,
         )
 
-    def get_thermal_displacement_matrices_dict(self):
+    def get_thermal_displacement_matrices_dict(self) -> dict:
         """Return thermal displacement matrices."""
         if self._thermal_displacement_matrices is None:
             msg = "run_thermal_displacement_matrices has to be done."
@@ -2982,11 +3018,11 @@ class Phonopy:
         tdm = self.get_thermal_displacement_matrices_dict()
         return (tdm["temperatures"], tdm["thermal_displacement_matrices"])
 
-    def write_yaml_thermal_displacement_matrices(self):
+    def write_yaml_thermal_displacement_matrices(self) -> None:
         """Write thermal displacement matrices in yaml format."""
         self._thermal_displacement_matrices.write_yaml()
 
-    def write_thermal_displacement_matrix_to_cif(self, temperature_index):
+    def write_thermal_displacement_matrix_to_cif(self, temperature_index) -> None:
         """Write thermal displacement matrices at a termperature in cif."""
         self._thermal_displacement_matrices.write_cif(
             self._primitive, temperature_index
@@ -3001,7 +3037,7 @@ class Phonopy:
         num_div=None,
         shift=None,
         filename=None,
-    ):
+    ) -> str:
         """Write atomic modulations in animation format.
 
         Returns
@@ -3060,7 +3096,7 @@ class Phonopy:
         delta_q=None,
         derivative_order=None,
         nac_q_direction=None,
-    ):
+    ) -> None:
         """Generate atomic displacements of phonon modes.
 
         The design of this feature is not very satisfactory, and thus API.
@@ -3103,7 +3139,7 @@ class Phonopy:
         )
         self._modulation.run()
 
-    def get_modulated_supercells(self):
+    def get_modulated_supercells(self) -> list[PhonopyAtoms]:
         """Return cells with atom modulations.
 
         list of PhonopyAtoms
@@ -3112,7 +3148,7 @@ class Phonopy:
         """
         return self._modulation.get_modulated_supercells()
 
-    def get_modulations_and_supercell(self):
+    def get_modulations_and_supercell(self) -> tuple[np.ndarray, PhonopyAtoms]:
         """Return atomic modulations and perfect supercell.
 
         (modulations, supercell)
@@ -3123,14 +3159,14 @@ class Phonopy:
         """
         return self._modulation.get_modulations_and_supercell()
 
-    def write_modulations(self, calculator=None, optional_structure_info=None):
+    def write_modulations(self, calculator=None, optional_structure_info=None) -> None:
         """Write modulated structures to MPOSCAR's."""
         self._modulation.write(
             interface_mode=calculator,
             optional_structure_info=optional_structure_info,
         )
 
-    def write_yaml_modulations(self):
+    def write_yaml_modulations(self) -> None:
         """Write atomic modulations in yaml format."""
         self._modulation.write_yaml()
 
@@ -3141,7 +3177,7 @@ class Phonopy:
         is_little_cogroup=False,
         nac_q_direction=None,
         degeneracy_tolerance=1e-4,
-    ):
+    ) -> None:
         """Identify ir-reps of phonon modes.
 
         The design of this API is not very satisfactory and is expceted
@@ -3175,13 +3211,17 @@ class Phonopy:
 
     def get_irreps(self):
         """Return Ir-reps."""
+        warnings.warn(
+            "Phonopy.get_irreps() is deprecated. " "Use Phonopy.irreps attribute.",
+            DeprecationWarning,
+        )
         return self._irreps
 
-    def show_irreps(self, show_irreps=False):
+    def show_irreps(self, show_irreps=False) -> None:
         """Show Ir-reps."""
         self._irreps.show(show_irreps=show_irreps)
 
-    def write_yaml_irreps(self, show_irreps=False):
+    def write_yaml_irreps(self, show_irreps=False) -> None:
         """Write Ir-reps in yaml format."""
         self._irreps.write_yaml(show_irreps=show_irreps)
 
@@ -3208,7 +3248,7 @@ class Phonopy:
         )
         return self._group_velocity.get_group_velocity()
 
-    def get_group_velocity_at_q(self, q_point):
+    def get_group_velocity_at_q(self, q_point) -> np.ndarray:
         """Return group velocity at a q-point."""
         if self._group_velocity is None:
             self._set_group_velocity()
@@ -3225,7 +3265,9 @@ class Phonopy:
         return self._band_structure.group_velocities
 
     # Moment
-    def run_moment(self, order=1, is_projection=False, freq_min=None, freq_max=None):
+    def run_moment(
+        self, order=1, is_projection=False, freq_min=None, freq_max=None
+    ) -> None:
         """Run moment calculation."""
         if self._mesh is None:
             msg = "run_mesh has to be done before run_moment."
@@ -3262,7 +3304,7 @@ class Phonopy:
             freq_max=freq_max,
         )
 
-    def get_moment(self):
+    def get_moment(self) -> Optional[float]:
         """Return moment."""
         return self._moment.moment
 
@@ -3274,7 +3316,7 @@ class Phonopy:
         scattering_lengths=None,
         freq_min=None,
         freq_max=None,
-    ):
+    ) -> None:
         """Initialize dynamic structure factor calculation.
 
         *******************************************************************
@@ -3348,7 +3390,7 @@ class Phonopy:
         scattering_lengths=None,
         freq_min=None,
         freq_max=None,
-    ):
+    ) -> None:
         """Run dynamic structure factor calculation.
 
         See the detail of parameters at
@@ -3400,7 +3442,7 @@ class Phonopy:
                 freq_max=freq_max,
             )
 
-    def get_dynamic_structure_factor(self):
+    def get_dynamic_structure_factor(self) -> tuple[np.ndarray, np.ndarray]:
         """Return dynamic structure factors."""
         return (
             self._dynamic_structure_factor.qpoints,
@@ -3412,7 +3454,7 @@ class Phonopy:
         dist_func: Optional[str] = None,
         cutoff_frequency: Optional[float] = None,
         max_distance: Optional[float] = None,
-    ):
+    ) -> None:
         """Initialize random displacements at finite temperature.
 
         dist_func : str or None, optional
@@ -3448,7 +3490,7 @@ class Phonopy:
         number_of_snapshots: int,
         is_plusminus: bool = False,
         random_seed: Optional[int] = None,
-    ):
+    ) -> np.ndarray:
         """Generate random displacements from phonon structure.
 
         Some more details are written at generate_displacements.
@@ -3485,7 +3527,9 @@ class Phonopy:
             )
         return d
 
-    def save(self, filename="phonopy_params.yaml", settings=None, hdf5_settings=None):
+    def save(
+        self, filename="phonopy_params.yaml", settings=None, hdf5_settings=None
+    ) -> None:
         """Save phonopy parameters into file.
 
         Parameters
@@ -3532,7 +3576,7 @@ class Phonopy:
         with open(filename, "w") as w:
             w.write(str(phpy_yaml))
 
-    def ph2ph(self, supercell_matrix, with_nac=False):
+    def ph2ph(self, supercell_matrix, with_nac=False) -> Phonopy:
         """Transform force constants in Phonopy class instance to other shape.
 
         Fourier interpolation of force constants is performed. This Phonopy
@@ -3636,7 +3680,7 @@ class Phonopy:
         fc_calculator=None,
         fc_calculator_options=None,
         decimals=None,
-    ):
+    ) -> None:
         if self._displacement_dataset is not None:
             if fc_calculator is not None:
                 disps, forces = get_displacements_and_forces(self._displacement_dataset)
@@ -3669,7 +3713,7 @@ class Phonopy:
                     decimals=decimals,
                 )
 
-    def _set_dynamical_matrix(self):
+    def _set_dynamical_matrix(self) -> None:
         import phonopy._phonopy as phonoc
 
         self._dynamical_matrix = None
@@ -3710,7 +3754,7 @@ class Phonopy:
         if self._group_velocity is not None:
             self._set_group_velocity()
 
-    def _set_group_velocity(self):
+    def _set_group_velocity(self) -> None:
         if self._dynamical_matrix is None:
             raise RuntimeError("Dynamical matrix has not yet built.")
 
@@ -3738,7 +3782,7 @@ class Phonopy:
             frequency_factor_to_THz=self._factor,
         )
 
-    def _search_symmetry(self):
+    def _search_symmetry(self) -> None:
         self._symmetry = Symmetry(
             self._supercell,
             self._symprec,
@@ -3746,7 +3790,7 @@ class Phonopy:
             s2p_map=self._primitive.s2p_map,
         )
 
-    def _search_primitive_symmetry(self):
+    def _search_primitive_symmetry(self) -> None:
         self._primitive_symmetry = Symmetry(
             self._primitive, self._symprec, self._is_symmetry
         )
@@ -3759,7 +3803,7 @@ class Phonopy:
                 "cell are different."
             )
 
-    def _build_supercell(self):
+    def _build_supercell(self) -> None:
         self._supercell = get_supercell(
             self._unitcell,
             self._supercell_matrix,
@@ -3767,7 +3811,7 @@ class Phonopy:
             symprec=self._symprec,
         )
 
-    def _build_supercells_with_displacements(self):
+    def _build_supercells_with_displacements(self) -> None:
         all_positions = []
         if "first_atoms" in self._displacement_dataset:  # type-1
             for disp in self._displacement_dataset["first_atoms"]:
@@ -3793,7 +3837,7 @@ class Phonopy:
             )
         self._supercells_with_displacements = supercells
 
-    def _build_primitive_cell(self):
+    def _build_primitive_cell(self) -> None:
         """Create primitive cell.
 
         primitive_matrix:
@@ -3824,12 +3868,14 @@ class Phonopy:
             )
             raise RuntimeError(msg)
 
-    def _set_primitive_matrix(self, primitive_matrix):
+    def _set_primitive_matrix(
+        self, primitive_matrix
+    ) -> Optional[Union[str, np.ndarray]]:
         pmat = get_primitive_matrix(primitive_matrix, symprec=self._symprec)
         if isinstance(pmat, str) and pmat == "auto":
             return guess_primitive_matrix(self._unitcell, symprec=self._symprec)
         else:
             return pmat
 
-    def _shape_supercell_matrix(self, smat):
+    def _shape_supercell_matrix(self, smat) -> np.ndarray:
         return shape_supercell_matrix(smat)

phonopy/api_qha.py

@@ -1,4 +1,5 @@
 """API for QHA calculation."""
+
 # Copyright (C) 2015 Atsushi Togo
 # All rights reserved.
 #

phonopy/cui/__init__.py

@@ -1,4 +1,5 @@
 """Routines for command user interface."""
+
 # Copyright (C) 2021 Atsushi Togo
 # All rights reserved.
 #

phonopy/cui/collect_cell_info.py

@@ -1,4 +1,5 @@
 """Routines to collect crystal structure information."""
+
 # Copyright (C) 2018 Atsushi Togo
 # All rights reserved.
 #

phonopy/cui/create_force_sets.py

@@ -1,4 +1,5 @@
 """Utilities to create force sets for main CUI script."""
+
 # Copyright (C) 2020 Atsushi Togo
 # All rights reserved.
 #

phonopy/cui/load.py

@@ -1,4 +1,5 @@
 """Phonopy loader."""
+
 # Copyright (C) 2018 Atsushi Togo
 # All rights reserved.
 #

phonopy/cui/load_helper.py

@@ -1,4 +1,5 @@
 """Helper methods of phonopy loader."""
+
 # Copyright (C) 2018 Atsushi Togo
 # All rights reserved.
 #

phonopy/cui/phonopy_argparse.py

@@ -1,4 +1,5 @@
 """Phonopy command line argument parser."""
+
 # Copyright (C) 2016 Atsushi Togo
 # All rights reserved.
 #
@@ -65,9 +66,7 @@ def show_deprecated_option_warnings(deprecated):
     print("")
 
 
-def get_parser(
-    fc_symmetry=False, is_nac=False, include_born=False, load_phonopy_yaml=False
-):
+def get_parser(fc_symmetry=False, is_nac=False, load_phonopy_yaml=False):
     """Return ArgumentParser instance."""
     deprecated = fix_deprecated_option_names(sys.argv)
     import argparse

phonopy/cui/phonopy_script.py

@@ -1,4 +1,5 @@
 """Phonopy command user interface."""
+
 # Copyright (C) 2020 Atsushi Togo
 # All rights reserved.
 #
@@ -873,7 +874,7 @@ def store_nac_params(
                 print("-" * 76)
 
 
-def run(phonon: Phonopy, settings, cell_info, plot_conf, log_level):
+def run_calculation(phonon: Phonopy, settings, plot_conf, log_level):
     """Run phonon calculations."""
     interface_mode = phonon.calculator
     physical_units = get_default_physical_units(interface_mode)
@@ -1491,15 +1492,7 @@ def start_phonopy(**argparse_control):
         if deprecated:
             show_deprecated_option_warnings(deprecated)
 
-    return (
-        args,
-        log_level,
-        {
-            "plot_graph": args.is_graph_plot,
-            "save_graph": args.is_graph_save,
-            "with_legend": args.is_legend,
-        },
-    )
+    return args, log_level
 
 
 def read_phonopy_settings(args, argparse_control, log_level):
@@ -1760,7 +1753,17 @@ def main(**argparse_control):
     ############################################
     load_phonopy_yaml = argparse_control.get("load_phonopy_yaml", False)
 
-    args, log_level, plot_conf = start_phonopy(**argparse_control)
+    if "args" in argparse_control:  # For pytest
+        args = argparse_control["args"]
+        log_level = args.log_level
+    else:
+        args, log_level = start_phonopy(**argparse_control)
+
+    plot_conf = {
+        "plot_graph": args.is_graph_plot,
+        "save_graph": args.is_graph_save,
+        "with_legend": args.is_legend,
+    }
 
     settings, confs, cell_filename = read_phonopy_settings(
         args, argparse_control, log_level
@@ -1988,7 +1991,7 @@ def main(**argparse_control):
             print(" - %s" % mode)
         print("-" * 76)
 
-    run(phonon, settings, cell_info, plot_conf, log_level)
+    run_calculation(phonon, settings, plot_conf, log_level)
 
     ########################
     # Phonopy finalization #

phonopy/cui/settings.py

@@ -1,4 +1,5 @@
 """Phonopy input and command option tools."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #
@@ -1639,9 +1640,9 @@ class PhonopyConfParser(ConfParser):
 
         if "rd_temperature" in arg_list:
             if self._args.rd_temperature is not None:
-                self._confs[
-                    "random_displacement_temperature"
-                ] = self._args.rd_temperature
+                self._confs["random_displacement_temperature"] = (
+                    self._args.rd_temperature
+                )
 
         if "temperature" in arg_list:
             if self._args.temperature is not None:
@@ -1660,11 +1661,7 @@ class PhonopyConfParser(ConfParser):
         if "random_seed" in arg_list:
             if self._args.random_seed:
                 seed = self._args.random_seed
-                if (
-                    np.issubdtype(type(seed), np.integer)
-                    and seed >= 0
-                    and seed < 2**32
-                ):
+                if np.issubdtype(type(seed), np.integer) and seed >= 0 and seed < 2**32:
                     self._confs["random_seed"] = seed
 
         if "include_fc" in arg_list:

phonopy/cui/show_symmetry.py

@@ -1,4 +1,5 @@
 """Show symmetry information invoked by --symmetry command option."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #
@@ -35,7 +36,6 @@
 
 import numpy as np
 import spglib
-from spglib import get_pointgroup
 
 from phonopy import Phonopy
 from phonopy.interface.calculator import (
@@ -56,9 +56,9 @@ def check_symmetry(phonon: Phonopy, optional_structure_info):
 
     if phonon.unitcell.magnetic_moments is None:
         base_fname = get_default_cell_filename(phonon.calculator)
-        symprec = phonon.primitive_symmetry.get_symmetry_tolerance()
+        symprec = phonon.primitive_symmetry.tolerance
         (bravais_lattice, bravais_pos, bravais_numbers) = spglib.refine_cell(
-            phonon.primitive, symprec
+            phonon.primitive.totuple(), symprec
         )
         bravais = PhonopyAtoms(
             numbers=bravais_numbers, scaled_positions=bravais_pos, cell=bravais_lattice
@@ -115,7 +115,7 @@ def _get_symmetry_yaml(cell: PhonopyAtoms, symmetry: Symmetry, phonopy_version=N
         spg_number = int(spg_number.replace("(", "").replace(")", ""))
         lines.append("space_group_type: '%s'" % spg_symbol)
         lines.append("space_group_number: %d" % spg_number)
-        lines.append("point_group_type: '%s'" % symmetry.get_pointgroup())
+        lines.append("point_group_type: '%s'" % symmetry.pointgroup_symbol)
     lines.append("space_group_operations:")
     for i, (r, t) in enumerate(zip(rotations, translations)):
         lines.append("- rotation: # %d" % (i + 1))
@@ -142,7 +142,7 @@ def _get_symmetry_yaml(cell: PhonopyAtoms, symmetry: Symmetry, phonopy_version=N
 
         if cell.magnetic_moments is None:
             lines.append("  Wyckoff: '%s'" % wyckoffs[i])
-        site_pointgroup = get_pointgroup(sitesym)
+        site_pointgroup = spglib.get_pointgroup(sitesym)
         lines.append("  site_point_group: '%s'" % site_pointgroup[0].strip())
         lines.append("  orientation:")
         for v in site_pointgroup[2]:

phonopy/exception.py

@@ -1,4 +1,5 @@
 """Phonopy exceptions."""
+
 # Copyright (C) 2022 Atsushi Togo
 # All rights reserved.
 #

phonopy/file_IO.py

@@ -1,4 +1,5 @@
 """File I/O related routines."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/gruneisen/__init__.py

@@ -1,4 +1,5 @@
 """Mode Grueneisen paramater calculation."""
+
 # Copyright (C) 2017 Atsushi Togo
 # All rights reserved.
 #

phonopy/gruneisen/band_structure.py

@@ -1,4 +1,5 @@
 """Mode Grueneisen parameter band structure calculation."""
+
 # Copyright (C) 2012 Atsushi Togo
 # All rights reserved.
 #

phonopy/gruneisen/core.py

@@ -1,4 +1,5 @@
 """Mode Grueneisen parameter calculation."""
+
 # Copyright (C) 2012 Atsushi Togo
 # All rights reserved.
 #

phonopy/gruneisen/mesh.py

@@ -1,4 +1,5 @@
 """Mode Grueneisen parameters calculation on sampling mesh."""
+
 # Copyright (C) 2012 Atsushi Togo
 # All rights reserved.
 #

phonopy/harmonic/__init__.py

@@ -1,4 +1,5 @@
 """Routines of pre-phonon calculations."""
+
 # Copyright (C) 2021 Atsushi Togo
 # All rights reserved.
 #

phonopy/harmonic/derivative_dynmat.py

@@ -1,4 +1,5 @@
 """Calculation of derivative of dynamical matrix with respect to q."""
+
 # Copyright (C) 2013 Atsushi Togo
 # All rights reserved.
 #

phonopy/harmonic/displacement.py

@@ -1,4 +1,5 @@
 """Routines to handle displacements in supercells."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #
@@ -245,16 +246,70 @@ def _determinant(a, b, c):
 
 def get_random_displacements_dataset(
     num_supercells: int,
-    distance: float,
     num_atoms: int,
+    distance: float,
     random_seed: Optional[int] = None,
     is_plusminus: bool = False,
+    is_random_distance: bool = False,
+    min_distance: Optional[float] = None,
 ) -> np.ndarray:
-    """Return random displacements at constant displacement distance."""
-    disps = (
-        _get_random_directions(num_atoms * num_supercells, random_seed=random_seed)
-        * distance
-    )
+    """Return random displacements at constant displacement distance.
+
+    num_supercells : int
+        Number of snapshots of supercells with random displacements. Random
+        displacements are generated displacing all atoms in random directions
+        with a fixed displacement distance specified by 'distance' parameter,
+        i.e., all atoms in supercell are displaced with the same displacement
+        distance in direct space.
+    num_atoms : int
+        Number of atoms in supercell.
+    distance : float
+        Displacement distance. Unit is the same as that used for crystal
+        structure.
+    random_seed : int or None, optional
+        Random seed for random displacements generation. Default is None.
+    is_plusminus : True, or False, optional
+        In addition to sets of usual random displacements for supercell, sets
+        of the opposite displacements for supercell are concatenated.
+        Therefore, total number of sets of displacements is `2 *
+        num_supercells`. Default is False.
+    is_random_distance : bool, optional
+        Random direction displacements are generated also with random
+        amplitudes. The maximum value is defined by `distance` and minimum value
+        is given by `min_distance`. Default is False. Random distance is given
+        by `sqrt(random(distance - min_distance) + min_distance)`.
+    min_distance : float or None, optional
+        In random direction displacements generation with random distance
+        (`is_random_distance=True`), the minimum distance is given by this
+        value.
+
+    """
+    if is_random_distance:
+        if min_distance is None:
+            _min_distance = 0.0
+        else:
+            _min_distance = min_distance
+
+    if np.issubdtype(type(random_seed), np.integer):
+        rng = np.random.default_rng(seed=random_seed)
+    else:
+        rng = np.random.default_rng()
+
+    if is_random_distance:
+        if distance < _min_distance:
+            raise RuntimeError(
+                "Random displacements generation failed. min_distance is too large."
+            )
+        directions = _get_random_directions(num_atoms * num_supercells, rng)
+        rand_dists = np.array([])
+        while len(rand_dists) < num_atoms * num_supercells:
+            rd = np.sqrt(rng.random(num_atoms * num_supercells)) * distance
+            rand_dists = np.r_[rand_dists, rd[rd > _min_distance]]
+        disps = rand_dists[: num_atoms * num_supercells, None] * directions
+    else:
+        directions = _get_random_directions(num_atoms * num_supercells, rng)
+        disps = directions * distance
+
     supercell_disps = np.array(
         disps.reshape(num_supercells, num_atoms, 3), dtype="double", order="C"
     )
@@ -267,15 +322,9 @@ def get_random_displacements_dataset(
     return supercell_disps
 
 
-def _get_random_directions(num_atoms, random_seed=None):
+def _get_random_directions(num_atoms: int, rng: np.random.Generator) -> np.ndarray:
     """Return random directions in sphere with radius 1."""
-    if (
-        np.issubdtype(type(random_seed), np.integer)
-        and random_seed >= 0
-        and random_seed < 2**32
-    ):
-        np.random.seed(random_seed)
-
-    xy = np.random.randn(3, num_atoms)
-    r = np.sqrt((xy**2).sum(axis=0))
-    return (xy / r).T
+    xy = rng.standard_normal(size=(3, num_atoms * 2))
+    r = np.linalg.norm(xy, axis=0)
+    condition = r > 1e-10
+    return (xy[:, condition][:, :num_atoms] / r[condition][:num_atoms]).T

phonopy/harmonic/dynamical_matrix.py

@@ -37,7 +37,7 @@
 import itertools
 import sys
 import warnings
-from typing import Type, Union
+from typing import Optional, Type, Union
 
 import numpy as np
 
@@ -205,7 +205,7 @@ class DynamicalMatrix:
         return self.force_constants
 
     @property
-    def dynamical_matrix(self):
+    def dynamical_matrix(self) -> Optional[np.ndarray]:
         """Return dynamcial matrix calculated at q.
 
         Returns
@@ -230,7 +230,7 @@ class DynamicalMatrix:
         warnings.warn(
             "DynamicalMatrix.get_get_dynamical_matrix() is "
             "deprecated."
-            "Use DynamicalMatrix.get_dynamical_matrix attribute.",
+            "Use DynamicalMatrix.dynamical_matrix attribute.",
             DeprecationWarning,
         )
         return self.dynamical_matrix

phonopy/harmonic/dynmat_to_fc.py

@@ -1,4 +1,5 @@
 """Transform dynamical matrix to force constants."""
+
 # Copyright (C) 2014 Atsushi Togo
 # All rights reserved.
 #

phonopy/harmonic/force_constants.py

@@ -1,4 +1,5 @@
 """Force constants calculation."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/__init__.py

@@ -1,4 +1,5 @@
 """Routines for user and calculator interfaces to phonopy."""
+
 # Copyright (C) 2019 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/abacus.py

@@ -1,4 +1,5 @@
 """ABACUS calculator interface."""
+
 # Copyright (C) 2022 Yuyang Ji
 # All rights reserved.
 #
@@ -313,11 +314,24 @@ def read_abacus_output(filename):
                 natom = int(re.search("[0-9]+", line).group())
             if re.search(r"TOTAL-FORCE \(eV/Angstrom\)", line):
                 force = np.zeros((natom, 3))
-                for i in range(4):
-                    file.readline()
-                for i in range(natom):
-                    _, fx, fy, fz = file.readline().split()
-                    force[i] = (float(fx), float(fy), float(fz))
+                _match_pattern = r"^(\s*)([A-Za-z]*[0-9]+)((\s*[+-]?"
+                _match_pattern += r"[0-9]+\.[0-9]+(e[+-][0-9]{2})?){3})(.*)$"
+                _match = re.match(_match_pattern, line)
+                while not _match:
+                    line = file.readline()
+                    _match = re.match(_match_pattern, line)
+                iatom = 0
+                while _match:
+                    print(_match.group(3).split())
+                    fx, fy, fz = (_match.group(3).split()[i].strip() for i in range(3))
+                    force[iatom] = (float(fx), float(fy), float(fz))
+                    iatom += 1
+                    if iatom == natom:
+                        break
+                    else:
+                        line = file.readline()
+                        _match = re.match(_match_pattern, line)
+
     if force is None:
         raise ValueError("Force data not found.")
 

phonopy/interface/abinit.py

@@ -1,4 +1,5 @@
 """Abinit calculator interface."""
+
 # Copyright (C) 2014 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/aims.py

@@ -1,4 +1,5 @@
 """FHIaims calculator interface."""
+
 # FHIaims.py - IO routines for phonopy-FHI-aims
 # methods compatible with the corresponding ones from ase.io.aims
 # only minimal subset of functionality required within phonopy context is implemented

phonopy/interface/alm.py

@@ -1,4 +1,5 @@
 """ALM force constants interface."""
+
 # Copyright (C) 2018 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/calculator.py

@@ -1,4 +1,5 @@
 """Routines to handle various calculator interfaces."""
+
 # Copyright (C) 2014 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/castep.py

@@ -1,4 +1,5 @@
 """CASTEP calculator interface."""
+
 # Copyright (C) 2014 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/cif.py

@@ -1,4 +1,5 @@
 """Tests for CIF tools."""
+
 # Copyright (C) 2016 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/cp2k.py

@@ -1,4 +1,5 @@
 """CP2K calculator interface."""
+
 # vim: set fileencoding=utf-8 :
 # Copyright (C) 2017-2019 Tiziano Müller
 # All rights reserved.

phonopy/interface/crystal.py

@@ -1,4 +1,5 @@
 """CRYSTAL calculator interface."""
+
 # Copyright (C) 2016 Antti J. Karttunen (antti.j.karttunen@iki.fi)
 # All rights reserved.
 #

phonopy/interface/dftbp.py

@@ -1,4 +1,5 @@
 """DFTB+ calculator interface."""
+
 # Copyright (C) 2015 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/elk.py

@@ -1,4 +1,5 @@
 """Elk calculator interface."""
+
 # Copyright (C) 2015 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/fc_calculator.py

@@ -1,4 +1,5 @@
 """Force constants calculator interface."""
+
 # Copyright (C) 2019 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/fleur.py

@@ -1,4 +1,5 @@
 """Fleur calculator interface."""
+
 # Copyright (C) 2021 Alexander Neukirchen
 # All rights reserved.
 #

phonopy/interface/hiphive_interface.py

@@ -1,4 +1,5 @@
 """hiPhive force constants calculator interface."""
+
 # Copyright (C) 2018 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/lammps.py

@@ -1,4 +1,5 @@
 """LAMMPS calculator interface."""
+
 # Copyright (C) 2023 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/phonopy_yaml.py

@@ -1,4 +1,5 @@
 """phonopy.yaml reader and writer."""
+
 # Copyright (C) 2018 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/qe.py

@@ -1,4 +1,5 @@
 """QE calculator interface."""
+
 # Copyright (C) 2014 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/siesta.py

@@ -1,4 +1,5 @@
 """SIESTA calculator interface."""
+
 # Copyright (C) 2015 Henrique Pereira Coutada Miranda
 # All rights reserved.
 #

phonopy/interface/turbomole.py

@@ -1,4 +1,5 @@
 """CRYSTAL calculator interface."""
+
 # Copyright (C) 2019 Antti J. Karttunen (antti.j.karttunen@iki.fi)
 # All rights reserved.
 #

phonopy/interface/vasp.py

@@ -1,4 +1,5 @@
 """VASP calculator interface."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/interface/wien2k.py

@@ -1,4 +1,5 @@
 """Wien2k calculator interface."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/__init__.py

@@ -1,4 +1,5 @@
 """Routines of phonon and post-phonon calculations."""
+
 # Copyright (C) 2021 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/animation.py

@@ -1,4 +1,5 @@
 """Convert phonon results to animation formats."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/character_table.py

@@ -1,4 +1,5 @@
 """Character table datasets of point groups."""
+
 # Copyright (C) 2020 Atsushi Togo
 # All rights reserved.
 #
@@ -464,7 +465,23 @@ character_table = {
                     ((1, 0, 0), (1, -1, 0), (0, 0, 1)),
                 ),
             },
-        }
+        },
+        {
+            "rotation_list": ("E", "C3", "sgv"),
+            "character_table": {"A1": (1, 1, 1), "A2": (1, 1, -1), "E": (2, -1, 0)},
+            "mapping_table": {
+                "E": (((1, 0, 0), (0, 1, 0), (0, 0, 1)),),
+                "C3": (
+                    ((0, -1, 0), (1, -1, 0), (0, 0, 1)),
+                    ((-1, 1, 0), (-1, 0, 0), (0, 0, 1)),
+                ),
+                "sgv": (
+                    ((0, 1, 0), (1, 0, 0), (0, 0, 1)),
+                    ((1, -1, 0), (0, -1, 0), (0, 0, 1)),
+                    ((-1, 0, 0), (-1, 1, 0), (0, 0, 1)),
+                ),
+            },
+        },
     ],
     # D3d (20)
     "-3m": [

phonopy/phonon/degeneracy.py

@@ -1,4 +1,5 @@
 """Utility routines to handle degeneracy."""
+
 # Copyright (C) 2014 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/dos.py

@@ -1,4 +1,5 @@
 """Calculation of density of states."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/group_velocity.py

@@ -1,4 +1,5 @@
 """Group velocity calculation."""
+
 # Copyright (C) 2013 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/irreps.py

@@ -1,4 +1,5 @@
 """Calculate irreducible representation from eigenvectors."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/mesh.py

@@ -1,4 +1,5 @@
 """Phonon calculation on sampling mesh."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/modulation.py

@@ -1,4 +1,5 @@
 """Create atomic displacements."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/moment.py

@@ -1,4 +1,5 @@
 """Calculate phonon state moments."""
+
 # Copyright (C) 2016 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/qpoints.py

@@ -1,4 +1,5 @@
 """Phonon calculation at specific q-points."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/random_displacements.py

@@ -1,4 +1,5 @@
 """Calculate random displacements from phonons at finite temperatures."""
+
 # Copyright (C) 2018 Atsushi Togo
 # All rights reserved.
 #
@@ -242,33 +243,42 @@ class RandomDisplacements:
         number_of_snapshots : int
             Number of snapshots to be generated.
         random_seed : int or None, optional
-            Random seed passed to np.random.seed. Default is None. Integer
-            number has to be positive.
+            Random seed passed to np.random.default_rng(seed). Default is None.
         randn : tuple
-            (randn_ii, randn_ij).
-            Used for testing purpose for the fixed random numbers of
-            np.random.normal that can depends on system.
+            (randn_ii, randn_ij). Used for testing purpose for the fixed random
+            numbers of random.Generator.standard_normal that can depends on
+            system.
 
         """
-        np.random.seed(seed=random_seed)
+        if np.issubdtype(type(random_seed), np.integer):
+            rng = np.random.default_rng(seed=random_seed)
+        else:
+            rng = np.random.default_rng()
 
         N = len(self._comm_points)
 
         # This randn is used only for testing purpose.
         if randn is None:
-            randn_ii = None
-            randn_ij = None
+            shape = (
+                len(self._eigvals_ii),
+                number_of_snapshots,
+                len(self._eigvals_ii[0]),
+            )
+            randn_ii = rng.standard_normal(size=shape)
+            shape = (
+                len(self._eigvals_ij),
+                2,
+                number_of_snapshots,
+                len(self._eigvals_ij[0]),
+            )
+            randn_ij = rng.standard_normal(size=shape)
         else:
             randn_ii = randn[0]
             randn_ij = randn[1]
 
-        u_ii, self._conditions_ii = self._solve_ii(
-            T, number_of_snapshots, randn=randn_ii
-        )
+        u_ii, self._conditions_ii = self._solve_ii(T, number_of_snapshots, randn_ii)
         if self._ij:
-            u_ij, self._conditions_ij = self._solve_ij(
-                T, number_of_snapshots, randn=randn_ij
-            )
+            u_ij, self._conditions_ij = self._solve_ij(T, number_of_snapshots, randn_ij)
         else:
             u_ij = 0
             self._conditions_ij = None
@@ -481,7 +491,12 @@ class RandomDisplacements:
         self._comm_points = get_commensurate_points_in_integers(smat)
         self._ii, self._ij = categorize_commensurate_points(self._comm_points)
 
-    def _solve_ii(self, T, number_of_snapshots, randn=None):
+    def _solve_ii(
+        self,
+        T: float,
+        number_of_snapshots: int,
+        randn: np.ndarray,
+    ):
         """Solve ii terms.
 
         randn parameter is used for the test.
@@ -490,14 +505,9 @@ class RandomDisplacements:
         natom = len(self._dynmat.supercell)
         u = np.zeros((number_of_snapshots, natom, 3), dtype="double")
 
-        shape = (len(self._eigvals_ii), number_of_snapshots, len(self._eigvals_ii[0]))
-        if randn is None:
-            _randn = np.random.normal(size=shape)
-        else:
-            _randn = randn
         sigmas, conditions = self._get_sigma(self._eigvals_ii, T)
         for norm_dist, sigma, eigvecs, phase in zip(
-            _randn, sigmas, self._eigvecs_ii, self._phase_ii
+            randn, sigmas, self._eigvecs_ii, self._phase_ii
         ):
             u_red = np.dot(norm_dist * sigma, eigvecs.T).reshape(
                 number_of_snapshots, -1, 3
@@ -508,7 +518,12 @@ class RandomDisplacements:
 
         return u, conditions
 
-    def _solve_ij(self, T, number_of_snapshots, randn=None):
+    def _solve_ij(
+        self,
+        T: float,
+        number_of_snapshots: int,
+        randn: Optional[np.ndarray] = None,
+    ):
         """Solve ij terms.
 
         randn parameter is used for the test.
@@ -516,19 +531,9 @@ class RandomDisplacements:
         """
         natom = len(self._dynmat.supercell)
         u = np.zeros((number_of_snapshots, natom, 3), dtype="double")
-        shape = (
-            len(self._eigvals_ij),
-            2,
-            number_of_snapshots,
-            len(self._eigvals_ij[0]),
-        )
-        if randn is None:
-            _randn = np.random.normal(size=shape)
-        else:
-            _randn = randn
         sigmas, conditions = self._get_sigma(self._eigvals_ij, T)
         for norm_dist, sigma, eigvecs, phase in zip(
-            _randn, sigmas, self._eigvecs_ij, self._phase_ij
+            randn, sigmas, self._eigvecs_ij, self._phase_ij
         ):
             u_red = np.dot(norm_dist * sigma, eigvecs.T).reshape(
                 2, number_of_snapshots, -1, 3

phonopy/phonon/tetrahedron_mesh.py

@@ -1,4 +1,5 @@
 """Linear tetrahedron method on regular grid."""
+
 # copyright (C) 2013 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/thermal_displacement.py

@@ -1,4 +1,5 @@
 """Calculations of thermal displacements."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/phonon/thermal_properties.py

@@ -1,4 +1,5 @@
 """Phonon thermal properties at constant volume."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/qha/__init__.py

@@ -1,4 +1,5 @@
 """Quasi harmonic approximation."""
+
 # Copyright (C) 2012 Atsushi Togo
 # All rights reserved.
 #

phonopy/qha/core.py

@@ -1,4 +1,5 @@
 """Core routines for QHA."""
+
 # Copyright (C) 2012 Atsushi Togo
 # All rights reserved.
 #
@@ -1123,9 +1124,7 @@ class QHA:
                     msg += ["At least 5 volume points are needed for the fitting."]
                 raise RuntimeError("\n".join(msg))
 
-            dsdv_t = np.dot(
-                parameters[:4], np.array([4 * x**3, 3 * x**2, 2 * x, 1])
-            )
+            dsdv_t = np.dot(parameters[:4], np.array([4 * x**3, 3 * x**2, 2 * x, 1]))
             self._volume_entropy_parameters.append(parameters)
 
             try:

phonopy/qha/electron.py

@@ -1,4 +1,5 @@
 """Calculation of free energy of one-electronic states."""
+
 # Copyright (C) 2018 Atsushi Togo
 # All rights reserved.
 #

phonopy/qha/eos.py

@@ -1,4 +1,5 @@
 """Equation of states and fitting routine."""
+
 # Copyright (C) 2012 Atsushi Togo
 # All rights reserved.
 #
@@ -80,11 +81,7 @@ def get_eos(eos):
         x = (v / p[3]) ** (1.0 / 3)
         xi = 3.0 / 2 * (p[2] - 1)
         return p[0] + (
-            9
-            * p[1]
-            * p[3]
-            / (xi**2)
-            * (1 + (xi * (1 - x) - 1) * np.exp(xi * (1 - x)))
+            9 * p[1] * p[3] / (xi**2) * (1 + (xi * (1 - x) - 1) * np.exp(xi * (1 - x)))
         )
 
     if eos == "murnaghan":
@@ -152,7 +149,7 @@ class EOSFit:
         except RuntimeError:
             logging.exception("Fitting to EOS failed.")
             raise
-        except (RuntimeWarning, scipy.optimize.optimize.OptimizeWarning):
+        except (RuntimeWarning, scipy.optimize.OptimizeWarning):
             logging.exception("Difficulty in fitting to EOS.")
             raise
         else:

phonopy/spectrum/__init__.py

@@ -1,4 +1,5 @@
 """Routines to calculate power spectrum from MD data."""
+
 # Copyright (C) 2021 Atsushi Togo
 # All rights reserved.
 #

phonopy/spectrum/dynamic_structure_factor.py

@@ -1,4 +1,5 @@
 """Calculate dynamic structure factor at harmonic level."""
+
 # Copyright (C) 2016 Atsushi Togo
 # All rights reserved.
 #

phonopy/spectrum/velocity.py

@@ -1,4 +1,5 @@
 """Routines to analyze MD velocity data."""
+
 # Copyright (C) 2016 Atsushi Togo
 # All rights reserved.
 #

phonopy/structure/__init__.py

@@ -1,4 +1,5 @@
 """Crystal structure routines in direct and reciprocal spaces."""
+
 # Copyright (C) 2021 Atsushi Togo
 # All rights reserved.
 #

phonopy/structure/atoms.py

@@ -1,4 +1,5 @@
 """PhonopyAtoms class and routines related to atoms."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/structure/brillouin_zone.py

@@ -1,4 +1,5 @@
 """Use first Brillouin zone (Wigner–Seitz cell) to locate q-points."""
+
 # Copyright (C) 2013 Atsushi Togo
 # All rights reserved.
 #

phonopy/structure/cells.py

@@ -1,4 +1,5 @@
 """Primitive cell and supercell, and related utilities."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #
@@ -480,10 +481,30 @@ class Primitive(PhonopyAtoms):
         )
         return self.p2p_map
 
-    def get_smallest_vectors(self):
+    def get_smallest_vectors(self) -> tuple[np.ndarray, np.ndarray]:
         """Return shortest vectors and multiplicities.
 
-        See the docstring of `Primitive_get_smallest_vectors()`.
+        See also the docstring of `ShortestPairs`. The older less densen format
+        is deprecated. The detailed explaination is found in `ShortestPairs`
+        class.
+
+        Returns
+        -------
+        tuple
+            shortest_vectors : np.ndarray
+                Shortest vectors of atomic pairs in supercell from an atom in
+                the primitive cell to an atom in the supercell. The vectors are
+                given in the coordinates with respect to the primitive cell
+                basis vectors. In the dense format
+                shape=(sum(multiplicities[:,:, 0], 3), dtype='double',
+                dtype='double', order='C'.
+            multiplicities : np.ndarray
+                Number of equidistance shortest vectors. The last dimension
+                indicates [0] multipliticy at the pair of atoms in the supercell
+                and primitive cell, and [1] integral of multiplicities to this
+                pair, i.e., which indicates address used in `shortest_vectors`.
+                In the dense format, shape=(size_super, size_prim, 2),
+                dtype='int_' dtype='intc', order='C'.
 
         """
         return self._smallest_vectors, self._multiplicity
@@ -592,7 +613,9 @@ class Primitive(PhonopyAtoms):
 
         return atomic_permutations
 
-    def _get_smallest_vectors(self, supercell):
+    def _get_smallest_vectors(
+        self, supercell: PhonopyAtoms
+    ) -> tuple[np.ndarray, np.ndarray]:
         """Find shortest vectors.
 
         See the docstring of `ShortestPairs`.
@@ -1621,17 +1644,24 @@ def determinant(m):
 
 def get_primitive_matrix(
     pmat: Optional[Union[str, np.ndarray, Sequence]] = None,
-    symprec=1e-5,
-):
+    symprec: float = 1e-5,
+) -> Optional[Union[str, np.ndarray]]:
     """Find primitive matrix from primitive cell.
 
     None is equivalent to "P" but None is returned.
 
-    ``pmat`` can be
+    Parameters
+    ----------
+    pmat : str, np.ndarray, Sequency, or None
+        symbol of centring type: "P", "F", "I", "A", "C", "R"
+        "auto" : estimates a centring type.
+        3x3 matrix (can be flattened, i.e., 9 elements)
+    symprec : float
+        Tolerance.
 
-    - a symbol of centring type: "P", "F", "I", "A", "C", "R"
-    - "auto" : estimates a centring type.
-    - 3x3 matrix (can be flattened, i.e., 9 elements)
+    Returns
+    -------
+    None or 3x3 np.ndarray representing transformation matrix to primitive cell.
 
     """
     if isinstance(pmat, str) and pmat in ("P", "F", "I", "A", "C", "R", "auto"):
@@ -1660,28 +1690,43 @@ def get_primitive_matrix(
     return _pmat
 
 
-def get_primitive_matrix_by_centring(centring):
+def get_primitive_matrix_by_centring(centring) -> Optional[np.ndarray]:
     """Return primitive matrix corresponding to centring."""
     if centring == "P":
-        return [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+        return np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype="double")
     elif centring == "F":
-        return [[0, 1.0 / 2, 1.0 / 2], [1.0 / 2, 0, 1.0 / 2], [1.0 / 2, 1.0 / 2, 0]]
+        return np.array(
+            [[0.0, 1.0 / 2, 1.0 / 2], [1.0 / 2, 0, 1.0 / 2], [1.0 / 2, 1.0 / 2, 0.0]],
+            dtype="double",
+        )
     elif centring == "I":
-        return [
-            [-1.0 / 2, 1.0 / 2, 1.0 / 2],
-            [1.0 / 2, -1.0 / 2, 1.0 / 2],
-            [1.0 / 2, 1.0 / 2, -1.0 / 2],
-        ]
+        return np.array(
+            [
+                [-1.0 / 2, 1.0 / 2, 1.0 / 2],
+                [1.0 / 2, -1.0 / 2, 1.0 / 2],
+                [1.0 / 2, 1.0 / 2, -1.0 / 2],
+            ],
+            dtype="double",
+        )
     elif centring == "A":
-        return [[1, 0, 0], [0, 1.0 / 2, -1.0 / 2], [0, 1.0 / 2, 1.0 / 2]]
+        return np.array(
+            [[1.0, 0.0, 0.0], [0.0, 1.0 / 2, -1.0 / 2], [0.0, 1.0 / 2, 1.0 / 2]],
+            dtype="double",
+        )
     elif centring == "C":
-        return [[1.0 / 2, 1.0 / 2, 0], [-1.0 / 2, 1.0 / 2, 0], [0, 0, 1]]
+        return np.array(
+            [[1.0 / 2, 1.0 / 2, 0], [-1.0 / 2, 1.0 / 2, 0], [0.0, 0.0, 1.0]],
+            dtype="double",
+        )
     elif centring == "R":
-        return [
-            [2.0 / 3, -1.0 / 3, -1.0 / 3],
-            [1.0 / 3, 1.0 / 3, -2.0 / 3],
-            [1.0 / 3, 1.0 / 3, 1.0 / 3],
-        ]
+        return np.array(
+            [
+                [2.0 / 3, -1.0 / 3, -1.0 / 3],
+                [1.0 / 3, 1.0 / 3, -2.0 / 3],
+                [1.0 / 3, 1.0 / 3, 1.0 / 3],
+            ],
+            dtype="double",
+        )
     else:
         return None
 
@@ -1699,7 +1744,7 @@ def guess_primitive_matrix(unitcell: PhonopyAtoms, symprec: float = 1e-5):
     return np.array(np.dot(np.linalg.inv(tmat), pmat), dtype="double", order="C")
 
 
-def shape_supercell_matrix(smat: Optional[Union[int, float, np.ndarray]]):
+def shape_supercell_matrix(smat: Optional[Union[Sequence, np.ndarray]]) -> np.ndarray:
     """Reshape supercell matrix."""
     if smat is None:
         _smat = np.eye(3, dtype="intc", order="C")

phonopy/structure/dataset.py

@@ -1,4 +1,5 @@
 """Tools to manage displacement dataset."""
+
 # Copyright (C) 2020 Atsushi Togo
 # All rights reserved.
 #

phonopy/structure/snf.py

@@ -1,4 +1,5 @@
 """Calculation of Smith normal form of 3x3 matrix."""
+
 # Copyright (C) 2020 Atsushi Togo
 # All rights reserved.
 #

phonopy/structure/symmetry.py

@@ -1,4 +1,5 @@
 """Crystal symmetry routines."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/structure/tetrahedron_method.py

@@ -1,4 +1,5 @@
 """Core routine of linear tetrahedon method on regular grid."""
+
 # Copyright (C) 2013 Atsushi Togo
 # All rights reserved.
 #

phonopy/unfolding/__init__.py

@@ -1,4 +1,5 @@
 """Band unfolding routines."""
+
 # Copyright (C) 2019 Atsushi Togo
 # All rights reserved.
 #

phonopy/unfolding/core.py

@@ -1,4 +1,5 @@
 """Band unfolding calculation."""
+
 # Copyright (C) 2015 Atsushi Togo
 # All rights reserved.
 #

phonopy/units.py

@@ -1,4 +1,5 @@
 """Collection of physical units."""
+
 # Copyright (C) 2011 Atsushi Togo
 # All rights reserved.
 #

phonopy/utils.py

@@ -1,4 +1,5 @@
 """Utility functions."""
+
 # Copyright (C) 2022 Atsushi Togo
 # All rights reserved.
 #

phonopy/version.py

@@ -1,4 +1,5 @@
 """Phonopy version number only."""
+
 # Copyright (C) 2013 Atsushi Togo
 # All rights reserved.
 #
@@ -33,4 +34,4 @@
 # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 # POSSIBILITY OF SUCH DAMAGE.
 
-__version__ = "2.21.0"
+__version__ = "2.21.1"

setup.py

@@ -7,6 +7,7 @@ Custom parameter setting can be written in site.cfg.
 Examples are written at _get_params_from_site_cfg().
 
 """
+
 import os
 import pathlib
 import shutil

test/conftest.py

@@ -1,4 +1,5 @@
 """Pytest configuration."""
+
 from __future__ import annotations
 
 from collections.abc import Callable