Set fc3-r0-average default

c/interaction.c

@@ -196,13 +196,6 @@ static void real_to_normal(double *fc3_normal_squared, const long (*g_pos)[4],
     }
     r2r_real_to_reciprocal(fc3_reciprocal, q_vecs, fc3, is_compact_fc3,
                            atom_triplets, openmp_at_bands);
-    if (atom_triplets->make_r0_average) {
-        for (i = 0; i < num_band * num_band * num_band; i++) {
-            fc3_reciprocal[i] = lapack_make_complex_double(
-                lapack_complex_double_real(fc3_reciprocal[i]) / 3,
-                lapack_complex_double_imag(fc3_reciprocal[i]) / 3);
-        }
-    }
 
 #ifdef MEASURE_R2N
     if (openmp_at_bands && num_triplets > 0) {

c/real_to_reciprocal.c

@@ -43,11 +43,18 @@
 #include "phonoc_array.h"
 #include "phonoc_const.h"
 
-static void real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
+static void real_to_reciprocal_legacy(lapack_complex_double *fc3_reciprocal,
-                               const double q_vecs[3][3], const double *fc3,
+                                      const double q_vecs[3][3],
-                               const long is_compact_fc3,
+                                      const double *fc3,
-                               const AtomTriplets *atom_triplets,
+                                      const long is_compact_fc3,
-                               const long openmp_at_bands);
+                                      const AtomTriplets *atom_triplets,
+                                      const long openmp_at_bands);
+static void real_to_reciprocal_r0_average(lapack_complex_double *fc3_reciprocal,
+                                          const double q_vecs[3][3],
+                                          const double *fc3,
+                                          const long is_compact_fc3,
+                                          const AtomTriplets *atom_triplets,
+                                          const long openmp_at_bands);
 static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
                                         const double q1[3], const double q2[3],
                                         const double *fc3,
@@ -70,21 +77,30 @@ void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
                             const long is_compact_fc3,
                             const AtomTriplets *atom_triplets,
                             const long openmp_at_bands) {
-    real_to_reciprocal(fc3_reciprocal, q_vecs, fc3, is_compact_fc3,
+    long i, num_band;
-                       atom_triplets, openmp_at_bands);
+
+    if (atom_triplets->make_r0_average) {
+        real_to_reciprocal_r0_average(fc3_reciprocal, q_vecs, fc3,
+                                      is_compact_fc3, atom_triplets,
+                                      openmp_at_bands);
+        num_band = atom_triplets->multi_dims[1] * 3;
+        for (i = 0; i < num_band * num_band * num_band; i++) {
+            fc3_reciprocal[i] = lapack_make_complex_double(
+                lapack_complex_double_real(fc3_reciprocal[i]) / 3,
+                lapack_complex_double_imag(fc3_reciprocal[i]) / 3);
+        }
+    } else {
+        real_to_reciprocal_legacy(fc3_reciprocal, q_vecs, fc3, is_compact_fc3,
+                                  atom_triplets, openmp_at_bands);
+    }
 }
 
