From 5bf7dc938229980ea4b7b2dd856d406d63316626 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Thu, 20 Aug 2026 14:30:01 +0900 Subject: [PATCH 1/4] Refactoring of FrequencyShift --- phono3py/phonon4/api_phono4py.py | 15 ++- phono3py/phonon4/frequency_shift.py | 121 +++++++++++++++---------- phono3py/phonon4/real_to_reciprocal.py | 17 ++-- test/phonon4/test_fc4.py | 4 +- test/phonon4/test_frequency_shift.py | 5 +- 5 files changed, 94 insertions(+), 68 deletions(-) diff --git a/phono3py/phonon4/api_phono4py.py b/phono3py/phonon4/api_phono4py.py index bbcfb647..4f036a30 100644 --- a/phono3py/phonon4/api_phono4py.py +++ b/phono3py/phonon4/api_phono4py.py @@ -32,6 +32,7 @@ import numpy as np from numpy.typing import NDArray from phonopy.api_phonopy import set_data_to_phonopy_yaml +from phonopy.harmonic.dynamical_matrix import NacParams from phonopy.structure.atoms import PhonopyAtoms from phonopy.structure.cells import Primitive from phonopy.structure.symmetry import Symmetry @@ -307,6 +308,7 @@ def run_frequency_shift( grid_points: NDArray[np.int64] | list[int] | None = None, temperatures: NDArray[np.double] | list[float] | None = None, band_indices: NDArray[np.int64] | list[int] | None = None, + nac_params: NacParams | None = None, cutoff_frequency: float = 1e-4, frequency_factor_to_THz: float | None = None, ) -> NDArray[np.double]: @@ -325,12 +327,18 @@ def run_frequency_shift( Harmonic force constants defined on :attr:`supercell` (the same supercell fc4 is solved in), for the dynamical matrix. grid_points : array_like, optional - Grid-point indices into the no-symmetry ``np.ndindex(*mesh)`` - addresses. Default is ``None`` (all grid points). + BZ-grid point indices, i.e. indices into + :attr:`FrequencyShift.grid_address`. Default is ``None``, which + uses :attr:`FrequencyShift.mesh_grid_points`, one point per mesh + point. temperatures : array_like, optional Temperatures in K. Default is ``[0.0]``. band_indices : array_like, optional Band indices to compute. Default is all bands. + nac_params : NacParams, optional + Non-analytical term correction parameters. Default is None, i.e. + no correction. The correction is dropped at Gamma, where it is + direction dependent and the grid gives no direction to take. cutoff_frequency : float, optional Frequencies at or below this (THz) are skipped. Default ``1e-4``. frequency_factor_to_THz : float, optional @@ -358,13 +366,14 @@ def run_frequency_shift( mesh, temperatures=temperatures, band_indices=band_indices, + nac_params=nac_params, is_compact_fc4=is_compact_fc4, frequency_factor_to_THz=frequency_factor_to_THz, cutoff_frequency=cutoff_frequency, lang=self._lang, ) if grid_points is None: - grid_points = range(int(np.prod(mesh))) + grid_points = self._frequency_shift.mesh_grid_points return np.array([self._frequency_shift.run(int(gp)) for gp in grid_points]) def save( diff --git a/phono3py/phonon4/frequency_shift.py b/phono3py/phonon4/frequency_shift.py index 28d65f4b..441dc502 100644 --- a/phono3py/phonon4/frequency_shift.py +++ b/phono3py/phonon4/frequency_shift.py @@ -22,30 +22,18 @@ import numpy as np from numpy.typing import NDArray -from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix +from phonopy.harmonic.dynamical_matrix import NacParams, get_dynamical_matrix +from phonopy.phonon.grid import BZGrid, get_qpoints_from_bz_grid_points from phonopy.physical_units import get_physical_units from phonopy.structure.atoms import PhonopyAtoms from phonopy.structure.cells import Primitive from phono3py._lang import resolve_lang +from phono3py.phonon.func import bose_einstein +from phono3py.phonon.solver import run_phonon_solver_c, run_phonon_solver_rust from phono3py.phonon4.real_to_reciprocal import RealToReciprocalFc4 -def _bose_einstein(frequencies: NDArray[np.double], t: float) -> NDArray[np.double]: - """Return Bose-Einstein occupations for frequencies (THz) at temperature t (K). - - Frequencies at or below zero get zero occupation. - """ - units = get_physical_units() - occ = np.zeros_like(frequencies) - if t <= 0: - return occ - mask = frequencies > 0 - x = units.THzToEv * frequencies[mask] / (units.KB * t) - occ[mask] = 1.0 / (np.exp(x) - 1.0) - return occ - - class ReciprocalToNormalFc4: """Contract a reciprocal-space fc4 with eigenvectors to normal coordinates.""" @@ -141,6 +129,7 @@ def __init__( mesh: NDArray[np.int64], temperatures: NDArray[np.double] | None = None, band_indices: NDArray[np.int64] | None = None, + nac_params: NacParams | None = None, is_compact_fc4: bool = False, frequency_factor_to_THz: float | None = None, cutoff_frequency: float = 1e-4, @@ -162,6 +151,10 @@ def __init__( Temperatures in K. Default is ``[0.0]``. band_indices : array_like, optional Band indices to compute. Default is all bands. + nac_params : NacParams, optional + Non-analytical term correction parameters. Default is None, i.e. + no correction. The correction is dropped at Gamma, where it is + direction dependent and the grid gives no direction to take. is_compact_fc4 : bool, optional Whether ``fc4`` is compact. Default is False. frequency_factor_to_THz : float, optional @@ -188,13 +181,13 @@ def __init__( self._cutoff_frequency = cutoff_frequency self._lang = resolve_lang(lang) - self._dm = get_dynamical_matrix(fc2, supercell, primitive) - self._r2r = RealToReciprocalFc4( - fc4, primitive, self._mesh, is_compact_fc4, lang=self._lang + self._dm = get_dynamical_matrix( + fc2, supercell, primitive, nac_params=nac_params ) + self._r2r = RealToReciprocalFc4(fc4, primitive, is_compact_fc4, lang=self._lang) - self._grid_address = np.array( - list(np.ndindex(*self._mesh.tolist())), dtype="int64" + self._bz_grid = BZGrid( + self._mesh, lattice=primitive.cell, store_dense_gp_map=True ) self._frequencies, self._eigenvectors = self._solve_phonons() self._r2n = ReciprocalToNormalFc4( @@ -215,57 +208,85 @@ def __init__( * units.EV / (2 * np.pi * units.THz) / 8 - / np.prod(self._mesh) + / np.prod(self._bz_grid.D_diag) ) def _solve_phonons(self) -> tuple[NDArray[np.double], NDArray[np.complex128]]: - num_grid = len(self._grid_address) + num_bzgp = len(self._bz_grid.addresses) num_band = len(self._primitive) * 3 - frequencies = np.zeros((num_grid, num_band), dtype="double") - eigenvectors = np.zeros((num_grid, num_band, num_band), dtype="complex128") - for gi, address in enumerate(self._grid_address): - q = address / self._mesh - self._dm.run(q) - eigvals, eigenvectors[gi] = np.linalg.eigh(self._dm.dynamical_matrix) - frequencies[gi] = np.sqrt(np.abs(eigvals)) * np.sign(eigvals) * self._factor + frequencies = np.zeros((num_bzgp, num_band), dtype="double") + eigenvectors = np.zeros((num_bzgp, num_band, num_band), dtype="complex128") + solver = run_phonon_solver_rust if self._lang == "Rust" else run_phonon_solver_c + solver( + self._dm, + frequencies, + eigenvectors, + np.zeros(num_bzgp, dtype="byte"), # phonon_done + np.arange(num_bzgp, dtype="int64"), + self._bz_grid.addresses, + self._bz_grid.QDinv, + self._factor, + None, # No q-direction, so the NAC is dropped at Gamma. + "L", + ) return frequencies, eigenvectors + @property + def bz_grid(self) -> BZGrid: + """Return the BZ-grid the phonons and the mesh sum are defined on.""" + return self._bz_grid + @property def grid_address(self) -> NDArray[np.int64]: - """Return the (no-symmetry) grid addresses.""" - return self._grid_address + """Return the BZ-grid addresses.""" + return self._bz_grid.addresses + + @property + def mesh_grid_points(self) -> NDArray[np.int64]: + """Return one BZ-grid point per mesh point. + + The BZ-grid holds every translationally equivalent address on the BZ + surface, so it has more points than the mesh. These are the unique + representatives the weight-one mesh sum runs over. + + """ + return self._bz_grid.grg2bzg @property def frequencies(self) -> NDArray[np.double]: - """Return phonon frequencies on the grid (THz).""" + """Return phonon frequencies on the BZ-grid (THz).""" return self._frequencies def run(self, grid_point: int) -> NDArray[np.double]: - """Return frequency shifts at a grid point, shape (n_temperatures, n_bands).""" - address0 = self._grid_address[grid_point] - num_grid = len(self._grid_address) + """Return frequency shifts at a BZ-grid point. + + Returns shape (n_temperatures, n_bands). + """ + q0 = get_qpoints_from_bz_grid_points(grid_point, self._bz_grid) + gps1 = self.mesh_grid_points + qpoints1 = get_qpoints_from_bz_grid_points(gps1, self._bz_grid) num_band = len(self._primitive) * 3 # fc4_normal[gp1, band_j, band'] for the requested bands. fc4_normal = np.zeros( - (num_grid, len(self._band_indices), num_band), dtype="complex128" + (len(gps1), len(self._band_indices), num_band), dtype="complex128" ) - for gi, address1 in enumerate(self._grid_address): - quartet = np.array([-address0, address0, address1, -address1]) - fc4_reciprocal = self._r2r.run(quartet) + for i, (gp1, q1) in enumerate(zip(gps1, qpoints1, strict=True)): + fc4_reciprocal = self._r2r.run(np.array([-q0, q0, q1, -q1])) for j, band_index in enumerate(self._band_indices): - fc4_normal[gi, j] = self._r2n.run( - fc4_reciprocal, grid_point, int(band_index), gi + fc4_normal[i, j] = self._r2n.run( + fc4_reciprocal, grid_point, int(band_index), int(gp1) ) + freqs1 = self._frequencies[gps1] shifts = np.zeros((len(self._temperatures), len(self._band_indices))) for i_t, temperature in enumerate(self._temperatures): - for j in range(len(self._band_indices)): - total = 0.0 + 0.0j - for gi in range(num_grid): - occ = _bose_einstein(self._frequencies[gi], temperature) - total += ( - fc4_normal[gi, j] * self._unit_conversion * (2 * occ + 1) - ).sum() - shifts[i_t, j] = total.real + occs = np.zeros_like(freqs1) + if temperature > 0: + mask = freqs1 > 0 + occs[mask] = bose_einstein(freqs1[mask], temperature) + shifts[i_t] = ( + self._unit_conversion + * np.einsum("gjb,gb->j", fc4_normal, 2 * occs + 1).real + ) return shifts diff --git a/phono3py/phonon4/real_to_reciprocal.py b/phono3py/phonon4/real_to_reciprocal.py index bd0e2d63..c725601a 100644 --- a/phono3py/phonon4/real_to_reciprocal.py +++ b/phono3py/phonon4/real_to_reciprocal.py @@ -29,7 +29,6 @@ def __init__( self, fc4: NDArray[np.double], primitive: Primitive, - mesh: NDArray[np.int64], is_compact_fc: bool = False, lang: Literal["C", "Rust"] = "Rust", ) -> None: @@ -42,8 +41,6 @@ def __init__( ``(n_patom, N, N, N, 3, 3, 3, 3)`` layout. primitive : Primitive Primitive cell (provides p2s/s2p maps and smallest vectors). - mesh : array_like - Reciprocal sampling mesh, shape ``(3,)``. is_compact_fc : bool, optional Whether ``fc4`` is in the compact layout. Default is False. lang : {"C", "Rust"}, optional @@ -53,7 +50,6 @@ def __init__( """ self._fc4 = np.array(fc4, dtype="double", order="C") self._primitive = primitive - self._mesh = np.array(mesh, dtype="int64") self._is_compact_fc = is_compact_fc self._lang = resolve_lang(lang) self._p2s_map = np.array(primitive.p2s_map, dtype="int64") @@ -64,16 +60,15 @@ def __init__( self._num_satom = len(self._s2p_map) self._num_patom = len(self._p2s_map) - def run(self, quartet: NDArray[np.int64]) -> NDArray[np.complex128]: + def run(self, q_points: NDArray[np.double]) -> NDArray[np.complex128]: """Return fc4 in reciprocal space at the given q-point quartet. Parameters ---------- - quartet : ndarray - Four grid addresses (integers), shape ``(4, 3)``. The fractional - q-points are ``quartet / mesh``; the first q-point is the reference - (only the last three enter the phase factors), mirroring the 2015 - convention. + q_points : ndarray + Four q-points in fractional coordinates without 2pi, shape + ``(4, 3)``. The first q-point is the reference (only the last three + enter the phase factors), mirroring the 2015 convention. Returns ------- @@ -83,7 +78,7 @@ def run(self, quartet: NDArray[np.int64]) -> NDArray[np.complex128]: dtype=complex128. """ - q = np.array(quartet, dtype="double", order="C") / self._mesh + q = np.array(q_points, dtype="double", order="C") n_patom = self._num_patom fc4_reciprocal = np.zeros((n_patom,) * 4 + (3,) * 4, dtype="complex128") if self._lang == "Rust": diff --git a/test/phonon4/test_fc4.py b/test/phonon4/test_fc4.py index 78b40002..e831c902 100644 --- a/test/phonon4/test_fc4.py +++ b/test/phonon4/test_fc4.py @@ -266,8 +266,8 @@ def test_real_to_reciprocal_gamma(ph3: Phono3py) -> None: n_satom = len(ph3.supercell) rng = np.random.default_rng(0) fc4 = rng.standard_normal((n_satom, n_satom, n_satom, n_satom, 3, 3, 3, 3)) - r2r = RealToReciprocalFc4(fc4, ph3.primitive, np.array([2, 2, 2])) - rec = r2r.run(np.zeros((4, 3), dtype=int)) + r2r = RealToReciprocalFc4(fc4, ph3.primitive) + rec = r2r.run(np.zeros((4, 3), dtype="double")) n_patom = len(ph3.primitive) assert rec.shape == (n_patom, n_patom, n_patom, n_patom, 3, 3, 3, 3) diff --git a/test/phonon4/test_frequency_shift.py b/test/phonon4/test_frequency_shift.py index 5539f983..7597d9b1 100644 --- a/test/phonon4/test_frequency_shift.py +++ b/test/phonon4/test_frequency_shift.py @@ -27,8 +27,9 @@ from phonopy.structure.atoms import PhonopyAtoms from phonopy.structure.cells import Primitive, Supercell +from phono3py.phonon.func import bose_einstein from phono3py.phonon4.fc4 import set_permutation_symmetry_fc4 -from phono3py.phonon4.frequency_shift import FrequencyShift, _bose_einstein +from phono3py.phonon4.frequency_shift import FrequencyShift pytestmark = pytest.mark.filterwarnings("ignore::UserWarning") @@ -88,7 +89,7 @@ def _renormalized_shift( for nu, e in zip(freqs, eigvecs.T, strict=True): if nu <= CUTOFF: continue - n = _bose_einstein(np.array([nu]), temperature)[0] + n = bose_einstein(np.array([nu]), temperature)[0] q2 = ( units.Hbar * units.EV From 02f5028984f227d8ce213b5e4186e28ff015c4f5 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Thu, 20 Aug 2026 14:30:28 +0900 Subject: [PATCH 2/4] Fix the inter-band heat capacity matrix at degenerate frequencies --- phono3py/phonon/func.py | 4 +- phono3py/phonon/heat_capacity_matrix.py | 44 ++++++++++++++++----- phono3py/phonon4/file_IO.py | 2 +- test/conductivity/njc23/test_kappa_njc23.py | 27 +++++++------ test/conductivity/smm19/test_kappa_smm19.py | 27 +++++++------ 5 files changed, 69 insertions(+), 35 deletions(-) diff --git a/phono3py/phonon/func.py b/phono3py/phonon/func.py index a72865a5..bd9da1f0 100644 --- a/phono3py/phonon/func.py +++ b/phono3py/phonon/func.py @@ -65,8 +65,8 @@ def bose_einstein(x: NDArray[np.double], T: float) -> NDArray[np.double]: Temperature in K """ - return 1.0 / ( - np.exp(get_physical_units().THzToEv * x / (get_physical_units().KB * T)) - 1 + return 1.0 / np.expm1( + get_physical_units().THzToEv * x / (get_physical_units().KB * T) ) diff --git a/phono3py/phonon/heat_capacity_matrix.py b/phono3py/phonon/heat_capacity_matrix.py index 12c03cc6..01ae6f27 100644 --- a/phono3py/phonon/heat_capacity_matrix.py +++ b/phono3py/phonon/heat_capacity_matrix.py @@ -41,17 +41,40 @@ from phonopy.phonon.thermal_properties import mode_cv from phonopy.physical_units import get_physical_units +# Threshold on x = E / (kB T) below which two modes are treated as +# degenerate and (n_j - n_j') / (x_j - x_j') is replaced by its analytic +# limit instead of the quotient. +_DEGENERACY_TOLERANCE = 1e-6 -def _bose_einstein( - freqs: NDArray[np.double], temps: NDArray[np.double] + +def _bose_einstein_difference_ratio( + x: NDArray[np.double], n: NDArray[np.double] ) -> NDArray[np.double]: - """Bose-Einstein distribution. + r"""Return (n_j - n_j') / (x_j - x_j') with x = E / (k_B T). + + Evaluated as a quotient, this is 0/0 for degenerate modes and loses all + significant digits when two frequencies differ only by round-off, which + happens at every q-point where bands are degenerate. The exact identity + + n_{\mathbf{q}j} - n_{\mathbf{q}j'} = + -n_{\mathbf{q}j} (n_{\mathbf{q}j'} + 1) + \left( e^{x_{\mathbf{q}j} - x_{\mathbf{q}j'}} - 1 \right) - Re-implemented here because of different physical units. + removes the cancellation and gives the finite limit + -n_{\mathbf{q}j} (n_{\mathbf{q}j} + 1) at x_{\mathbf{q}j} = x_{\mathbf{q}j'}. + It is used for the nearly degenerate pairs, for which the exponential is + safely expanded to first order. """ - x = np.divide.outer(freqs, get_physical_units().KB * temps).T - return 1.0 / (np.exp(x) - 1) + d = x[:, :, None] - x[:, None, :] + n_j = n[:, :, None] + n_jp = n[:, None, :] + is_deg = np.abs(d) < _DEGENERACY_TOLERANCE + return np.where( + is_deg, + -n_j * (n_jp + 1) * (1 + d / 2), + (n_j - n_jp) / np.where(is_deg, 1, d), + ) def _mode_cv_matrix( @@ -85,10 +108,11 @@ def _mode_cv_matrix( shape=(num_temps, num_band, num_band), dtype='double', order='C'. """ - n = _bose_einstein(freqs, temps) - f_sub = np.subtract.outer(freqs, freqs) - n_sub = n[:, :, None] - n[:, None, :] - cvm = -prefactor[None, :, :] / temps[:, None, None] * n_sub / f_sub[None, :, :] + kb_temps = get_physical_units().KB * temps + x = np.divide.outer(freqs, kb_temps).T + n = 1.0 / np.expm1(x) + ratio = _bose_einstein_difference_ratio(x, n) + cvm = -prefactor[None, :, :] * ratio / (temps * kb_temps)[:, None, None] return np.ascontiguousarray(cvm) diff --git a/phono3py/phonon4/file_IO.py b/phono3py/phonon4/file_IO.py index 5dad6267..276f173e 100644 --- a/phono3py/phonon4/file_IO.py +++ b/phono3py/phonon4/file_IO.py @@ -13,9 +13,9 @@ import h5py # type: ignore[import-untyped] import numpy as np from numpy.typing import NDArray -from phono3py._version import __version__ from phonopy.file_IO import check_force_constants_indices, get_io_module_to_decompress +from phono3py._version import __version__ from phono3py.phonon4.dataset import ( count_supercells_fc4, get_displacements_and_forces_fc4, diff --git a/test/conductivity/njc23/test_kappa_njc23.py b/test/conductivity/njc23/test_kappa_njc23.py index 4e7a1295..c356f0a9 100644 --- a/test/conductivity/njc23/test_kappa_njc23.py +++ b/test/conductivity/njc23/test_kappa_njc23.py @@ -4,7 +4,10 @@ from phono3py import Phono3py -TOLERANCE = 0.25 +TOLERANCE = 0.05 +# Isotope scattering is built from the eigenvectors of degenerate bands, whose +# basis is not fixed by the eigensolver, so it varies more among architectures. +TOLERANCE_ISO = 0.3 def test_kappa_njc23_si(si_pbesol: Phono3py): @@ -25,10 +28,10 @@ def test_kappa_njc23_si_with_sigma(si_pbesol: Phono3py): ref_kappa_inter = [0.587, 0.587, 0.587, 0.0, 0.0, 0.0] si_pbesol.sigmas = [0.1] tc = _run_njc23_rta(si_pbesol, [9, 9, 9]) + si_pbesol.sigmas = None np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) - si_pbesol.sigmas = None def test_kappa_njc23_si_iso(si_pbesol: Phono3py): @@ -37,16 +40,18 @@ def test_kappa_njc23_si_iso(si_pbesol: Phono3py): ref_kappa_intra = [97.213, 97.213, 97.213, 0.0, 0.0, 0.0] ref_kappa_inter = [0.540, 0.540, 0.540, 0.0, 0.0, 0.0] tc = _run_njc23_rta(si_pbesol, [9, 9, 9], is_isotope=True) - np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) + np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE_ISO) + np.testing.assert_allclose( + ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE_ISO + ) np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) def test_kappa_njc23_nacl(nacl_pbe: Phono3py): """Test NJC23-RTA by NaCl.""" - ref_kappa = [7.93, 7.93, 7.93, 0.0, 0.0, 0.0] - ref_kappa_intra = [7.881, 7.881, 7.881, 0.0, 0.0, 0.0] - ref_kappa_inter = [0.049, 0.049, 0.049, 0.0, 0.0, 0.0] + ref_kappa = [7.956, 7.956, 7.956, 0.0, 0.0, 0.0] + ref_kappa_intra = [7.862, 7.862, 7.862, 0.0, 0.0, 0.0] + ref_kappa_inter = [0.094, 0.094, 0.094, 0.0, 0.0, 0.0] tc = _run_njc23_rta(nacl_pbe, [9, 9, 9]) np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) @@ -55,17 +60,17 @@ def test_kappa_njc23_nacl(nacl_pbe: Phono3py): def test_kappa_njc23_nacl_with_sigma(nacl_pbe: Phono3py): """Test NJC23-RTA with smearing method by NaCl.""" - ref_kappa = [7.944, 7.944, 7.944, 0.0, 0.0, 0.0] + ref_kappa = [7.988, 7.988, 7.988, 0.0, 0.0, 0.0] ref_kappa_intra = [7.895, 7.895, 7.895, 0.0, 0.0, 0.0] - ref_kappa_inter = [0.049, 0.049, 0.049, 0.0, 0.0, 0.0] + ref_kappa_inter = [0.094, 0.094, 0.094, 0.0, 0.0, 0.0] nacl_pbe.sigmas = [0.1] nacl_pbe.sigma_cutoff = 3 tc = _run_njc23_rta(nacl_pbe, [9, 9, 9]) + nacl_pbe.sigmas = None + nacl_pbe.sigma_cutoff = None np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) - nacl_pbe.sigmas = None - nacl_pbe.sigma_cutoff = None def _run_njc23_rta(ph3: Phono3py, mesh, is_isotope: bool = False): diff --git a/test/conductivity/smm19/test_kappa_smm19.py b/test/conductivity/smm19/test_kappa_smm19.py index 93381186..4da67107 100644 --- a/test/conductivity/smm19/test_kappa_smm19.py +++ b/test/conductivity/smm19/test_kappa_smm19.py @@ -4,7 +4,10 @@ from phono3py import Phono3py -TOLERANCE = 0.25 +TOLERANCE = 0.05 +# Isotope scattering is built from the eigenvectors of degenerate bands, whose +# basis is not fixed by the eigensolver, so it varies more among architectures. +TOLERANCE_ISO = 0.3 def test_kappa_smm19_si(si_pbesol: Phono3py): @@ -25,10 +28,10 @@ def test_kappa_smm19_si_with_sigma(si_pbesol: Phono3py): ref_kappa_inter = [0.592, 0.592, 0.592, 0.0, 0.0, 0.0] si_pbesol.sigmas = [0.1] tc = _run_smm19_rta(si_pbesol, [9, 9, 9]) + si_pbesol.sigmas = None np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) - si_pbesol.sigmas = None def test_kappa_smm19_si_iso(si_pbesol: Phono3py): @@ -37,16 +40,18 @@ def test_kappa_smm19_si_iso(si_pbesol: Phono3py): ref_kappa_intra = [97.213, 97.213, 97.213, 0.0, 0.0, 0.0] ref_kappa_inter = [0.545, 0.545, 0.545, 0.0, 0.0, 0.0] tc = _run_smm19_rta(si_pbesol, [9, 9, 9], is_isotope=True) - np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) + np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE_ISO) + np.testing.assert_allclose( + ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE_ISO + ) np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) def test_kappa_smm19_nacl(nacl_pbe: Phono3py): """Test SMM19-RTA by NaCl.""" - ref_kappa = [7.929, 7.929, 7.929, 0.0, 0.0, 0.0] - ref_kappa_intra = [7.881, 7.881, 7.881, 0.0, 0.0, 0.0] - ref_kappa_inter = [0.048, 0.048, 0.048, 0.0, 0.0, 0.0] + ref_kappa = [7.956, 7.956, 7.956, 0.0, 0.0, 0.0] + ref_kappa_intra = [7.862, 7.862, 7.862, 0.0, 0.0, 0.0] + ref_kappa_inter = [0.094, 0.094, 0.094, 0.0, 0.0, 0.0] tc = _run_smm19_rta(nacl_pbe, [9, 9, 9]) np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) @@ -55,17 +60,17 @@ def test_kappa_smm19_nacl(nacl_pbe: Phono3py): def test_kappa_smm19_nacl_with_sigma(nacl_pbe: Phono3py): """Test SMM19-RTA with smearing method by NaCl.""" - ref_kappa = [7.943, 7.943, 7.943, 0.0, 0.0, 0.0] + ref_kappa = [7.988, 7.988, 7.988, 0.0, 0.0, 0.0] ref_kappa_intra = [7.895, 7.895, 7.895, 0.0, 0.0, 0.0] - ref_kappa_inter = [0.049, 0.049, 0.049, 0.0, 0.0, 0.0] + ref_kappa_inter = [0.094, 0.094, 0.094, 0.0, 0.0, 0.0] nacl_pbe.sigmas = [0.1] nacl_pbe.sigma_cutoff = 3 tc = _run_smm19_rta(nacl_pbe, [9, 9, 9]) + nacl_pbe.sigmas = None + nacl_pbe.sigma_cutoff = None np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) - nacl_pbe.sigmas = None - nacl_pbe.sigma_cutoff = None def _run_smm19_rta(ph3: Phono3py, mesh, is_isotope: bool = False): From 07499aa5c80218ccd9621f66d88930357d61ca45 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 20 Aug 2026 05:31:01 +0000 Subject: [PATCH 3/4] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- phono3py/phonon4/file_IO.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phono3py/phonon4/file_IO.py b/phono3py/phonon4/file_IO.py index 276f173e..5dad6267 100644 --- a/phono3py/phonon4/file_IO.py +++ b/phono3py/phonon4/file_IO.py @@ -13,9 +13,9 @@ import h5py # type: ignore[import-untyped] import numpy as np from numpy.typing import NDArray +from phono3py._version import __version__ from phonopy.file_IO import check_force_constants_indices, get_io_module_to_decompress -from phono3py._version import __version__ from phono3py.phonon4.dataset import ( count_supercells_fc4, get_displacements_and_forces_fc4, From dc80da7542b7e8ab0bdc7b0009906866ea9d9d0c Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Thu, 20 Aug 2026 15:06:59 +0900 Subject: [PATCH 4/4] Increase test tolerance for inter-band transport --- test/conductivity/njc23/test_kappa_njc23.py | 30 ++++++++++++++++----- test/conductivity/smm19/test_kappa_smm19.py | 30 ++++++++++++++++----- 2 files changed, 46 insertions(+), 14 deletions(-) diff --git a/test/conductivity/njc23/test_kappa_njc23.py b/test/conductivity/njc23/test_kappa_njc23.py index c356f0a9..fd8ab6d7 100644 --- a/test/conductivity/njc23/test_kappa_njc23.py +++ b/test/conductivity/njc23/test_kappa_njc23.py @@ -4,10 +4,16 @@ from phono3py import Phono3py -TOLERANCE = 0.05 +# kappa and kappa_intra follow the RTA solution, which varies among +# architectures by up to ~0.06 W/m-K for Si with the tetrahedron method. 0.25 is +# the value calibrated for the conda Windows build. +TOLERANCE = 0.25 # Isotope scattering is built from the eigenvectors of degenerate bands, whose -# basis is not fixed by the eigensolver, so it varies more among architectures. +# basis is not fixed by the eigensolver, so it varies more. TOLERANCE_ISO = 0.3 +# The inter-band part is what these tests exist to pin down, and it is stable +# among architectures, so it is checked tightly. +TOLERANCE_INTER = 0.05 def test_kappa_njc23_si(si_pbesol: Phono3py): @@ -18,7 +24,9 @@ def test_kappa_njc23_si(si_pbesol: Phono3py): tc = _run_njc23_rta(si_pbesol, [9, 9, 9]) np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def test_kappa_njc23_si_with_sigma(si_pbesol: Phono3py): @@ -31,7 +39,9 @@ def test_kappa_njc23_si_with_sigma(si_pbesol: Phono3py): si_pbesol.sigmas = None np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def test_kappa_njc23_si_iso(si_pbesol: Phono3py): @@ -44,7 +54,9 @@ def test_kappa_njc23_si_iso(si_pbesol: Phono3py): np.testing.assert_allclose( ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE_ISO ) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def test_kappa_njc23_nacl(nacl_pbe: Phono3py): @@ -55,7 +67,9 @@ def test_kappa_njc23_nacl(nacl_pbe: Phono3py): tc = _run_njc23_rta(nacl_pbe, [9, 9, 9]) np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def test_kappa_njc23_nacl_with_sigma(nacl_pbe: Phono3py): @@ -70,7 +84,9 @@ def test_kappa_njc23_nacl_with_sigma(nacl_pbe: Phono3py): nacl_pbe.sigma_cutoff = None np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def _run_njc23_rta(ph3: Phono3py, mesh, is_isotope: bool = False): diff --git a/test/conductivity/smm19/test_kappa_smm19.py b/test/conductivity/smm19/test_kappa_smm19.py index 4da67107..040c2b70 100644 --- a/test/conductivity/smm19/test_kappa_smm19.py +++ b/test/conductivity/smm19/test_kappa_smm19.py @@ -4,10 +4,16 @@ from phono3py import Phono3py -TOLERANCE = 0.05 +# kappa and kappa_intra follow the RTA solution, which varies among +# architectures by up to ~0.06 W/m-K for Si with the tetrahedron method. 0.25 is +# the value calibrated for the conda Windows build. +TOLERANCE = 0.25 # Isotope scattering is built from the eigenvectors of degenerate bands, whose -# basis is not fixed by the eigensolver, so it varies more among architectures. +# basis is not fixed by the eigensolver, so it varies more. TOLERANCE_ISO = 0.3 +# The inter-band part is what these tests exist to pin down, and it is stable +# among architectures, so it is checked tightly. +TOLERANCE_INTER = 0.05 def test_kappa_smm19_si(si_pbesol: Phono3py): @@ -18,7 +24,9 @@ def test_kappa_smm19_si(si_pbesol: Phono3py): tc = _run_smm19_rta(si_pbesol, [9, 9, 9]) np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def test_kappa_smm19_si_with_sigma(si_pbesol: Phono3py): @@ -31,7 +39,9 @@ def test_kappa_smm19_si_with_sigma(si_pbesol: Phono3py): si_pbesol.sigmas = None np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def test_kappa_smm19_si_iso(si_pbesol: Phono3py): @@ -44,7 +54,9 @@ def test_kappa_smm19_si_iso(si_pbesol: Phono3py): np.testing.assert_allclose( ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE_ISO ) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def test_kappa_smm19_nacl(nacl_pbe: Phono3py): @@ -55,7 +67,9 @@ def test_kappa_smm19_nacl(nacl_pbe: Phono3py): tc = _run_smm19_rta(nacl_pbe, [9, 9, 9]) np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def test_kappa_smm19_nacl_with_sigma(nacl_pbe: Phono3py): @@ -70,7 +84,9 @@ def test_kappa_smm19_nacl_with_sigma(nacl_pbe: Phono3py): nacl_pbe.sigma_cutoff = None np.testing.assert_allclose(ref_kappa, tc.kappa.ravel(), atol=TOLERANCE) np.testing.assert_allclose(ref_kappa_intra, tc.kappa_intra.ravel(), atol=TOLERANCE) - np.testing.assert_allclose(ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE) + np.testing.assert_allclose( + ref_kappa_inter, tc.kappa_inter.ravel(), atol=TOLERANCE_INTER + ) def _run_smm19_rta(ph3: Phono3py, mesh, is_isotope: bool = False):