-// Summations are performed with respect to three different lattice reference
+static void real_to_reciprocal_legacy(lapack_complex_double *fc3_reciprocal,
-// point for the index of real space fc3 when make_r0_average=True. For cubic
+                                      const double q_vecs[3][3],
-// case, these three are roughly equivalent but small difference comes from the
+                                      const double *fc3,
-// q-points in triplets used for summation implemented in
+                                      const long is_compact_fc3,
-// real_to_reciprocal_elements().
+                                      const AtomTriplets *atom_triplets,
-// --sym-fc3q makes them almost equivalent.
+                                      const long openmp_at_bands) {
-static void real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
-                               const double q_vecs[3][3], const double *fc3,
-                               const long is_compact_fc3,
-                               const AtomTriplets *atom_triplets,
-                               const long openmp_at_bands) {
     long i, j, k, l, m, n, ijk;
     long num_patom, num_band;
     lapack_complex_double pre_phase_factor, fc3_rec_elem[27], fc3_rec;
@@ -104,8 +120,55 @@ static void real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
         pre_phase_factor = get_pre_phase_factor(i, q_vecs, atom_triplets);
 
         real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[2], fc3,
-                                    is_compact_fc3, atom_triplets, i, j, k,
+                                    is_compact_fc3, atom_triplets, i, j, k, 0);
-                                    atom_triplets->make_r0_average);
+        for (l = 0; l < 3; l++) {
+            for (m = 0; m < 3; m++) {
+                for (n = 0; n < 3; n++) {
+                    fc3_rec = phonoc_complex_prod(
+                        fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factor);
+                    fc3_reciprocal[(i * 3 + l) * num_band * num_band +
+                                   (j * 3 + m) * num_band + k * 3 + n] =
+                        sum_lapack_complex_double(
+                            fc3_reciprocal[(i * 3 + l) * num_band * num_band +
+                                           (j * 3 + m) * num_band + k * 3 + n],
+                            fc3_rec);
+                }
+            }
+        }
+    }
+}
+
+// Summations are performed with respect to three different lattice reference
+// point for the index of real space fc3 when make_r0_average=True. For cubic
+// case, these three are roughly equivalent but small difference comes from the
+// q-points in triplets used for summation implemented in
+// real_to_reciprocal_elements().
+// --sym-fc3q makes them almost equivalent.
+static void real_to_reciprocal_r0_average(lapack_complex_double *fc3_reciprocal,
+                                          const double q_vecs[3][3],
+                                          const double *fc3,
+                                          const long is_compact_fc3,
+                                          const AtomTriplets *atom_triplets,
+                                          const long openmp_at_bands) {
+    long i, j, k, l, m, n, ijk;
+    long num_patom, num_band;
+    lapack_complex_double pre_phase_factor, fc3_rec_elem[27], fc3_rec;
+
+    num_patom = atom_triplets->multi_dims[1];
+    num_band = num_patom * 3;
+
+#ifdef _OPENMP
+#pragma omp parallel for private(i, j, k, l, m, n, fc3_rec_elem, fc3_rec, \
+                                     pre_phase_factor) if (openmp_at_bands)
+#endif
+    for (ijk = 0; ijk < num_patom * num_patom * num_patom; ijk++) {
+        i = ijk / (num_patom * num_patom);
+        j = (ijk - (i * num_patom * num_patom)) / num_patom;
+        k = ijk % num_patom;
+
+        pre_phase_factor = get_pre_phase_factor(i, q_vecs, atom_triplets);
+        real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[2], fc3,
+                                    is_compact_fc3, atom_triplets, i, j, k, 1);
         for (l = 0; l < 3; l++) {
             for (m = 0; m < 3; m++) {
                 for (n = 0; n < 3; n++) {
@@ -121,46 +184,40 @@ static void real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
             }
         }
 
-        if (atom_triplets->make_r0_average) {
+        // fc3_rec is stored in a way swapping jm <-> il.
-            real_to_reciprocal_elements(fc3_rec_elem, q_vecs[0], q_vecs[2], fc3,
+        pre_phase_factor = get_pre_phase_factor(j, q_vecs, atom_triplets);
-                                        is_compact_fc3, atom_triplets, i, j, k,
+        real_to_reciprocal_elements(fc3_rec_elem, q_vecs[0], q_vecs[2], fc3,
-                                        2);
+                                    is_compact_fc3, atom_triplets, j, i, k, 2);
-            for (l = 0; l < 3; l++) {
+        for (l = 0; l < 3; l++) {
-                for (m = 0; m < 3; m++) {
+            for (m = 0; m < 3; m++) {
-                    for (n = 0; n < 3; n++) {
+                for (n = 0; n < 3; n++) {
-                        // fc3_rec is stored in a way swapping jm <-> il.
+                    fc3_rec = phonoc_complex_prod(
-                        fc3_rec = phonoc_complex_prod(
+                        fc3_rec_elem[m * 9 + l * 3 + n], pre_phase_factor);
-                            fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factor);
+                    fc3_reciprocal[(i * 3 + l) * num_band * num_band +
-                        fc3_reciprocal[(j * 3 + m) * num_band * num_band +
+                                   (j * 3 + m) * num_band + k * 3 + n] =
-                                       (i * 3 + l) * num_band + k * 3 + n] =
+                        sum_lapack_complex_double(
-                            sum_lapack_complex_double(
+                            fc3_reciprocal[(i * 3 + l) * num_band * num_band +
-                                fc3_reciprocal[(j * 3 + m) * num_band *
+                                           (j * 3 + m) * num_band + k * 3 + n],
-                                                   num_band +
+                            fc3_rec);
-                                               (i * 3 + l) * num_band + k * 3 +
-                                               n],
-                                fc3_rec);
-                    }
                 }
             }
+        }
 
-            real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[0], fc3,
+        // fc3_rec is stored in a way swapping kn <-> il.
-                                        is_compact_fc3, atom_triplets, i, j, k,
+        pre_phase_factor = get_pre_phase_factor(k, q_vecs, atom_triplets);
-                                        3);
+        real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[0], fc3,
-            for (l = 0; l < 3; l++) {
+                                    is_compact_fc3, atom_triplets, k, j, i, 3);
-                for (m = 0; m < 3; m++) {
+        for (l = 0; l < 3; l++) {
-                    for (n = 0; n < 3; n++) {
+            for (m = 0; m < 3; m++) {
-                        // fc3_rec is stored in a way swapping kn <-> il.
+                for (n = 0; n < 3; n++) {
-                        fc3_rec = phonoc_complex_prod(
+                    fc3_rec = phonoc_complex_prod(
-                            fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factor);
+                        fc3_rec_elem[n * 9 + m * 3 + l], pre_phase_factor);
-                        fc3_reciprocal[(k * 3 + n) * num_band * num_band +
+                    fc3_reciprocal[(i * 3 + l) * num_band * num_band +
-                                       (j * 3 + m) * num_band + i * 3 + l] =
+                                   (j * 3 + m) * num_band + k * 3 + n] =
-                            sum_lapack_complex_double(
+                        sum_lapack_complex_double(
-                                fc3_reciprocal[(k * 3 + n) * num_band *
+                            fc3_reciprocal[(i * 3 + l) * num_band * num_band +
-                                                   num_band +
+                                           (j * 3 + m) * num_band + k * 3 + n],
-                                               (j * 3 + m) * num_band + i * 3 +
+                            fc3_rec);
-                                               l],
-                                fc3_rec);
-                    }
                 }
             }
         }

phono3py/api_phono3py.py

@@ -231,12 +231,6 @@ class Phono3py:
         self._use_grg = use_grg
         self._SNF_coordinates = SNF_coordinates
 
-        if not make_r0_average:
-            warnings.warn(
-                "Phono3py init parameter of make_r0_average is deprecated. "
-                "This is always True but exists for backward compatibility.",
-                DeprecationWarning,
-            )
         self._make_r0_average = make_r0_average
 
         self._cutoff_frequency = cutoff_frequency

phono3py/cui/load.py

@@ -42,6 +42,7 @@ import numpy as np
 import phonopy.cui.load_helper as load_helper
 from phonopy.harmonic.force_constants import show_drift_force_constants
 from phonopy.interface.calculator import get_default_physical_units
+from phonopy.structure.atoms import PhonopyAtoms
 from phonopy.structure.cells import determinant
 
 from phono3py import Phono3py
@@ -52,33 +53,36 @@ from phono3py.phonon3.fc3 import show_drift_fc3
 
 
 def load(
-    phono3py_yaml=None,  # phono3py.yaml-like must be the first argument.
+    phono3py_yaml: Optional[
-    supercell_matrix=None,
+        Union[str, bytes, os.PathLike]
-    primitive_matrix=None,
+    ] = None,  # phono3py.yaml-like must be the first argument.
-    phonon_supercell_matrix=None,
+    supercell_matrix: Optional[Union[Sequence, np.ndarray]] = None,
-    is_nac=True,
+    primitive_matrix: Optional[Union[Sequence, np.ndarray]] = None,
-    calculator=None,
+    phonon_supercell_matrix: Optional[Union[Sequence, np.ndarray]] = None,
-    unitcell=None,
+    is_nac: bool = True,
-    supercell=None,
+    calculator: Optional[str] = None,
-    nac_params=None,
+    unitcell: Optional[PhonopyAtoms] = None,
-    unitcell_filename=None,
+    supercell: Optional[PhonopyAtoms] = None,
-    supercell_filename=None,
+    nac_params: Optional[dict] = None,
-    born_filename=None,
+    unitcell_filename: Optional[Union[str, bytes, os.PathLike]] = None,
-    forces_fc3_filename: Optional[Union[os.PathLike, Sequence]] = None,
+    supercell_filename: Optional[Union[str, bytes, os.PathLike]] = None,
-    forces_fc2_filename: Optional[Union[os.PathLike, Sequence]] = None,
+    born_filename: Optional[Union[str, bytes, os.PathLike]] = None,
-    fc3_filename=None,
+    forces_fc3_filename: Optional[Union[str, bytes, os.PathLike]] = None,
-    fc2_filename=None,
+    forces_fc2_filename: Optional[Union[str, bytes, os.PathLike]] = None,
-    fc_calculator=None,
+    fc3_filename: Optional[Union[str, bytes, os.PathLike]] = None,
-    fc_calculator_options=None,
+    fc2_filename: Optional[Union[str, bytes, os.PathLike]] = None,
-    factor=None,
+    fc_calculator: Optional[str] = None,
-    produce_fc=True,
+    fc_calculator_options: Optional[str] = None,
-    is_symmetry=True,
+    factor: Optional[float] = None,
-    symmetrize_fc=True,
+    produce_fc: bool = True,
-    is_mesh_symmetry=True,
+    is_symmetry: bool = True,
-    is_compact_fc=False,
+    symmetrize_fc: bool = True,
-    use_grg=False,
+    is_mesh_symmetry: bool = True,
-    symprec=1e-5,
+    is_compact_fc: bool = False,
-    log_level=0,
+    use_grg: bool = False,
+    make_r0_average: bool = True,
+    symprec: float = 1e-5,
+    log_level: int = 0,
 ) -> Phono3py:
     """Create Phono3py instance from parameters and/or input files.
 
@@ -227,6 +231,11 @@ def load(
         cells. Default is False.
     use_grg : bool, optional
         Use generalized regular grid when True. Default is False.
+    make_r0_average : bool, optional
+        fc3 transformation from real to reciprocal space is done
+        around three atoms and averaged when True. Default is False, i.e.,
+        only around the first atom. Setting False is for rough compatibility
+        with v2.x. Default is True.
     symprec : float, optional
         Tolerance used to find crystal symmetry. Default is 1e-5.
     log_level : int, optional
@@ -297,6 +306,7 @@ def load(
         is_symmetry=is_symmetry,
         is_mesh_symmetry=is_mesh_symmetry,
         use_grg=use_grg,
+        make_r0_average=make_r0_average,
         calculator=calculator,
         log_level=log_level,
     )

test/conductivity/test_kappa_LBTE.py

@@ -3,15 +3,13 @@ import numpy as np
 
 from phono3py.api_phono3py import Phono3py
 
-si_pbesol_kappa_LBTE = [111.117, 111.117, 111.117, 0, 0, 0]
-si_pbesol_kappa_LBTE_redcol = [63.019, 63.019, 63.019, 0, 0, 0]
-aln_lda_kappa_LBTE = [2.313066e02, 2.313066e02, 2.483627e02, 0, 0, 0]
-aln_lda_kappa_LBTE_with_sigma = [2.500303e02, 2.500303e02, 2.694047e02, 0, 0, 0]
-aln_lda_kappa_LBTE_with_r0_ave = [2.342499e02, 2.342499e02, 2.540009e02, 0, 0, 0]
-
 
 def test_kappa_LBTE(si_pbesol: Phono3py):
     """Test for symmetry reduced collision matrix."""
+    if si_pbesol._make_r0_average:
+        ref_kappa = [110.896, 110.896, 110.896, 0, 0, 0]
+    else:
+        ref_kappa = [111.149, 111.149, 111.149, 0, 0, 0]
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
     si_pbesol.run_thermal_conductivity(
@@ -21,11 +19,16 @@ def test_kappa_LBTE(si_pbesol: Phono3py):
         ],
     )
     kappa = si_pbesol.thermal_conductivity.kappa.ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_LBTE, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa, kappa, atol=0.3)
 
 
 def test_kappa_LBTE_full_colmat(si_pbesol: Phono3py):
     """Test for full collision matrix."""
+    if si_pbesol._make_r0_average:
+        ref_kappa = [62.497, 62.497, 62.497, 0, 0, 0]
+    else:
+        ref_kappa = [62.777, 62.777, 62.777, 0, 0, 0]
+
     si_pbesol.mesh_numbers = [5, 5, 5]
     si_pbesol.init_phph_interaction()
     si_pbesol.run_thermal_conductivity(
@@ -36,11 +39,16 @@ def test_kappa_LBTE_full_colmat(si_pbesol: Phono3py):
         is_reducible_collision_matrix=True,
     )
     kappa = si_pbesol.thermal_conductivity.kappa.ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_LBTE_redcol, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa, kappa, atol=0.3)
 
 
 def test_kappa_LBTE_aln(aln_lda: Phono3py):
     """Test direct solution by AlN."""
+    if aln_lda._make_r0_average:
+        ref_kappa = [234.141, 234.141, 254.006, 0, 0, 0]
+    else:
+        ref_kappa = [231.191, 231.191, 248.367, 0, 0, 0]
+
     aln_lda.mesh_numbers = [7, 7, 5]
     aln_lda.init_phph_interaction()
     aln_lda.run_thermal_conductivity(
@@ -51,29 +59,15 @@ def test_kappa_LBTE_aln(aln_lda: Phono3py):
     )
     kappa = aln_lda.thermal_conductivity.kappa.ravel()
     # print(", ".join([f"{k:e}" for k in kappa]))
-    np.testing.assert_allclose(aln_lda_kappa_LBTE, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa, kappa, atol=0.3)
-
-
-def test_kappa_LBTE_aln_with_r0_ave(aln_lda: Phono3py):
-    """Test direct solution by AlN."""
-    aln_lda.mesh_numbers = [7, 7, 5]
-    make_r0_average_orig = aln_lda._make_r0_average
-    aln_lda._make_r0_average = True
-    aln_lda.init_phph_interaction()
-    aln_lda._make_r0_average = make_r0_average_orig
-    aln_lda.run_thermal_conductivity(
-        is_LBTE=True,
-        temperatures=[
-            300,
-        ],
-    )
-    kappa = aln_lda.thermal_conductivity.kappa.ravel()
-    # print(", ".join([f"{k:e}" for k in kappa]))
-    np.testing.assert_allclose(aln_lda_kappa_LBTE_with_r0_ave, kappa, atol=0.5)
 
 
 def test_kappa_LBTE_aln_with_sigma(aln_lda: Phono3py):
     """Test direct solution by AlN."""
+    if aln_lda._make_r0_average:
+        ref_kappa = [254.111, 254.111, 271.406, 0, 0, 0]
+    else:
+        ref_kappa = [250.030, 250.030, 269.405, 0, 0, 0]
     aln_lda.sigmas = [
         0.1,
     ]
@@ -91,4 +85,4 @@ def test_kappa_LBTE_aln_with_sigma(aln_lda: Phono3py):
     aln_lda.sigmas = None
     aln_lda.sigma_cutoff = None
     # print(", ".join([f"{k:e}" for k in kappa]))
-    np.testing.assert_allclose(aln_lda_kappa_LBTE_with_sigma, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa, kappa, atol=0.3)

test/conductivity/test_kappa_LBTE_Wigner.py

@@ -1,19 +1,19 @@
 """Tests for direct solution of LBTE."""
 import numpy as np
+import pytest
 
 from phono3py.api_phono3py import Phono3py
 
-si_pbesol_kappa_P_LBTE = [111.123, 111.123, 111.123, 0, 0, 0]  # old value 111.117
-si_pbesol_kappa_C = [0.167, 0.167, 0.167, 0.000, 0.000, 0.000]
-
-si_pbesol_kappa_P_LBTE_redcol = [62.783, 62.783, 62.783, 0, 0, 0]  # old value 63.019
-si_pbesol_kappa_C_redcol = (
-    -1
-)  # coherences conductivity is not implemented for is_reducible_collision_matrix=True,
-
 
 def test_kappa_LBTE(si_pbesol: Phono3py):
     """Test for symmetry reduced collision matrix."""
+    if si_pbesol._make_r0_average:
+        ref_kappa_P_LBTE = [110.896, 110.896, 110.896, 0, 0, 0]
+        ref_kappa_C = [0.166, 0.166, 0.166, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_P_LBTE = [111.149, 111.149, 111.149, 0, 0, 0]
+        ref_kappa_C = [0.166, 0.166, 0.166, 0.000, 0.000, 0.000]
+
     si_pbesol.mesh_numbers = [9, 9, 9]
     si_pbesol.init_phph_interaction()
     si_pbesol.run_thermal_conductivity(
@@ -26,14 +26,20 @@ def test_kappa_LBTE(si_pbesol: Phono3py):
     # kappa = si_pbesol.thermal_conductivity.kappa.ravel()
     kappa_P = si_pbesol.thermal_conductivity.kappa_P_exact.ravel()
     kappa_C = si_pbesol.thermal_conductivity.kappa_C.ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_P_LBTE, kappa_P, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_P_LBTE, kappa_P, atol=0.5)
-    np.testing.assert_allclose(si_pbesol_kappa_C, kappa_C, atol=0.02)
+    np.testing.assert_allclose(ref_kappa_C, kappa_C, atol=0.02)
 
 
-'''
+@pytest.mark.skip(
-#coherences conductivity is not implemented for is_reducible_collision_matrix=True,
+    reason=(
+        "coherences conductivity is not implemented for "
+        "is_reducible_collision_matrix=True"
+    )
+)
 def test_kappa_LBTE_full_colmat(si_pbesol: Phono3py):
     """Test for full collision matrix."""
+    si_pbesol_kappa_P_LBTE_redcol = [62.783, 62.783, 62.783, 0, 0, 0]
+    si_pbesol_kappa_C_redcol = -1
     si_pbesol.mesh_numbers = [5, 5, 5]
     si_pbesol.init_phph_interaction()
     si_pbesol.run_thermal_conductivity(
@@ -47,4 +53,3 @@ def test_kappa_LBTE_full_colmat(si_pbesol: Phono3py):
     kappa_C = si_pbesol.thermal_conductivity.kappa_C.ravel()
     np.testing.assert_allclose(si_pbesol_kappa_P_LBTE_redcol, kappa_P, atol=0.5)
     np.testing.assert_allclose(si_pbesol_kappa_C_redcol, kappa_C, atol=0.02)
-'''

test/conductivity/test_kappa_RTA.py

@@ -1,131 +1,132 @@
 """Test for Conductivity_RTA.py."""
+import itertools
+
 import numpy as np
 import pytest
 
 from phono3py import Phono3py
 
-si_pbesol_kappa_RTA = [107.991, 107.991, 107.991, 0, 0, 0]
-si_pbesol_kappa_RTA_with_sigmas = [109.6985, 109.6985, 109.6985, 0, 0, 0]
-si_pbesol_kappa_RTA_iso = [96.92419, 96.92419, 96.92419, 0, 0, 0]
-si_pbesol_kappa_RTA_with_sigmas_iso = [96.03248, 96.03248, 96.03248, 0, 0, 0]
-si_pbesol_kappa_RTA_si_nosym = [
-    38.242347,
-    38.700219,
-    39.198018,
-    0.3216,
-    0.207731,
-    0.283,
-]
-si_pbesol_kappa_RTA_si_nomeshsym = [81.31304, 81.31304, 81.31304, 0, 0, 0]
-si_pbesol_kappa_RTA_grg = [93.99526, 93.99526, 93.99526, 0, 0, 0]
-si_pbesol_kappa_RTA_grg_iso = [104.281556, 104.281556, 104.281556, 0, 0, 0]
-si_pbesol_kappa_RTA_grg_sigma_iso = [107.2834, 107.2834, 107.2834, 0, 0, 0]
-nacl_pbe_kappa_RTA = [7.72798252, 7.72798252, 7.72798252, 0, 0, 0]
-nacl_pbe_kappa_RTA_with_sigma = [7.71913708, 7.71913708, 7.71913708, 0, 0, 0]
 
-
+@pytest.mark.parametrize(
-aln_lda_kappa_RTA = [203.304059, 203.304059, 213.003125, 0, 0, 0]
+    "openmp_per_triplets,is_full_pp,is_compact_fc",
-aln_lda_kappa_RTA_r0_ave = [2.06489355e02, 2.06489355e02, 2.19821864e02, 0, 0, 0]
+    itertools.product([False, True], [False, True], [False, True]),
-aln_lda_kappa_RTA_with_sigmas = [213.820000, 213.820000, 224.800121, 0, 0, 0]
+)
-aln_lda_kappa_RTA_with_sigmas_r0_ave = [
-    2.17597885e02,
-    2.17597885e02,
-    2.30098660e02,
-    0,
-    0,
-    0,
-]
-
-
-@pytest.mark.parametrize("openmp_per_triplets", [True, False])
 def test_kappa_RTA_si(
     si_pbesol: Phono3py,
+    si_pbesol_compact_fc: Phono3py,
     openmp_per_triplets: bool,
+    is_full_pp: bool,
+    is_compact_fc: bool,
 ):
     """Test RTA by Si."""
+    if is_compact_fc:
+        ph3 = si_pbesol_compact_fc
+    else:
+        ph3 = si_pbesol
+    if ph3._make_r0_average:
+        if is_compact_fc:
+            ref_kappa_RTA = [107.794, 107.794, 107.794, 0, 0, 0]
+        else:
+            ref_kappa_RTA = [107.694, 107.694, 107.694, 0, 0, 0]
+    else:
+        if is_compact_fc:
+            ref_kappa_RTA = [107.956, 107.956, 107.956, 0, 0, 0]
+        else:
+            ref_kappa_RTA = [107.844, 107.844, 107.844, 0, 0, 0]
+
     kappa = _get_kappa(
-        si_pbesol,
+        ph3,
         [9, 9, 9],
+        is_full_pp=is_full_pp,
         openmp_per_triplets=openmp_per_triplets,
     ).ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_RTA, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA, kappa, atol=0.5)
-
-
-@pytest.mark.parametrize("openmp_per_triplets", [True, False])
-def test_kappa_RTA_si_full_pp(si_pbesol: Phono3py, openmp_per_triplets: bool):
-    """Test RTA with full-pp by Si."""
-    kappa = _get_kappa(
-        si_pbesol, [9, 9, 9], is_full_pp=True, openmp_per_triplets=openmp_per_triplets
-    ).ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_RTA, kappa, atol=0.5)
 
 
 def test_kappa_RTA_si_iso(si_pbesol: Phono3py):
     """Test RTA with isotope scattering by Si."""
+    if si_pbesol._make_r0_average:
+        ref_kappa_RTA_iso = [97.296, 97.296, 97.296, 0, 0, 0]
+    else:
+        ref_kappa_RTA_iso = [97.346, 97.346, 97.346, 0, 0, 0]
+
     kappa = _get_kappa(si_pbesol, [9, 9, 9], is_isotope=True).ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_RTA_iso, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_iso, kappa, atol=0.5)
 
 
-@pytest.mark.parametrize("openmp_per_triplets", [True, False])
+@pytest.mark.parametrize(
-def test_kappa_RTA_si_with_sigma(si_pbesol: Phono3py, openmp_per_triplets: bool):
+    "openmp_per_triplets,is_full_pp", itertools.product([False, True], [False, True])
+)
+def test_kappa_RTA_si_with_sigma(
+    si_pbesol: Phono3py, openmp_per_triplets: bool, is_full_pp: bool
+):
     """Test RTA with smearing method by Si."""
+    if si_pbesol._make_r0_average:
+        ref_kappa_RTA_with_sigmas = [109.999, 109.999, 109.999, 0, 0, 0]
+    else:
+        ref_kappa_RTA_with_sigmas = [109.699, 109.699, 109.699, 0, 0, 0]
+
     si_pbesol.sigmas = [
         0.1,
     ]
     kappa = _get_kappa(
-        si_pbesol, [9, 9, 9], openmp_per_triplets=openmp_per_triplets
+        si_pbesol,
+        [9, 9, 9],
+        is_full_pp=is_full_pp,
+        openmp_per_triplets=openmp_per_triplets,
     ).ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_RTA_with_sigmas, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_with_sigmas, kappa, atol=0.5)
-    si_pbesol.sigmas = None
-
-
-def test_kappa_RTA_si_with_sigma_full_pp(si_pbesol: Phono3py):
-    """Test RTA with smearing method and full-pp by Si."""
-    si_pbesol.sigmas = [
-        0.1,
-    ]
-    kappa = _get_kappa(si_pbesol, [9, 9, 9], is_full_pp=True).ravel()
-    print(kappa)
-    np.testing.assert_allclose(si_pbesol_kappa_RTA_with_sigmas, kappa, atol=0.5)
     si_pbesol.sigmas = None
 
 
 def test_kappa_RTA_si_with_sigma_iso(si_pbesol: Phono3py):
     """Test RTA with smearing method and isotope scattering by Si."""
+    if si_pbesol._make_r0_average:
+        ref_kappa_RTA_with_sigmas_iso = [96.368, 96.368, 96.368, 0, 0, 0]
+    else:
+        ref_kappa_RTA_with_sigmas_iso = [96.032, 96.032, 96.032, 0, 0, 0]
+
     si_pbesol.sigmas = [
         0.1,
     ]
     kappa = _get_kappa(si_pbesol, [9, 9, 9], is_isotope=True).ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_RTA_with_sigmas_iso, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_with_sigmas_iso, kappa, atol=0.5)
     si_pbesol.sigmas = None
 
 
-def test_kappa_RTA_si_compact_fc(si_pbesol_compact_fc: Phono3py):
-    """Test RTA with compact-fc by Si."""
-    kappa = _get_kappa(si_pbesol_compact_fc, [9, 9, 9]).ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_RTA, kappa, atol=0.5)
-
-
 def test_kappa_RTA_si_nosym(si_pbesol: Phono3py, si_pbesol_nosym: Phono3py):
     """Test RTA without considering symmetry by Si."""
+    if si_pbesol_nosym._make_r0_average:
+        ref_kappa_RTA_si_nosym = [38.315, 38.616, 39.093, 0.221, 0.166, 0.284]
+    else:
+        ref_kappa_RTA_si_nosym = [38.342, 38.650, 39.105, 0.224, 0.170, 0.288]
+
     si_pbesol_nosym.fc2 = si_pbesol.fc2
     si_pbesol_nosym.fc3 = si_pbesol.fc3
     kappa = _get_kappa(si_pbesol_nosym, [4, 4, 4]).reshape(-1, 3).sum(axis=1)
-    kappa_ref = np.reshape(si_pbesol_kappa_RTA_si_nosym, (-1, 3)).sum(axis=1)
+    kappa_ref = np.reshape(ref_kappa_RTA_si_nosym, (-1, 3)).sum(axis=1)
     np.testing.assert_allclose(kappa_ref / 3, kappa / 3, atol=0.5)
 
 
 def test_kappa_RTA_si_nomeshsym(si_pbesol: Phono3py, si_pbesol_nomeshsym: Phono3py):
     """Test RTA without considering mesh symmetry by Si."""
+
+    if si_pbesol_nomeshsym._make_r0_average:
+        ref_kappa_RTA_si_nomeshsym = [81.147, 81.147, 81.147, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_RTA_si_nomeshsym = [81.263, 81.263, 81.263, 0.000, 0.000, 0.000]
     si_pbesol_nomeshsym.fc2 = si_pbesol.fc2
     si_pbesol_nomeshsym.fc3 = si_pbesol.fc3
     kappa = _get_kappa(si_pbesol_nomeshsym, [7, 7, 7]).ravel()
-    kappa_ref = si_pbesol_kappa_RTA_si_nomeshsym
+    np.testing.assert_allclose(ref_kappa_RTA_si_nomeshsym, kappa, atol=0.5)
-    np.testing.assert_allclose(kappa_ref, kappa, atol=0.5)
 
 
 def test_kappa_RTA_si_grg(si_pbesol_grg: Phono3py):
     """Test RTA by Si with GR-grid."""
+    if si_pbesol_grg._make_r0_average:
+        ref_kappa_RTA_grg = [94.293, 94.293, 94.293, 0, 0, 0]
+    else:
+        ref_kappa_RTA_grg = [94.306, 94.306, 94.306, 0, 0, 0]
     mesh = 20
     ph3 = si_pbesol_grg
     ph3.mesh_numbers = mesh
@@ -150,11 +151,16 @@ def test_kappa_RTA_si_grg(si_pbesol_grg: Phono3py):
     Q = ph3.grid.Q
     np.testing.assert_equal(np.dot(P, np.dot(A, Q)), np.diag(D_diag))
 
-    np.testing.assert_allclose(si_pbesol_kappa_RTA_grg, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_grg, kappa, atol=0.5)
 
 
 def test_kappa_RTA_si_grg_iso(si_pbesol_grg: Phono3py):
     """Test RTA with isotope scattering by Si with GR-grid.."""
+    if si_pbesol_grg._make_r0_average:
+        ref_kappa_RTA_grg_iso = [104.290, 104.290, 104.290, 0, 0, 0]
+    else:
+        ref_kappa_RTA_grg_iso = [104.425, 104.425, 104.425, 0, 0, 0]
+
     mesh = 30
     ph3 = si_pbesol_grg
     ph3.mesh_numbers = mesh
@@ -166,12 +172,16 @@ def test_kappa_RTA_si_grg_iso(si_pbesol_grg: Phono3py):
         is_isotope=True,
     )
     kappa = ph3.thermal_conductivity.kappa.ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_RTA_grg_iso, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_grg_iso, kappa, atol=0.5)
     np.testing.assert_equal(ph3.grid.grid_matrix, [[-6, 6, 6], [6, -6, 6], [6, 6, -6]])
 
 
 def test_kappa_RTA_si_grg_sigma_iso(si_pbesol_grg: Phono3py):
     """Test RTA with isotope scattering by Si with GR-grid.."""
+    if si_pbesol_grg._make_r0_average:
+        ref_kappa_RTA_grg_sigma_iso = [107.264, 107.264, 107.264, 0, 0, 0]
+    else:
+        ref_kappa_RTA_grg_sigma_iso = [107.283, 107.283, 107.283, 0, 0, 0]
     mesh = 30
     ph3 = si_pbesol_grg
     ph3.sigmas = [
@@ -186,7 +196,7 @@ def test_kappa_RTA_si_grg_sigma_iso(si_pbesol_grg: Phono3py):
         is_isotope=True,
     )
     kappa = ph3.thermal_conductivity.kappa.ravel()
-    np.testing.assert_allclose(si_pbesol_kappa_RTA_grg_sigma_iso, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_grg_sigma_iso, kappa, atol=0.5)
     ph3.sigmas = None
 
 
@@ -204,166 +214,112 @@ def test_kappa_RTA_si_N_U(si_pbesol):
         is_N_U=is_N_U,
     )
     gN, gU = ph3.thermal_conductivity.get_gamma_N_U()
+    # for line in gN.reshape(-1, 6):
+    #     print("[", ",".join([f"{val:.8f}" for val in line]), "],")
+    # for line in gU.reshape(-1, 6):
+    #     print("[", ",".join([f"{val:.8f}" for val in line]), "],")
 
-    gN_ref = [
+    if si_pbesol._make_r0_average:
-        0.00000000,
+        gN_ref = [
-        0.00000000,
+            [0.00000000, 0.00000000, 0.00000000, 0.07898606, 0.07898606, 0.07898606],
-        0.00000000,
+            [0.00079647, 0.00079647, 0.00913611, 0.01911102, 0.04553001, 0.04553001],
-        0.07402084,
+            [0.00173868, 0.00173868, 0.01404937, 0.00201732, 0.03354033, 0.03354033],
-        0.07402084,
+            [0.00223616, 0.00223616, 0.01039331, 0.02860916, 0.02860916, 0.01987485],
-        0.07402084,
+            [0.00291788, 0.00356241, 0.02858543, 0.00367742, 0.02065990, 0.01533763],
-        0.00078535,
+            [0.00146333, 0.00343175, 0.01596851, 0.00626596, 0.02431620, 0.01091592],
-        0.00078535,
+            [0.00396766, 0.00396766, 0.00159161, 0.00159161, 0.01479018, 0.01479018],
-        0.00917995,
+            [0.00682740, 0.00682740, 0.03983399, 0.03983399, 0.02728522, 0.02728522],
-        0.02178049,
+        ]
-        0.04470075,
+        gU_ref = [
-        0.04470075,
+            [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
-        0.00173337,
+            [0.00015184, 0.00015184, 0.00075965, 0.00736940, 0.00114177, 0.00114177],
-        0.00173337,
+            [0.00022400, 0.00022400, 0.00072237, 0.00000112, 0.00022016, 0.00022016],
-        0.01240191,
+            [0.00079188, 0.00079188, 0.00106579, 0.00418717, 0.00418717, 0.00712761],
-        0.00198981,
+            [0.00219252, 0.00262840, 0.01927670, 0.00491388, 0.01254730, 0.00519414],
-        0.03165195,
+            [0.00146999, 0.00168024, 0.01596274, 0.00641979, 0.00597353, 0.00859841],
-        0.03165195,
+            [0.00307881, 0.00307881, 0.00036554, 0.00036554, 0.01176737, 0.01176737],
-        0.00224713,
+            [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
-        0.00224713,
+        ]
-        0.00860026,
+    else:
-        0.03083611,
+        gN_ref = [
-        0.03083611,
+            [0.00000000, 0.00000000, 0.00000000, 0.07832198, 0.07832198, 0.07832198],
-        0.02142118,
+            [0.00079578, 0.00079578, 0.00909025, 0.01917012, 0.04557656, 0.04557656],
-        0.00277534,
+            [0.00176235, 0.00176235, 0.01414436, 0.00204092, 0.03361112, 0.03361112],
-        0.00330170,
+            [0.00221919, 0.00221919, 0.01020133, 0.02889554, 0.02889554, 0.01995543],
-        0.02727451,
+            [0.00292189, 0.00356099, 0.02855954, 0.00370530, 0.02071850, 0.01533334],
-        0.00356415,
+            [0.00147656, 0.00342335, 0.01589430, 0.00630792, 0.02427768, 0.01099287],
-        0.01847744,
+            [0.00400675, 0.00400675, 0.00162186, 0.00162186, 0.01478489, 0.01478489],
-        0.01320643,
+            [0.00676576, 0.00676576, 0.03984290, 0.03984290, 0.02715102, 0.02715102],
-        0.00155072,
+        ]
-        0.00365611,
+        gU_ref = [
-        0.01641919,
+            [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
-        0.00650083,
+            [0.00015178, 0.00015178, 0.00076936, 0.00727539, 0.00113112, 0.00113112],
-        0.02576069,
+            [0.00022696, 0.00022696, 0.00072558, 0.00000108, 0.00021968, 0.00021968],
-        0.01161589,
+            [0.00079397, 0.00079397, 0.00111068, 0.00424761, 0.00424761, 0.00697760],
-        0.00411969,
+            [0.00219456, 0.00261878, 0.01928629, 0.00490046, 0.01249235, 0.00517685],
-        0.00411969,
+            [0.00149539, 0.00161230, 0.01594274, 0.00653088, 0.00593572, 0.00849890],
-        0.00168211,
+            [0.00311169, 0.00311169, 0.00036610, 0.00036610, 0.01171667, 0.01171667],
-        0.00168211,
+            [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
-        0.01560092,
+        ]
-        0.01560092,
-        0.00620091,
-        0.00620091,
-        0.03764912,
-        0.03764912,
-        0.02668523,
-        0.02668523,
-    ]
-    gU_ref = [
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00015178,
-        0.00015178,
-        0.00076936,
-        0.00727539,
-        0.00113112,
-        0.00113112,
-        0.00022696,
-        0.00022696,
-        0.00072558,
-        0.00000108,
-        0.00021968,
-        0.00021968,
-        0.00079397,
-        0.00079397,
-        0.00111068,
-        0.00424761,
-        0.00424761,
-        0.00697760,
-        0.00221593,
-        0.00259510,
-        0.01996296,
-        0.00498962,
-        0.01258375,
-        0.00513825,
-        0.00148802,
-        0.00161955,
-        0.01589219,
-        0.00646134,
-        0.00577275,
-        0.00849711,
-        0.00313208,
-        0.00313208,
-        0.00036610,
-        0.00036610,
-        0.01135335,
-        0.01135335,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-    ]
 
     # print(np.sum(gN), np.sum(gU))
     np.testing.assert_allclose(
-        np.sum([gN_ref, gU_ref], axis=0), gN.ravel() + gU.ravel(), atol=1e-2
+        np.sum([np.ravel(gN_ref), np.ravel(gU_ref)], axis=0),
+        gN.ravel() + gU.ravel(),
+        atol=1e-2,
     )
-    np.testing.assert_allclose(gN_ref, gN.ravel(), atol=1e-2)
+    np.testing.assert_allclose(np.ravel(gN_ref), gN.ravel(), atol=1e-3)
-    np.testing.assert_allclose(gU_ref, gU.ravel(), atol=1e-2)
+    np.testing.assert_allclose(np.ravel(gU_ref), gU.ravel(), atol=1e-3)
 
 
 def test_kappa_RTA_nacl(nacl_pbe: Phono3py):
     """Test RTA by NaCl."""
+    if nacl_pbe._make_r0_average:
+        ref_kappa_RTA = [7.881, 7.881, 7.881, 0, 0, 0]
+    else:
+        ref_kappa_RTA = [7.741, 7.741, 7.741, 0, 0, 0]
     kappa = _get_kappa(nacl_pbe, [9, 9, 9]).ravel()
-    np.testing.assert_allclose(nacl_pbe_kappa_RTA, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA, kappa, atol=0.5)
 
 
 def test_kappa_RTA_nacl_with_sigma(nacl_pbe: Phono3py):
     """Test RTA with smearing method by NaCl."""
+    if nacl_pbe._make_r0_average:
+        ref_kappa_RTA_with_sigma = [7.895, 7.895, 7.895, 0, 0, 0]
+    else:
+        ref_kappa_RTA_with_sigma = [7.719, 7.719, 7.719, 0, 0, 0]
     nacl_pbe.sigmas = [
         0.1,
     ]
     nacl_pbe.sigma_cutoff = 3
     kappa = _get_kappa(nacl_pbe, [9, 9, 9]).ravel()
-    np.testing.assert_allclose(nacl_pbe_kappa_RTA_with_sigma, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_with_sigma, kappa, atol=0.5)
     nacl_pbe.sigmas = None
     nacl_pbe.sigma_cutoff = None
 
 
 def test_kappa_RTA_aln(aln_lda: Phono3py):
     """Test RTA by AlN."""
+    if aln_lda._make_r0_average:
+        ref_kappa_RTA = [206.379, 206.379, 219.786, 0, 0, 0]
+    else:
+        ref_kappa_RTA = [203.278, 203.278, 212.965, 0, 0, 0]
     kappa = _get_kappa(aln_lda, [7, 7, 5]).ravel()
-    np.testing.assert_allclose(aln_lda_kappa_RTA, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA, kappa, atol=0.5)
-
-
-def test_kappa_RTA_aln_with_r0_ave(aln_lda: Phono3py):
-    """Test RTA by AlN."""
-    kappa = _get_kappa(aln_lda, [7, 7, 5], make_r0_average=True).ravel()
-    np.testing.assert_allclose(aln_lda_kappa_RTA_r0_ave, kappa, atol=0.5)
 
 
 def test_kappa_RTA_aln_with_sigma(aln_lda: Phono3py):
     """Test RTA with smearing method by AlN."""
+    if aln_lda._make_r0_average:
+        ref_kappa_RTA_with_sigmas = [217.598, 217.598, 230.099, 0, 0, 0]
+    else:
+        ref_kappa_RTA_with_sigmas = [213.820, 213.820, 224.800, 0, 0, 0]
     aln_lda.sigmas = [
         0.1,
     ]
     aln_lda.sigma_cutoff = 3
     kappa = _get_kappa(aln_lda, [7, 7, 5]).ravel()
-    np.testing.assert_allclose(aln_lda_kappa_RTA_with_sigmas, kappa, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_with_sigmas, kappa, atol=0.5)
-    aln_lda.sigmas = None
-    aln_lda.sigma_cutoff = None
-
-
-def test_kappa_RTA_aln_with_sigma_and_r0_ave(aln_lda: Phono3py):
-    """Test RTA with smearing method by AlN."""
-    aln_lda.sigmas = [
-        0.1,
-    ]
-    aln_lda.sigma_cutoff = 3
-    kappa = _get_kappa(aln_lda, [7, 7, 5], make_r0_average=True).ravel()
-    np.testing.assert_allclose(aln_lda_kappa_RTA_with_sigmas_r0_ave, kappa, atol=0.5)
     aln_lda.sigmas = None
     aln_lda.sigma_cutoff = None
 
@@ -374,13 +330,9 @@ def _get_kappa(
     is_isotope=False,
     is_full_pp=False,
     openmp_per_triplets=None,
-    make_r0_average=False,
 ):
     ph3.mesh_numbers = mesh
-    make_r0_average_orig = ph3._make_r0_average
-    ph3._make_r0_average = make_r0_average
     ph3.init_phph_interaction(openmp_per_triplets=openmp_per_triplets)
-    ph3._make_r0_average = make_r0_average_orig
     ph3.run_thermal_conductivity(
         temperatures=[
             300,

test/conductivity/test_kappa_RTA_Wigner.py

@@ -1,127 +1,132 @@
 """Test for Conductivity_RTA.py."""
+import itertools
+
 import numpy as np
+import pytest
 
-# first list is k_P, second list is k_C
+from phono3py import Phono3py
-si_pbesol_kappa_P_RTA = [108.723, 108.723, 108.723, 0.000, 0.000, 0.000]
-si_pbesol_kappa_C = [0.167, 0.167, 0.167, 0.000, 0.000, 0.000]
-
-si_pbesol_kappa_P_RTA_iso = [98.008, 98.008, 98.008, 0.000, 0.000, 0.000]
-si_pbesol_kappa_C_iso = [0.177, 0.177, 0.177, 0.000, 0.000, 0.000]
-
-si_pbesol_kappa_P_RTA_with_sigmas = [110.534, 110.534, 110.534, 0, 0, 0]
-si_pbesol_kappa_C_with_sigmas = [0.163, 0.163, 0.163, 0.000, 0.000, 0.000]
-
-si_pbesol_kappa_P_RTA_with_sigmas_iso = [97.268, 97.268, 97.268, 0, 0, 0]
-si_pbesol_kappa_C_with_sigmas_iso = [0.179, 0.179, 0.179, 0.000, 0.000, 0.000]
-
-si_pbesol_kappa_P_RTA_si_nosym = [39.325, 39.323, 39.496, -0.004, 0.020, 0.018]
-
-si_pbesol_kappa_C_si_nosym = [0.009, 0.009, 0.009, 0.000, 0.000, 0.000]
-
-si_pbesol_kappa_P_RTA_si_nomeshsym = [39.411, 39.411, 39.411, 0, 0, 0]
-si_pbesol_kappa_C_si_nomeshsym = [0.009, 0.009, 0.009, 0.000, 0.000, 0.000]
-
-nacl_pbe_kappa_P_RTA = [7.753, 7.753, 7.753, 0.000, 0.000, 0.000]
-nacl_pbe_kappa_C = [0.081, 0.081, 0.081, 0.000, 0.000, 0.000]
-
-nacl_pbe_kappa_RTA_with_sigma = [7.742, 7.742, 7.742, 0, 0, 0]  # old value 7.71913708
-nacl_pbe_kappa_C_with_sigma = [0.081, 0.081, 0.081, 0.000, 0.000, 0.000]
-
-aln_lda_kappa_P_RTA = [203.304, 203.304, 213.003, 0, 0, 0]
-aln_lda_kappa_C = [0.084, 0.084, 0.037, 0, 0, 0]
-
-aln_lda_kappa_P_RTA_with_sigmas = [213.820000, 213.820000, 224.800000, 0, 0, 0]
-aln_lda_kappa_C_with_sigmas = [0.084, 0.084, 0.036, 0, 0, 0]
 
 
-def test_kappa_RTA_si(si_pbesol):
+@pytest.mark.parametrize(
+    "is_full_pp,is_compact_fc", itertools.product([False, True], [False, True])
+)
+def test_kappa_RTA_si(
+    si_pbesol: Phono3py, si_pbesol_compact_fc: Phono3py, is_full_pp: bool, is_compact_fc
+):
     """Test RTA by Si."""
-    kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol, [9, 9, 9])
+    if is_compact_fc:
-    np.testing.assert_allclose(si_pbesol_kappa_P_RTA, kappa_P_RTA, atol=0.5)
+        ph3 = si_pbesol_compact_fc
-    np.testing.assert_allclose(si_pbesol_kappa_C, kappa_C, atol=0.02)
+    else:
+        ph3 = si_pbesol
+    if ph3._make_r0_average:
+        if is_compact_fc:
+            ref_kappa_P_RTA = [108.428, 108.428, 108.428, 0.000, 0.000, 0.000]
+            ref_kappa_C = [0.167, 0.167, 0.167, 0.000, 0.000, 0.000]
+        else:
+            ref_kappa_P_RTA = [108.527, 108.527, 108.527, 0.000, 0.000, 0.000]
+            ref_kappa_C = [0.166, 0.166, 0.166, 0.000, 0.000, 0.000]
+    else:
+        if is_compact_fc:
+            ref_kappa_P_RTA = [108.573, 108.573, 108.573, 0.000, 0.000, 0.000]
+            ref_kappa_C = [0.166, 0.166, 0.166, 0.000, 0.000, 0.000]
+        else:
+            ref_kappa_P_RTA = [108.684, 108.684, 108.684, 0.000, 0.000, 0.000]
+            ref_kappa_C = [0.166, 0.166, 0.166, 0.000, 0.000, 0.000]
+
+    kappa_P_RTA, kappa_C = _get_kappa_RTA(ph3, [9, 9, 9], is_full_pp=is_full_pp)
+    np.testing.assert_allclose(ref_kappa_P_RTA, kappa_P_RTA, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_C, kappa_C, atol=0.02)
 
 
-def test_kappa_RTA_si_full_pp(si_pbesol):
+def test_kappa_RTA_si_iso(si_pbesol: Phono3py):
-    """Test RTA with full-pp by Si."""
-    kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol, [9, 9, 9], is_full_pp=True)
-    np.testing.assert_allclose(si_pbesol_kappa_P_RTA, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(si_pbesol_kappa_C, kappa_C, atol=0.02)
-
-
-def test_kappa_RTA_si_iso(si_pbesol):
     """Test RTA with isotope scattering by Si."""
+    if si_pbesol._make_r0_average:
+        ref_kappa_P_RTA_iso = [98.026, 98.026, 98.026, 0.000, 0.000, 0.000]
+        ref_kappa_C_iso = [0.177, 0.177, 0.177, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_P_RTA_iso = [98.072, 98.072, 98.072, 0.000, 0.000, 0.000]
+        ref_kappa_C_iso = [0.177, 0.177, 0.177, 0.000, 0.000, 0.000]
+
     kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol, [9, 9, 9], is_isotope=True)
-    np.testing.assert_allclose(si_pbesol_kappa_P_RTA_iso, kappa_P_RTA, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_P_RTA_iso, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(si_pbesol_kappa_C_iso, kappa_C, atol=0.02)
+    np.testing.assert_allclose(ref_kappa_C_iso, kappa_C, atol=0.02)
 
 
-def test_kappa_RTA_si_with_sigma(si_pbesol):
+@pytest.mark.parametrize("is_full_pp", [False, True])
+def test_kappa_RTA_si_with_sigma(si_pbesol: Phono3py, is_full_pp: bool):
     """Test RTA with smearing method by Si."""
+    if si_pbesol._make_r0_average:
+        ref_kappa_P_RTA_with_sigmas = [110.857, 110.857, 110.857, 0, 0, 0]
+        ref_kappa_C_with_sigmas = [0.163, 0.163, 0.163, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_P_RTA_with_sigmas = [110.534, 110.534, 110.534, 0, 0, 0]
+        ref_kappa_C_with_sigmas = [0.162, 0.162, 0.162, 0.000, 0.000, 0.000]
+
     si_pbesol.sigmas = [
         0.1,
     ]
-    kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol, [9, 9, 9])
+    kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol, [9, 9, 9], is_full_pp=is_full_pp)
-    np.testing.assert_allclose(si_pbesol_kappa_P_RTA_with_sigmas, kappa_P_RTA, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_P_RTA_with_sigmas, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(si_pbesol_kappa_C_with_sigmas, kappa_C, atol=0.02)
+    np.testing.assert_allclose(ref_kappa_C_with_sigmas, kappa_C, atol=0.02)
     si_pbesol.sigmas = None
 
 
-def test_kappa_RTA_si_with_sigma_full_pp(si_pbesol):
+def test_kappa_RTA_si_with_sigma_iso(si_pbesol: Phono3py):
-    """Test RTA with smearing method and full-pp by Si."""
-    si_pbesol.sigmas = [
-        0.1,
-    ]
-    kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol, [9, 9, 9], is_full_pp=True)
-    np.testing.assert_allclose(si_pbesol_kappa_P_RTA_with_sigmas, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(si_pbesol_kappa_C_with_sigmas, kappa_C, atol=0.02)
-    si_pbesol.sigmas = None
-
-
-def test_kappa_RTA_si_with_sigma_iso(si_pbesol):
     """Test RTA with smearing method and isotope scattering by Si."""
+    if si_pbesol._make_r0_average:
+        ref_kappa_P_RTA_with_sigmas_iso = [97.203, 97.203, 97.203, 0, 0, 0]
+        ref_kappa_C_with_sigmas_iso = [0.176, 0.176, 0.176, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_P_RTA_with_sigmas_iso = [96.847, 96.847, 96.847, 0, 0, 0]
+        ref_kappa_C_with_sigmas_iso = [0.176, 0.176, 0.176, 0.000, 0.000, 0.000]
+
     si_pbesol.sigmas = [
         0.1,
     ]
     kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol, [9, 9, 9], is_isotope=True)
-    np.testing.assert_allclose(
+    np.testing.assert_allclose(ref_kappa_P_RTA_with_sigmas_iso, kappa_P_RTA, atol=0.5)
-        si_pbesol_kappa_P_RTA_with_sigmas_iso, kappa_P_RTA, atol=0.5
+    np.testing.assert_allclose(ref_kappa_C_with_sigmas_iso, kappa_C, atol=0.02)
-    )
-    np.testing.assert_allclose(si_pbesol_kappa_C_with_sigmas_iso, kappa_C, atol=0.02)
     si_pbesol.sigmas = None
 
 
-def test_kappa_RTA_si_compact_fc(si_pbesol_compact_fc):
+def test_kappa_RTA_si_nosym(si_pbesol: Phono3py, si_pbesol_nosym: Phono3py):
-    """Test RTA with compact-fc by Si."""
-    kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol_compact_fc, [9, 9, 9])
-    np.testing.assert_allclose(si_pbesol_kappa_P_RTA, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(si_pbesol_kappa_C, kappa_C, atol=0.02)
-
-
-def test_kappa_RTA_si_nosym(si_pbesol, si_pbesol_nosym):
     """Test RTA without considering symmetry by Si."""
+
+    if si_pbesol_nosym._make_r0_average:
+        ref_kappa_P_RTA_si_nosym = [39.396, 39.222, 39.368, -0.096, -0.022, 0.026]
+        ref_kappa_C_si_nosym = [0.009, 0.009, 0.009, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_P_RTA_si_nosym = [39.430, 39.259, 39.381, -0.096, -0.019, 0.028]
+        ref_kappa_C_si_nosym = [0.009, 0.009, 0.009, 0.000, 0.000, 0.000]
+
     si_pbesol_nosym.fc2 = si_pbesol.fc2
     si_pbesol_nosym.fc3 = si_pbesol.fc3
     kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol_nosym, [4, 4, 4])
     kappa_P_RTA_r = kappa_P_RTA.reshape(-1, 3).sum(axis=1)
     kappa_C_r = kappa_C.reshape(-1, 3).sum(axis=1)
-    kappa_P_ref = np.reshape(si_pbesol_kappa_P_RTA_si_nosym, (-1, 3)).sum(axis=1)
+    kappa_P_ref = np.reshape(ref_kappa_P_RTA_si_nosym, (-1, 3)).sum(axis=1)
-    kappa_C_ref = np.reshape(si_pbesol_kappa_C_si_nosym, (-1, 3)).sum(axis=1)
+    kappa_C_ref = np.reshape(ref_kappa_C_si_nosym, (-1, 3)).sum(axis=1)
     np.testing.assert_allclose(kappa_P_ref / 3, kappa_P_RTA_r / 3, atol=0.8)
     np.testing.assert_allclose(kappa_C_ref / 3, kappa_C_r / 3, atol=0.02)
 
 
-def test_kappa_RTA_si_nomeshsym(si_pbesol, si_pbesol_nomeshsym):
+def test_kappa_RTA_si_nomeshsym(si_pbesol: Phono3py, si_pbesol_nomeshsym: Phono3py):
     """Test RTA without considering mesh symmetry by Si."""
+    if si_pbesol_nomeshsym._make_r0_average:
+        ref_kappa_P_RTA_si_nomeshsym = [39.321, 39.321, 39.321, 0, 0, 0]
+        ref_kappa_C_si_nomeshsym = [0.009, 0.009, 0.009, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_P_RTA_si_nomeshsym = [39.332, 39.332, 39.332, 0, 0, 0]
+        ref_kappa_C_si_nomeshsym = [0.009, 0.009, 0.009, 0.000, 0.000, 0.000]
+
     si_pbesol_nomeshsym.fc2 = si_pbesol.fc2
     si_pbesol_nomeshsym.fc3 = si_pbesol.fc3
     kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol_nomeshsym, [4, 4, 4])
-    np.testing.assert_allclose(
+    np.testing.assert_allclose(ref_kappa_P_RTA_si_nomeshsym, kappa_P_RTA, atol=1.0)
-        si_pbesol_kappa_P_RTA_si_nomeshsym, kappa_P_RTA, atol=1.0
+    np.testing.assert_allclose(ref_kappa_C_si_nomeshsym, kappa_C, atol=0.02)
-    )
-    np.testing.assert_allclose(si_pbesol_kappa_C_si_nomeshsym, kappa_C, atol=0.02)
 
 
-def test_kappa_RTA_si_N_U(si_pbesol):
+def test_kappa_RTA_si_N_U(si_pbesol: Phono3py):
     """Test RTA with N and U scatterings by Si."""
     ph3 = si_pbesol
     mesh = [4, 4, 4]
@@ -137,156 +142,130 @@ def test_kappa_RTA_si_N_U(si_pbesol):
     )
     gN, gU = ph3.thermal_conductivity.get_gamma_N_U()
 
-    gN_ref = [
+    # for line in gU.reshape(-1, 6):
-        0.00000000,
+    #     print("[", ",".join([f"{val:.8f}" for val in line]), "],")
-        0.00000000,
+    if si_pbesol._make_r0_average:
-        0.00000000,
+        gN_ref = [
-        0.07402084,
+            [0.0, 0.0, 0.0, 0.07898606, 0.07898606, 0.07898606],
-        0.07402084,
+            [0.00079647, 0.00079647, 0.00913611, 0.01911102, 0.04553001, 0.04553001],
-        0.07402084,
+            [0.00173868, 0.00173868, 0.01404937, 0.00201732, 0.03354033, 0.03354033],
-        0.00078535,
+            [0.00223616, 0.00223616, 0.01039331, 0.02860916, 0.02860916, 0.01987485],
-        0.00078535,
+            [0.00291788, 0.00356241, 0.02858543, 0.00367742, 0.0206599, 0.01533763],
-        0.00917995,
+            [0.00146333, 0.00343175, 0.01596851, 0.00626596, 0.0243162, 0.01091592],
-        0.02178049,
+            [0.00396766, 0.00396766, 0.00159161, 0.00159161, 0.01479018, 0.01479018],
-        0.04470075,
+            [0.0068274, 0.0068274, 0.03983399, 0.03983399, 0.02728522, 0.02728522],
-        0.04470075,
+        ]
-        0.00173337,
+        gU_ref = [
-        0.00173337,
+            [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
-        0.01240191,
+            [0.00015184, 0.00015184, 0.00075965, 0.00736940, 0.00114177, 0.00114177],
-        0.00198981,
+            [0.00022400, 0.00022400, 0.00072237, 0.00000112, 0.00022016, 0.00022016],
-        0.03165195,
+            [0.00079188, 0.00079188, 0.00106579, 0.00418717, 0.00418717, 0.00712761],
-        0.03165195,
+            [0.00219252, 0.00262840, 0.01927670, 0.00491388, 0.01254730, 0.00519414],
-        0.00224713,
+            [0.00146999, 0.00168024, 0.01596274, 0.00641979, 0.00597353, 0.00859841],
-        0.00224713,
+            [0.00307881, 0.00307881, 0.00036554, 0.00036554, 0.01176737, 0.01176737],
-        0.00860026,
+            [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
-        0.03083611,
+        ]
-        0.03083611,
+    else:
-        0.02142118,
+        gN_ref = [
-        0.00277534,
+            [0.0, 0.0, 0.0, 0.07832198, 0.07832198, 0.07832198],
-        0.00330170,
+            [0.00079578, 0.00079578, 0.00909025, 0.01917012, 0.04557656, 0.04557656],
-        0.02727451,
+            [0.00176235, 0.00176235, 0.01414436, 0.00204092, 0.03361112, 0.03361112],
-        0.00356415,
+            [0.00221919, 0.00221919, 0.01020133, 0.02889554, 0.02889554, 0.01995543],
-        0.01847744,
+            [0.00292189, 0.00356099, 0.02855954, 0.0037053, 0.0207185, 0.01533334],
-        0.01320643,
+            [0.00147656, 0.00342335, 0.0158943, 0.00630792, 0.02427768, 0.01099287],
-        0.00155072,
+            [0.00400675, 0.00400675, 0.00162186, 0.00162186, 0.01478489, 0.01478489],
-        0.00365611,
+            [0.00676576, 0.00676576, 0.0398429, 0.0398429, 0.02715102, 0.02715102],
-        0.01641919,
+        ]
-        0.00650083,
+        gU_ref = [
-        0.02576069,
+            [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
-        0.01161589,
+            [0.00015178, 0.00015178, 0.00076936, 0.00727539, 0.00113112, 0.00113112],
-        0.00411969,
+            [0.00022696, 0.00022696, 0.00072558, 0.00000108, 0.00021968, 0.00021968],
-        0.00411969,
+            [0.00079397, 0.00079397, 0.00111068, 0.00424761, 0.00424761, 0.00697760],
-        0.00168211,
+            [0.00219456, 0.00261878, 0.01928629, 0.00490046, 0.01249235, 0.00517685],
-        0.00168211,
+            [0.00149539, 0.00161230, 0.01594274, 0.00653088, 0.00593572, 0.00849890],
-        0.01560092,
+            [0.00311169, 0.00311169, 0.00036610, 0.00036610, 0.01171667, 0.01171667],
-        0.01560092,
+            [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
-        0.00620091,
+        ]
-        0.00620091,
-        0.03764912,
-        0.03764912,
-        0.02668523,
-        0.02668523,
-    ]
-    gU_ref = [
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00015178,
-        0.00015178,
-        0.00076936,
-        0.00727539,
-        0.00113112,
-        0.00113112,
-        0.00022696,
-        0.00022696,
-        0.00072558,
-        0.00000108,
-        0.00021968,
-        0.00021968,
-        0.00079397,
-        0.00079397,
-        0.00111068,
-        0.00424761,
-        0.00424761,
-        0.00697760,
-        0.00221593,
-        0.00259510,
-        0.01996296,
-        0.00498962,
-        0.01258375,
-        0.00513825,
-        0.00148802,
-        0.00161955,
-        0.01589219,
-        0.00646134,
-        0.00577275,
-        0.00849711,
-        0.00313208,
-        0.00313208,
-        0.00036610,
-        0.00036610,
-        0.01135335,
-        0.01135335,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-        0.00000000,
-    ]
 
     # print(np.sum(gN), np.sum(gU))
     np.testing.assert_allclose(
-        np.sum([gN_ref, gU_ref], axis=0), gN.ravel() + gU.ravel(), atol=1e-2
+        np.sum([np.ravel(gN_ref), np.ravel(gU_ref)], axis=0),
+        gN.ravel() + gU.ravel(),
+        atol=1e-2,
     )
-    np.testing.assert_allclose(gN_ref, gN.ravel(), atol=1e-2)
+    np.testing.assert_allclose(np.ravel(gN_ref), gN.ravel(), atol=1e-3)
-    np.testing.assert_allclose(gU_ref, gU.ravel(), atol=1e-2)
+    np.testing.assert_allclose(np.ravel(gU_ref), gU.ravel(), atol=1e-3)
 
 
-def test_kappa_RTA_nacl(nacl_pbe):
+def test_kappa_RTA_nacl(nacl_pbe: Phono3py):
     """Test RTA by NaCl."""
+    if nacl_pbe._make_r0_average:
+        ref_kappa_P_RTA = [7.907, 7.907, 7.907, 0.000, 0.000, 0.000]
+        ref_kappa_C = [0.080, 0.080, 0.080, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_P_RTA = [7.766, 7.766, 7.766, 0.000, 0.000, 0.000]
+        ref_kappa_C = [0.080, 0.080, 0.080, 0.000, 0.000, 0.000]
+
     kappa_P_RTA, kappa_C = _get_kappa_RTA(nacl_pbe, [9, 9, 9])
-    np.testing.assert_allclose(nacl_pbe_kappa_P_RTA, kappa_P_RTA, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_P_RTA, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(nacl_pbe_kappa_C, kappa_C, atol=0.02)
+    np.testing.assert_allclose(ref_kappa_C, kappa_C, atol=0.02)
 
 
-def test_kappa_RTA_nacl_with_sigma(nacl_pbe):
+def test_kappa_RTA_nacl_with_sigma(nacl_pbe: Phono3py):
     """Test RTA with smearing method by NaCl."""
+    if nacl_pbe._make_r0_average:
+        ref_kappa_RTA_with_sigma = [7.918, 7.918, 7.918, 0, 0, 0]
+        ref_kappa_C_with_sigma = [0.081, 0.081, 0.081, 0.000, 0.000, 0.000]
+    else:
+        ref_kappa_RTA_with_sigma = [7.742, 7.742, 7.742, 0, 0, 0]
+        ref_kappa_C_with_sigma = [0.081, 0.081, 0.081, 0.000, 0.000, 0.000]
+
     nacl_pbe.sigmas = [
         0.1,
     ]
     nacl_pbe.sigma_cutoff = 3
     kappa_P_RTA, kappa_C = _get_kappa_RTA(nacl_pbe, [9, 9, 9])
-    np.testing.assert_allclose(nacl_pbe_kappa_RTA_with_sigma, kappa_P_RTA, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_RTA_with_sigma, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(nacl_pbe_kappa_C_with_sigma, kappa_C, atol=0.02)
+    np.testing.assert_allclose(ref_kappa_C_with_sigma, kappa_C, atol=0.02)
     nacl_pbe.sigmas = None
     nacl_pbe.sigma_cutoff = None
 
 
-def test_kappa_RTA_aln(aln_lda):
+def test_kappa_RTA_aln(aln_lda: Phono3py):
     """Test RTA by AlN."""
+    if aln_lda._make_r0_average:
+        ref_kappa_P_RTA = [206.379, 206.379, 219.786, 0, 0, 0]
+        ref_kappa_C = [0.083, 0.083, 0.037, 0, 0, 0]
+    else:
+        ref_kappa_P_RTA = [203.278, 203.278, 212.965, 0, 0, 0]
+        ref_kappa_C = [0.083, 0.083, 0.037, 0, 0, 0]
+
     kappa_P_RTA, kappa_C = _get_kappa_RTA(aln_lda, [7, 7, 5])
-    np.testing.assert_allclose(aln_lda_kappa_P_RTA, kappa_P_RTA, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_P_RTA, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(aln_lda_kappa_C, kappa_C, atol=0.02)
+    np.testing.assert_allclose(ref_kappa_C, kappa_C, atol=0.02)
 
 
-def test_kappa_RTA_aln_with_sigma(aln_lda):
+def test_kappa_RTA_aln_with_sigma(aln_lda: Phono3py):
     """Test RTA with smearing method by AlN."""
+    if aln_lda._make_r0_average:
+        ref_kappa_P_RTA_with_sigmas = [217.598, 217.598, 230.099, 0, 0, 0]
+        ref_kappa_C_with_sigmas = [0.084, 0.084, 0.036, 0, 0, 0]
+    else:
+        ref_kappa_P_RTA_with_sigmas = [213.820, 213.820, 224.800, 0, 0, 0]
+        ref_kappa_C_with_sigmas = [0.084, 0.084, 0.036, 0, 0, 0]
+
     aln_lda.sigmas = [
         0.1,
     ]
     aln_lda.sigma_cutoff = 3
     kappa_P_RTA, kappa_C = _get_kappa_RTA(aln_lda, [7, 7, 5])
-    np.testing.assert_allclose(aln_lda_kappa_P_RTA_with_sigmas, kappa_P_RTA, atol=0.5)
+    np.testing.assert_allclose(ref_kappa_P_RTA_with_sigmas, kappa_P_RTA, atol=0.5)
-    np.testing.assert_allclose(aln_lda_kappa_C_with_sigmas, kappa_C, atol=0.02)
+    np.testing.assert_allclose(ref_kappa_C_with_sigmas, kappa_C, atol=0.02)
     aln_lda.sigmas = None
     aln_lda.sigma_cutoff = None
 
 
-def _get_kappa_RTA(ph3, mesh, is_isotope=False, is_full_pp=False):
+def _get_kappa_RTA(ph3: Phono3py, mesh, is_isotope=False, is_full_pp=False):
     ph3.mesh_numbers = mesh
     ph3.init_phph_interaction()
     ph3.run_thermal_conductivity(

test/conftest.py

@@ -58,10 +58,11 @@ def si_pbesol(request) -> Phono3py:
     """
     yaml_filename = cwd / "phono3py_si_pbesol.yaml"
     forces_fc3_filename = cwd / "FORCES_FC3_si_pbesol"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         forces_fc3_filename=forces_fc3_filename,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -77,11 +78,12 @@ def si_pbesol_grg(request) -> Phono3py:
     """
     yaml_filename = cwd / "phono3py_si_pbesol.yaml"
     forces_fc3_filename = cwd / "FORCES_FC3_si_pbesol"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         forces_fc3_filename=forces_fc3_filename,
         use_grg=True,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -96,12 +98,13 @@ def si_pbesol_nosym(request) -> Phono3py:
     """
     yaml_filename = cwd / "phono3py_si_pbesol.yaml"
     forces_fc3_filename = cwd / "FORCES_FC3_si_pbesol"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         forces_fc3_filename=forces_fc3_filename,
         is_symmetry=False,
         produce_fc=False,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -116,12 +119,13 @@ def si_pbesol_nomeshsym(request) -> Phono3py:
     """
     yaml_filename = cwd / "phono3py_si_pbesol.yaml"
     forces_fc3_filename = cwd / "FORCES_FC3_si_pbesol"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         forces_fc3_filename=forces_fc3_filename,
         is_mesh_symmetry=False,
         produce_fc=False,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -136,11 +140,12 @@ def si_pbesol_compact_fc(request) -> Phono3py:
     """
     yaml_filename = cwd / "phono3py_si_pbesol.yaml"
     forces_fc3_filename = cwd / "FORCES_FC3_si_pbesol"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         forces_fc3_filename=forces_fc3_filename,
         is_compact_fc=True,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -154,9 +159,10 @@ def si_pbesol_111(request) -> Phono3py:
 
     """
     yaml_filename = cwd / "phono3py_params_Si111.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -173,10 +179,11 @@ def si_pbesol_111_alm(request) -> Phono3py:
     pytest.importorskip("alm")
 
     yaml_filename = cwd / "phono3py_params_Si111.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         fc_calculator="alm",
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -203,9 +210,10 @@ def si_pbesol_111_222_fd(request) -> Phono3py:
 
     """
     yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -222,10 +230,11 @@ def si_pbesol_111_222_alm(request) -> Phono3py:
     pytest.importorskip("alm")
 
     yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         fc_calculator="alm",
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -242,10 +251,11 @@ def si_pbesol_111_222_alm_fd(request) -> Phono3py:
     pytest.importorskip("alm")
 
     yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         fc_calculator="alm|",
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -262,10 +272,11 @@ def si_pbesol_111_222_fd_alm(request) -> Phono3py:
     pytest.importorskip("alm")
 
     yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         fc_calculator="|alm",
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -283,11 +294,12 @@ def si_pbesol_111_222_alm_cutoff(request) -> Phono3py:
     pytest.importorskip("alm")
 
     yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         fc_calculator="alm",
         fc_calculator_options="cutoff = 3",
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -305,11 +317,12 @@ def si_pbesol_111_222_alm_cutoff_fc2(request) -> Phono3py:
     pytest.importorskip("alm")
 
     yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         fc_calculator="alm",
         fc_calculator_options="cutoff = 3|",
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -327,11 +340,12 @@ def si_pbesol_111_222_alm_cutoff_fc3(request) -> Phono3py:
     pytest.importorskip("alm")
 
     yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         fc_calculator="alm",
         fc_calculator_options="|cutoff = 3",
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -345,9 +359,10 @@ def nacl_pbe(request) -> Phono3py:
 
     """
     yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -361,10 +376,11 @@ def nacl_pbe_compact_fc(request) -> Phono3py:
 
     """
     yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
         is_compact_fc=True,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 
@@ -377,10 +393,11 @@ def nacl_pbe_cutoff_fc3(request) -> Phono3py:
 
     """
     yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     ph3 = phono3py.load(
         yaml_filename,
         produce_fc=False,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
     forces = ph3.forces
@@ -409,10 +426,11 @@ def nacl_pbe_cutoff_fc3_all_forces(request) -> Phono3py:
 
     """
     yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     ph3 = phono3py.load(
         yaml_filename,
         produce_fc=False,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
     forces = ph3.forces
@@ -432,10 +450,11 @@ def nacl_pbe_cutoff_fc3_compact_fc(request) -> Phono3py:
 
     """
     yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     ph3 = phono3py.load(
         yaml_filename,
         produce_fc=False,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
     forces = ph3.forces
@@ -454,9 +473,10 @@ def aln_lda(request) -> Phono3py:
 
     """
     yaml_filename = cwd / "phono3py_params_AlN332.yaml.xz"
-    # enable_v2 = request.config.getoption("--v2")
+    enable_v2 = request.config.getoption("--v2")
     return phono3py.load(
         yaml_filename,
+        make_r0_average=not enable_v2,
         log_level=1,
     )
 

test/other/test_kaccum.py

@@ -5,199 +5,6 @@ from phono3py import Phono3py
 from phono3py.other.kaccum import GammaDOSsmearing, KappaDOS, get_mfp
 from phono3py.phonon.grid import get_ir_grid_points
 
-kappados_si = [
-    -0.0000002,
-    0.0000000,
-    0.0000000,
-    1.6966400,
-    2.1977566,
-    5.1814323,
-    3.3932803,
-    25.8022392,
-    15.5096766,
-    5.0899206,
-    56.6994259,
-    19.4995156,
-    6.7865608,
-    68.7759426,
-    3.2465477,
-    8.4832011,
-    72.8398965,
-    1.6583881,
-    10.1798413,
-    74.8143686,
-    0.7945952,
-    11.8764816,
-    77.2489625,
-    5.4385183,
-    13.5731219,
-    80.9162245,
-    0.5998735,
-    15.2697621,
-    81.4303646,
-    0.0000000,
-]
-mfpdos_si = [
-    0.0000000,
-    0.0000000,
-    0.0000000,
-    806.8089241,
-    33.7703552,
-    0.0225548,
-    1613.6178483,
-    45.0137786,
-    0.0103479,
-    2420.4267724,
-    53.3456168,
-    0.0106724,
-    3227.2356966,
-    62.4915811,
-    0.0107850,
-    4034.0446207,
-    69.8839011,
-    0.0075919,
-    4840.8535449,
-    74.8662085,
-    0.0049228,
-    5647.6624690,
-    78.2273252,
-    0.0035758,
-    6454.4713932,
-    80.5493065,
-    0.0020836,
-    7261.2803173,
-    81.4303646,
-    0.0000000,
-]
-gammados_si = [
-    -0.0000002,
-    0.0000000,
-    0.0000000,
-    1.6966400,
-    0.0000063,
-    0.0000149,
-    3.3932803,
-    0.0004133,
-    0.0012312,
-    5.0899206,
-    0.0071709,
-    0.0057356,
-    6.7865608,
-    0.0099381,
-    0.0006492,
-    8.4832011,
-    0.0133390,
-    0.0049604,
-    10.1798413,
-    0.0394030,
-    0.0198106,
-    11.8764816,
-    0.0495160,
-    0.0044113,
-    13.5731219,
-    0.0560223,
-    0.0050103,
-    15.2697621,
-    0.1300596,
-    0.0000000,
-]
-kappados_nacl = [
-    -0.0000002,
-    0.0000000,
-    0.0000000,
-    0.8051732,
-    0.0366488,
-    0.1820668,
-    1.6103466,
-    0.7748514,
-    1.5172957,
-    2.4155199,
-    2.0165794,
-    2.0077744,
-    3.2206933,
-    4.6670801,
-    2.8357892,
-    4.0258667,
-    6.6123781,
-    32.8560281,
-    4.8310401,
-    7.7105916,
-    0.6136893,
-    5.6362134,
-    7.9112790,
-    0.2391300,
-    6.4413868,
-    8.0272187,
-    0.0604842,
-    7.2465602,
-    8.0430831,
-    0.0000000,
-]
-mfpdos_nacl = [
-    0.0000000,
-    0.0000000,
-    0.0000000,
-    117.4892903,
-    3.1983595,
-    0.0266514,
-    234.9785806,
-    5.7974129,
-    0.0153383,
-    352.4678709,
-    7.2012603,
-    0.0075057,
-    469.9571612,
-    7.5964440,
-    0.0017477,
-    587.4464515,
-    7.7823291,
-    0.0013915,
-    704.9357418,
-    7.9195460,
-    0.0009363,
-    822.4250321,
-    8.0024702,
-    0.0004844,
-    939.9143223,
-    8.0375053,
-    0.0001382,
-    1057.4036126,
-    8.0430831,
-    0.0000000,
-]
-gammados_nacl = [
-    -0.0000002,
-    0.0000000,
-    0.0000000,
-    0.8051732,
-    0.0000822,
-    0.0004081,
-    1.6103466,
-    0.0018975,
-    0.0053389,
-    2.4155199,
-    0.0114668,
-    0.0182495,
-    3.2206933,
-    0.0353621,
-    0.0329440,
-    4.0258667,
-    0.0604996,
-    0.1138884,
-    4.8310401,
-    0.1038315,
-    0.0716216,
-    5.6362134,
-    0.1481243,
-    0.0468421,
-    6.4413868,
-    0.1982823,
-    0.0662494,
-    7.2465602,
-    0.2429551,
-    0.0000000,
-]
-
 
 def test_kappados_si(si_pbesol: Phono3py):
     """Test KappaDOS class with Si.
@@ -207,6 +14,81 @@ def test_kappados_si(si_pbesol: Phono3py):
     * kappa vs mean free path
 
     """
+    if si_pbesol._make_r0_average:
+        kappados_si = [
+            [-0.0000002, 0.0000000, 0.0000000],
+            [1.6966400, 2.1916229, 5.1669722],
+            [3.3932803, 25.7283368, 15.5208121],
+            [5.0899206, 56.6273812, 19.4749436],
+            [6.7865608, 68.6447676, 3.2126609],
+            [8.4832011, 72.6556224, 1.6322844],
+            [10.1798413, 74.6095011, 0.7885669],
+            [11.8764816, 77.0212439, 5.3742100],
+            [13.5731219, 80.7020498, 0.6098605],
+            [15.2697621, 81.2167777, 0.0000000],
+        ]
+        gammados_si = [
+            [-0.0000002, 0.0000000, 0.0000000],
+            [1.6966400, 0.0000009, 0.0000022],
+            [3.3932803, 0.0000344, 0.0000772],
+            [5.0899206, 0.0003875, 0.0002461],
+            [6.7865608, 0.0005180, 0.0000450],
+            [8.4832011, 0.0007156, 0.0002556],
+            [10.1798413, 0.0017878, 0.0008881],
+            [11.8764816, 0.0023096, 0.0001782],
+            [13.5731219, 0.0026589, 0.0003556],
+            [15.2697621, 0.0077743, 0.0000000],
+        ]
+        mfpdos_si = [
+            [0.0000000, 0.0000000, 0.0000000],
+            [806.8089241, 33.7286816, 0.0222694],
+            [1613.6178483, 44.9685295, 0.0104179],
+            [2420.4267724, 53.3416745, 0.0106935],
+            [3227.2356966, 62.5115915, 0.0108348],
+            [4034.0446207, 69.8830824, 0.0075011],
+            [4840.8535449, 74.7736977, 0.0048188],
+            [5647.6624690, 78.0965121, 0.0035467],
+            [6454.4713932, 80.3863740, 0.0020324],
+            [7261.2803173, 81.2167777, 0.0000000],
+        ]
+    else:
+        kappados_si = [
+            [-0.0000002, 0.0000000, 0.0000000],
+            [1.6966400, 2.1929621, 5.1701294],
+            [3.3932803, 25.7483415, 15.5091053],
+            [5.0899206, 56.6694055, 19.5050829],
+            [6.7865608, 68.7310303, 3.2377297],
+            [8.4832011, 72.7739143, 1.6408582],
+            [10.1798413, 74.7329367, 0.7889027],
+            [11.8764816, 77.1441825, 5.3645613],
+            [13.5731219, 80.8235276, 0.6098150],
+            [15.2697621, 81.3384416, 0.0000000],
+        ]
+        gammados_si = [
+            [-0.0000002, 0.0000000, 0.0000000],
+            [1.6966400, 0.0000009, 0.0000022],
+            [3.3932803, 0.0000346, 0.0000776],
+            [5.0899206, 0.0003887, 0.0002460],
+            [6.7865608, 0.0005188, 0.0000447],
+            [8.4832011, 0.0007153, 0.0002547],
+            [10.1798413, 0.0017871, 0.0008887],
+            [11.8764816, 0.0023084, 0.0001784],
+            [13.5731219, 0.0026578, 0.0003562],
+            [15.2697621, 0.0077689, 0.0000000],
+        ]
+        mfpdos_si = [
+            [0.0000000, 0.0000000, 0.0000000],
+            [806.8089241, 33.7483611, 0.0223642],
+            [1613.6178483, 44.9984189, 0.0104116],
+            [2420.4267724, 53.3597126, 0.0106858],
+            [3227.2356966, 62.5301968, 0.0108511],
+            [4034.0446207, 69.9297685, 0.0075505],
+            [4840.8535449, 74.8603881, 0.0048576],
+            [5647.6624690, 78.1954706, 0.0035546],
+            [6454.4713932, 80.4941074, 0.0020465],
+            [7261.2803173, 81.3384416, 0.0000000],
+        ]
+
     ph3 = si_pbesol
     ph3.mesh_numbers = [7, 7, 7]
     ph3.init_phph_interaction()
@@ -220,25 +102,27 @@ def test_kappados_si(si_pbesol: Phono3py):
     freq_points, kdos = _calculate_kappados(
         ph3, tc.mode_kappa[0], freq_points=freq_points_in
     )
-    for f, (jval, ival) in zip(freq_points, kdos):
+    # for f, (jval, ival) in zip(freq_points, kdos):
-        print("%.7f, %.7f, %.7f," % (f, jval, ival))
+    #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
     np.testing.assert_allclose(
-        kappados_si, np.vstack((freq_points, kdos.T)).T.ravel(), rtol=0, atol=0.5
+        kappados_si, np.vstack((freq_points, kdos.T)).T, rtol=0, atol=0.5
     )
 
     freq_points, kdos = _calculate_kappados(
         ph3, tc.gamma[0, :, :, :, None], freq_points=freq_points_in
     )
+    # for f, (jval, ival) in zip(freq_points, kdos):
+    #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
     np.testing.assert_allclose(
-        gammados_si, np.vstack((freq_points, kdos.T)).T.ravel(), rtol=0, atol=0.5
+        gammados_si, np.vstack((freq_points, kdos.T)).T, rtol=0, atol=1e-4
     )
 
     mfp_points_in = np.array(mfpdos_si).reshape(-1, 3)[:, 0]
     mfp_points, mfpdos = _calculate_mfpdos(ph3, mfp_points_in)
-    # for f, (jval, ival) in zip(freq_points, mfpdos):
+    # for f, (jval, ival) in zip(mfp_points, mfpdos):
-    #     print("%.7f, %.7f, %.7f," % (f, jval, ival))
+    #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
     np.testing.assert_allclose(
-        mfpdos_si, np.vstack((mfp_points, mfpdos.T)).T.ravel(), rtol=0, atol=0.5
+        mfpdos_si, np.vstack((mfp_points, mfpdos.T)).T, rtol=0, atol=0.5
     )
 
 
@@ -250,6 +134,81 @@ def test_kappados_nacl(nacl_pbe: Phono3py):
     * kappa vs mean free path
 
     """
+    if nacl_pbe._make_r0_average:
+        kappados_nacl = [
+            [-0.0000002, 0.0000000, 0.0000000],
+            [0.8051732, 0.0399444, 0.1984390],
+            [1.6103466, 0.8500862, 1.6651565],
+            [2.4155199, 2.1611612, 2.0462826],
+            [3.2206933, 4.8252014, 2.7906917],
+            [4.0258667, 6.7455774, 32.0221250],
+            [4.8310401, 7.8342086, 0.6232244],
+            [5.6362134, 8.0342122, 0.2370738],
+            [6.4413868, 8.1491279, 0.0600325],
+            [7.2465602, 8.1649079, 0.0000000],
+        ]
+        gammados_nacl = [
+            [-0.0000002, 0.0000000, 0.0000000],
+            [0.8051732, 0.0000106, 0.0000528],
+            [1.6103466, 0.0002046, 0.0004709],
+            [2.4155199, 0.0009472, 0.0012819],
+            [3.2206933, 0.0022622, 0.0016645],
+            [4.0258667, 0.0034103, 0.0054783],
+            [4.8310401, 0.0061284, 0.0029336],
+            [5.6362134, 0.0080135, 0.0019550],
+            [6.4413868, 0.0106651, 0.0046371],
+            [7.2465602, 0.0151994, 0.0000000],
+        ]
+        mfpdos_nacl = [
+            [0.0000000, 0.0000000, 0.0000000],
+            [117.4892903, 3.1996975, 0.0260716],
+            [234.9785806, 5.7553155, 0.0153949],
+            [352.4678709, 7.1540040, 0.0077163],
+            [469.9571612, 7.5793413, 0.0018178],
+            [587.4464515, 7.7799946, 0.0015717],
+            [704.9357418, 7.9439439, 0.0012052],
+            [822.4250321, 8.0613451, 0.0007915],
+            [939.9143223, 8.1309598, 0.0004040],
+            [1057.4036126, 8.1601561, 0.0001157],
+        ]
+    else:
+        kappados_nacl = [
+            [-0.0000002, 0.0000000, 0.0000000],
+            [0.8051732, 0.0367419, 0.1825292],
+            [1.6103466, 0.7836072, 1.5469421],
+            [2.4155199, 2.0449081, 2.0062821],
+            [3.2206933, 4.6731679, 2.7888186],
+            [4.0258667, 6.6041834, 32.7256000],
+            [4.8310401, 7.6993258, 0.6289821],
+            [5.6362134, 7.8997102, 0.2365916],
+            [6.4413868, 8.0146450, 0.0603293],
+            [7.2465602, 8.0305633, 0.0000000],
+        ]
+        gammados_nacl = [
+            [-0.0000002, 0.0000000, 0.0000000],
+            [0.8051732, 0.0000106, 0.0000524],
+            [1.6103466, 0.0002041, 0.0004715],
+            [2.4155199, 0.0009495, 0.0012874],
+            [3.2206933, 0.0022743, 0.0016787],
+            [4.0258667, 0.0034299, 0.0053880],
+            [4.8310401, 0.0061710, 0.0029085],
+            [5.6362134, 0.0080493, 0.0019511],
+            [6.4413868, 0.0106809, 0.0045989],
+            [7.2465602, 0.0151809, 0.0000000],
+        ]
+        mfpdos_nacl = [
+            [0.0000000, 0.0000000, 0.0000000],
+            [117.4892903, 3.2044884, 0.0265249],
+            [234.9785806, 5.8068154, 0.0153182],
+            [352.4678709, 7.1822717, 0.0071674],
+            [469.9571612, 7.5691736, 0.0017935],
+            [587.4464515, 7.7601125, 0.0014313],
+            [704.9357418, 7.9015132, 0.0009674],
+            [822.4250321, 7.9875088, 0.0005054],
+            [939.9143223, 8.0243816, 0.0001485],
+            [1057.4036126, 8.0305631, 0.0000001],
+        ]
+
     ph3 = nacl_pbe
     ph3.mesh_numbers = [7, 7, 7]
     ph3.init_phph_interaction()
@@ -264,24 +223,26 @@ def test_kappados_nacl(nacl_pbe: Phono3py):
         ph3, tc.mode_kappa[0], freq_points=freq_points_in
     )
     # for f, (jval, ival) in zip(freq_points, kdos):
-    #     print("%.7f, %.7f, %.7f," % (f, jval, ival))
+    #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
     np.testing.assert_allclose(
-        kappados_nacl, np.vstack((freq_points, kdos.T)).T.ravel(), rtol=0, atol=0.5
+        kappados_nacl, np.vstack((freq_points, kdos.T)).T, rtol=0, atol=0.5
     )
 
     freq_points, kdos = _calculate_kappados(
         ph3, tc.gamma[0, :, :, :, None], freq_points=freq_points_in
     )
+    for f, (jval, ival) in zip(freq_points, kdos):
+        print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
     np.testing.assert_allclose(
-        gammados_nacl, np.vstack((freq_points, kdos.T)).T.ravel(), rtol=0, atol=0.5
+        gammados_nacl, np.vstack((freq_points, kdos.T)).T.ravel(), rtol=0, atol=1e-4
     )
 
     mfp_points_in = np.array(mfpdos_nacl).reshape(-1, 3)[:, 0]
     mfp_points, mfpdos = _calculate_mfpdos(ph3, mfp_points_in)
-    # for f, (jval, ival) in zip(freq_points, mfpdos):
+    # for f, (jval, ival) in zip(mfp_points, mfpdos):
-    #     print("%.7f, %.7f, %.7f," % (f, jval, ival))
+    #     print("[%.7f, %.7f, %.7f]," % (f, jval, ival))
     np.testing.assert_allclose(
-        mfpdos_nacl, np.vstack((mfp_points, mfpdos.T)).T.ravel(), rtol=0, atol=0.5
+        mfpdos_nacl, np.vstack((mfp_points, mfpdos.T)).T, rtol=0, atol=0.5
     )
 
 
@@ -300,17 +261,18 @@ def test_GammaDOSsmearing(nacl_pbe: Phono3py):
         phonon_states, ir_frequencies, ir_weights, num_sampling_points=10
     )
     fpoints, gdos_vals = gdos.get_gdos()
+
     gdos_ref = [
-        [-1.4312845953710325e-07, 0.001748289450006],
+        [-1.30357573e-07, 1.74828946e-03],
-        [0.8213328685698041, 0.04545825822129761],
+        [8.21332876e-01, 4.54582590e-02],
-        [1.6426658802680678, 0.2533557541451728],
+        [1.64266588e00, 2.53356134e-01],
-        [2.463998891966331, 0.9005575010964907],
+        [2.46399889e00, 9.00558131e-01],
-        [3.285331903664595, 1.6202936411038107],
+        [3.28533190e00, 1.62029335e00],
-        [4.106664915362859, 1.9916061367478763],
+        [4.10666490e00, 1.99160666e00],
-        [4.9279979270611225, 2.5977728237205513],
+        [4.92799791e00, 2.59777233e00],
-        [5.749330938759386, 0.4504707799027985],
+        [5.74933092e00, 4.50470780e-01],
-        [6.57066395045765, 0.2936475034684396],
+        [6.57066392e00, 2.93647488e-01],
-        [7.391996962155914, 0.02869983288053483],
+        [7.39199693e00, 2.86997789e-02],
     ]
 
     np.testing.assert_allclose(