From 1529da23d1106366592d85efe00c1204d0ef5b1f Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Mon, 25 Dec 2023 15:19:55 +0900 Subject: [PATCH] Set fc3-r0-average default --- c/interaction.c | 7 - c/real_to_reciprocal.c | 171 ++++++--- phono3py/api_phono3py.py | 6 - phono3py/cui/load.py | 64 ++-- test/conductivity/test_kappa_LBTE.py | 50 ++- test/conductivity/test_kappa_LBTE_Wigner.py | 31 +- test/conductivity/test_kappa_RTA.py | 356 ++++++++---------- test/conductivity/test_kappa_RTA_Wigner.py | 373 +++++++++--------- test/conftest.py | 60 ++- test/other/test_kaccum.py | 394 +++++++++----------- 10 files changed, 739 insertions(+), 773 deletions(-) diff --git a/c/interaction.c b/c/interaction.c index 57ba1e15..bce45e9b 100644 --- a/c/interaction.c +++ b/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) { diff --git a/c/real_to_reciprocal.c b/c/real_to_reciprocal.c index 90a75f33..db54dd53 100644 --- a/c/real_to_reciprocal.c +++ b/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, - 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_legacy(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_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, - atom_triplets, openmp_at_bands); + long i, num_band; + + 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 -// 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(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_legacy(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, - atom_triplets->make_r0_average); + is_compact_fc3, atom_triplets, i, j, k, 0); + 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) { - real_to_reciprocal_elements(fc3_rec_elem, q_vecs[0], q_vecs[2], fc3, - is_compact_fc3, atom_triplets, i, j, k, - 2); - for (l = 0; l < 3; l++) { - for (m = 0; m < 3; m++) { - for (n = 0; n < 3; n++) { - // fc3_rec is stored in a way swapping jm <-> il. - fc3_rec = phonoc_complex_prod( - fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factor); - fc3_reciprocal[(j * 3 + m) * num_band * num_band + - (i * 3 + l) * num_band + k * 3 + n] = - sum_lapack_complex_double( - fc3_reciprocal[(j * 3 + m) * num_band * - num_band + - (i * 3 + l) * num_band + k * 3 + - n], - fc3_rec); - } + // fc3_rec is stored in a way swapping jm <-> il. + pre_phase_factor = get_pre_phase_factor(j, q_vecs, atom_triplets); + real_to_reciprocal_elements(fc3_rec_elem, q_vecs[0], q_vecs[2], fc3, + is_compact_fc3, atom_triplets, j, i, k, 2); + 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[m * 9 + l * 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); } } + } - real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[0], fc3, - is_compact_fc3, atom_triplets, i, j, k, - 3); - for (l = 0; l < 3; l++) { - for (m = 0; m < 3; m++) { - for (n = 0; n < 3; n++) { - // fc3_rec is stored in a way swapping kn <-> il. - fc3_rec = phonoc_complex_prod( - fc3_rec_elem[l * 9 + m * 3 + n], pre_phase_factor); - fc3_reciprocal[(k * 3 + n) * num_band * num_band + - (j * 3 + m) * num_band + i * 3 + l] = - sum_lapack_complex_double( - fc3_reciprocal[(k * 3 + n) * num_band * - num_band + - (j * 3 + m) * num_band + i * 3 + - l], - fc3_rec); - } + // fc3_rec is stored in a way swapping kn <-> il. + pre_phase_factor = get_pre_phase_factor(k, q_vecs, atom_triplets); + real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[0], fc3, + is_compact_fc3, atom_triplets, k, j, i, 3); + 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[n * 9 + m * 3 + l], 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); } } } diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index b1c0bd50..dbb4d467 100644 --- a/phono3py/api_phono3py.py +++ b/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 diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py index cffeca85..04a42847 100644 --- a/phono3py/cui/load.py +++ b/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. - supercell_matrix=None, - primitive_matrix=None, - phonon_supercell_matrix=None, - is_nac=True, - calculator=None, - unitcell=None, - supercell=None, - nac_params=None, - unitcell_filename=None, - supercell_filename=None, - born_filename=None, - forces_fc3_filename: Optional[Union[os.PathLike, Sequence]] = None, - forces_fc2_filename: Optional[Union[os.PathLike, Sequence]] = None, - fc3_filename=None, - fc2_filename=None, - fc_calculator=None, - fc_calculator_options=None, - factor=None, - produce_fc=True, - is_symmetry=True, - symmetrize_fc=True, - is_mesh_symmetry=True, - is_compact_fc=False, - use_grg=False, - symprec=1e-5, - log_level=0, + phono3py_yaml: Optional[ + Union[str, bytes, os.PathLike] + ] = None, # phono3py.yaml-like must be the first argument. + supercell_matrix: Optional[Union[Sequence, np.ndarray]] = None, + primitive_matrix: Optional[Union[Sequence, np.ndarray]] = None, + phonon_supercell_matrix: Optional[Union[Sequence, np.ndarray]] = None, + is_nac: bool = True, + calculator: Optional[str] = None, + unitcell: Optional[PhonopyAtoms] = None, + supercell: Optional[PhonopyAtoms] = None, + nac_params: Optional[dict] = None, + unitcell_filename: Optional[Union[str, bytes, os.PathLike]] = None, + supercell_filename: Optional[Union[str, bytes, os.PathLike]] = None, + born_filename: Optional[Union[str, bytes, os.PathLike]] = None, + forces_fc3_filename: Optional[Union[str, bytes, os.PathLike]] = None, + forces_fc2_filename: Optional[Union[str, bytes, os.PathLike]] = None, + fc3_filename: Optional[Union[str, bytes, os.PathLike]] = None, + fc2_filename: Optional[Union[str, bytes, os.PathLike]] = None, + fc_calculator: Optional[str] = None, + fc_calculator_options: Optional[str] = None, + factor: Optional[float] = None, + produce_fc: bool = True, + is_symmetry: bool = True, + symmetrize_fc: bool = True, + is_mesh_symmetry: bool = True, + is_compact_fc: bool = False, + 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, ) diff --git a/test/conductivity/test_kappa_LBTE.py b/test/conductivity/test_kappa_LBTE.py index 02132e8f..b7715791 100644 --- a/test/conductivity/test_kappa_LBTE.py +++ b/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) - - -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) + np.testing.assert_allclose(ref_kappa, kappa, atol=0.3) 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) diff --git a/test/conductivity/test_kappa_LBTE_Wigner.py b/test/conductivity/test_kappa_LBTE_Wigner.py index 73375593..f4b805fe 100644 --- a/test/conductivity/test_kappa_LBTE_Wigner.py +++ b/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(si_pbesol_kappa_C, kappa_C, atol=0.02) + np.testing.assert_allclose(ref_kappa_P_LBTE, kappa_P, atol=0.5) + np.testing.assert_allclose(ref_kappa_C, kappa_C, atol=0.02) -''' -#coherences conductivity is not implemented for is_reducible_collision_matrix=True, +@pytest.mark.skip( + 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) -''' diff --git a/test/conductivity/test_kappa_RTA.py b/test/conductivity/test_kappa_RTA.py index 07dcc56c..3c667d98 100644 --- a/test/conductivity/test_kappa_RTA.py +++ b/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] - -aln_lda_kappa_RTA = [203.304059, 203.304059, 213.003125, 0, 0, 0] -aln_lda_kappa_RTA_r0_ave = [2.06489355e02, 2.06489355e02, 2.19821864e02, 0, 0, 0] -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]) +@pytest.mark.parametrize( + "openmp_per_triplets,is_full_pp,is_compact_fc", + itertools.product([False, True], [False, True], [False, True]), +) 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) - - -@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) + np.testing.assert_allclose(ref_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]) -def test_kappa_RTA_si_with_sigma(si_pbesol: Phono3py, openmp_per_triplets: bool): +@pytest.mark.parametrize( + "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) - 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) + np.testing.assert_allclose(ref_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(kappa_ref, kappa, atol=0.5) + np.testing.assert_allclose(ref_kappa_RTA_si_nomeshsym, 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 = [ - 0.00000000, - 0.00000000, - 0.00000000, - 0.07402084, - 0.07402084, - 0.07402084, - 0.00078535, - 0.00078535, - 0.00917995, - 0.02178049, - 0.04470075, - 0.04470075, - 0.00173337, - 0.00173337, - 0.01240191, - 0.00198981, - 0.03165195, - 0.03165195, - 0.00224713, - 0.00224713, - 0.00860026, - 0.03083611, - 0.03083611, - 0.02142118, - 0.00277534, - 0.00330170, - 0.02727451, - 0.00356415, - 0.01847744, - 0.01320643, - 0.00155072, - 0.00365611, - 0.01641919, - 0.00650083, - 0.02576069, - 0.01161589, - 0.00411969, - 0.00411969, - 0.00168211, - 0.00168211, - 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, - ] + if si_pbesol._make_r0_average: + gN_ref = [ + [0.00000000, 0.00000000, 0.00000000, 0.07898606, 0.07898606, 0.07898606], + [0.00079647, 0.00079647, 0.00913611, 0.01911102, 0.04553001, 0.04553001], + [0.00173868, 0.00173868, 0.01404937, 0.00201732, 0.03354033, 0.03354033], + [0.00223616, 0.00223616, 0.01039331, 0.02860916, 0.02860916, 0.01987485], + [0.00291788, 0.00356241, 0.02858543, 0.00367742, 0.02065990, 0.01533763], + [0.00146333, 0.00343175, 0.01596851, 0.00626596, 0.02431620, 0.01091592], + [0.00396766, 0.00396766, 0.00159161, 0.00159161, 0.01479018, 0.01479018], + [0.00682740, 0.00682740, 0.03983399, 0.03983399, 0.02728522, 0.02728522], + ] + gU_ref = [ + [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000], + [0.00015184, 0.00015184, 0.00075965, 0.00736940, 0.00114177, 0.00114177], + [0.00022400, 0.00022400, 0.00072237, 0.00000112, 0.00022016, 0.00022016], + [0.00079188, 0.00079188, 0.00106579, 0.00418717, 0.00418717, 0.00712761], + [0.00219252, 0.00262840, 0.01927670, 0.00491388, 0.01254730, 0.00519414], + [0.00146999, 0.00168024, 0.01596274, 0.00641979, 0.00597353, 0.00859841], + [0.00307881, 0.00307881, 0.00036554, 0.00036554, 0.01176737, 0.01176737], + [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000], + ] + else: + gN_ref = [ + [0.00000000, 0.00000000, 0.00000000, 0.07832198, 0.07832198, 0.07832198], + [0.00079578, 0.00079578, 0.00909025, 0.01917012, 0.04557656, 0.04557656], + [0.00176235, 0.00176235, 0.01414436, 0.00204092, 0.03361112, 0.03361112], + [0.00221919, 0.00221919, 0.01020133, 0.02889554, 0.02889554, 0.01995543], + [0.00292189, 0.00356099, 0.02855954, 0.00370530, 0.02071850, 0.01533334], + [0.00147656, 0.00342335, 0.01589430, 0.00630792, 0.02427768, 0.01099287], + [0.00400675, 0.00400675, 0.00162186, 0.00162186, 0.01478489, 0.01478489], + [0.00676576, 0.00676576, 0.03984290, 0.03984290, 0.02715102, 0.02715102], + ] + 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.00219456, 0.00261878, 0.01928629, 0.00490046, 0.01249235, 0.00517685], + [0.00149539, 0.00161230, 0.01594274, 0.00653088, 0.00593572, 0.00849890], + [0.00311169, 0.00311169, 0.00036610, 0.00036610, 0.01171667, 0.01171667], + [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(gU_ref, gU.ravel(), atol=1e-2) + np.testing.assert_allclose(np.ravel(gN_ref), gN.ravel(), atol=1e-3) + 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) - - -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) + np.testing.assert_allclose(ref_kappa_RTA, 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) - 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) + np.testing.assert_allclose(ref_kappa_RTA_with_sigmas, 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, diff --git a/test/conductivity/test_kappa_RTA_Wigner.py b/test/conductivity/test_kappa_RTA_Wigner.py index 7e27ffac..38881d8a 100644 --- a/test/conductivity/test_kappa_RTA_Wigner.py +++ b/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 -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] +from phono3py import Phono3py -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]) - 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) + if is_compact_fc: + ph3 = si_pbesol_compact_fc + 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): - """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): +def test_kappa_RTA_si_iso(si_pbesol: Phono3py): """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(si_pbesol_kappa_C_iso, kappa_C, atol=0.02) + np.testing.assert_allclose(ref_kappa_P_RTA_iso, kappa_P_RTA, atol=0.5) + 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]) - 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) + kappa_P_RTA, kappa_C = _get_kappa_RTA(si_pbesol, [9, 9, 9], is_full_pp=is_full_pp) + np.testing.assert_allclose(ref_kappa_P_RTA_with_sigmas, kappa_P_RTA, atol=0.5) + 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): - """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): +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_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( - si_pbesol_kappa_P_RTA_with_sigmas_iso, kappa_P_RTA, atol=0.5 - ) - np.testing.assert_allclose(si_pbesol_kappa_C_with_sigmas_iso, kappa_C, atol=0.02) + np.testing.assert_allclose(ref_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) si_pbesol.sigmas = None -def test_kappa_RTA_si_compact_fc(si_pbesol_compact_fc): - """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): +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_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_C_ref = np.reshape(si_pbesol_kappa_C_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(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( - si_pbesol_kappa_P_RTA_si_nomeshsym, kappa_P_RTA, atol=1.0 - ) - np.testing.assert_allclose(si_pbesol_kappa_C_si_nomeshsym, kappa_C, atol=0.02) + np.testing.assert_allclose(ref_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) -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 = [ - 0.00000000, - 0.00000000, - 0.00000000, - 0.07402084, - 0.07402084, - 0.07402084, - 0.00078535, - 0.00078535, - 0.00917995, - 0.02178049, - 0.04470075, - 0.04470075, - 0.00173337, - 0.00173337, - 0.01240191, - 0.00198981, - 0.03165195, - 0.03165195, - 0.00224713, - 0.00224713, - 0.00860026, - 0.03083611, - 0.03083611, - 0.02142118, - 0.00277534, - 0.00330170, - 0.02727451, - 0.00356415, - 0.01847744, - 0.01320643, - 0.00155072, - 0.00365611, - 0.01641919, - 0.00650083, - 0.02576069, - 0.01161589, - 0.00411969, - 0.00411969, - 0.00168211, - 0.00168211, - 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, - ] + # for line in gU.reshape(-1, 6): + # print("[", ",".join([f"{val:.8f}" for val in line]), "],") + if si_pbesol._make_r0_average: + gN_ref = [ + [0.0, 0.0, 0.0, 0.07898606, 0.07898606, 0.07898606], + [0.00079647, 0.00079647, 0.00913611, 0.01911102, 0.04553001, 0.04553001], + [0.00173868, 0.00173868, 0.01404937, 0.00201732, 0.03354033, 0.03354033], + [0.00223616, 0.00223616, 0.01039331, 0.02860916, 0.02860916, 0.01987485], + [0.00291788, 0.00356241, 0.02858543, 0.00367742, 0.0206599, 0.01533763], + [0.00146333, 0.00343175, 0.01596851, 0.00626596, 0.0243162, 0.01091592], + [0.00396766, 0.00396766, 0.00159161, 0.00159161, 0.01479018, 0.01479018], + [0.0068274, 0.0068274, 0.03983399, 0.03983399, 0.02728522, 0.02728522], + ] + gU_ref = [ + [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000], + [0.00015184, 0.00015184, 0.00075965, 0.00736940, 0.00114177, 0.00114177], + [0.00022400, 0.00022400, 0.00072237, 0.00000112, 0.00022016, 0.00022016], + [0.00079188, 0.00079188, 0.00106579, 0.00418717, 0.00418717, 0.00712761], + [0.00219252, 0.00262840, 0.01927670, 0.00491388, 0.01254730, 0.00519414], + [0.00146999, 0.00168024, 0.01596274, 0.00641979, 0.00597353, 0.00859841], + [0.00307881, 0.00307881, 0.00036554, 0.00036554, 0.01176737, 0.01176737], + [0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000], + ] + else: + gN_ref = [ + [0.0, 0.0, 0.0, 0.07832198, 0.07832198, 0.07832198], + [0.00079578, 0.00079578, 0.00909025, 0.01917012, 0.04557656, 0.04557656], + [0.00176235, 0.00176235, 0.01414436, 0.00204092, 0.03361112, 0.03361112], + [0.00221919, 0.00221919, 0.01020133, 0.02889554, 0.02889554, 0.01995543], + [0.00292189, 0.00356099, 0.02855954, 0.0037053, 0.0207185, 0.01533334], + [0.00147656, 0.00342335, 0.0158943, 0.00630792, 0.02427768, 0.01099287], + [0.00400675, 0.00400675, 0.00162186, 0.00162186, 0.01478489, 0.01478489], + [0.00676576, 0.00676576, 0.0398429, 0.0398429, 0.02715102, 0.02715102], + ] + 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.00219456, 0.00261878, 0.01928629, 0.00490046, 0.01249235, 0.00517685], + [0.00149539, 0.00161230, 0.01594274, 0.00653088, 0.00593572, 0.00849890], + [0.00311169, 0.00311169, 0.00036610, 0.00036610, 0.01171667, 0.01171667], + [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(gU_ref, gU.ravel(), atol=1e-2) + np.testing.assert_allclose(np.ravel(gN_ref), gN.ravel(), atol=1e-3) + 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(nacl_pbe_kappa_C, kappa_C, atol=0.02) + 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_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(nacl_pbe_kappa_C_with_sigma, kappa_C, atol=0.02) + np.testing.assert_allclose(ref_kappa_RTA_with_sigma, kappa_P_RTA, atol=0.5) + 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(aln_lda_kappa_C, kappa_C, atol=0.02) + 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_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(aln_lda_kappa_C_with_sigmas, kappa_C, atol=0.02) + np.testing.assert_allclose(ref_kappa_P_RTA_with_sigmas, kappa_P_RTA, atol=0.5) + 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( diff --git a/test/conftest.py b/test/conftest.py index cca29787..3bb8515d 100644 --- a/test/conftest.py +++ b/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, ) diff --git a/test/other/test_kaccum.py b/test/other/test_kaccum.py index 7ab33c9f..7b0f1586 100644 --- a/test/other/test_kaccum.py +++ b/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): - print("%.7f, %.7f, %.7f," % (f, jval, ival)) + # for f, (jval, ival) in zip(freq_points, kdos): + # 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): - # print("%.7f, %.7f, %.7f," % (f, jval, ival)) + # for f, (jval, ival) in zip(mfp_points, mfpdos): + # 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): - # print("%.7f, %.7f, %.7f," % (f, jval, ival)) + # for f, (jval, ival) in zip(mfp_points, mfpdos): + # 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], - [0.8213328685698041, 0.04545825822129761], - [1.6426658802680678, 0.2533557541451728], - [2.463998891966331, 0.9005575010964907], - [3.285331903664595, 1.6202936411038107], - [4.106664915362859, 1.9916061367478763], - [4.9279979270611225, 2.5977728237205513], - [5.749330938759386, 0.4504707799027985], - [6.57066395045765, 0.2936475034684396], - [7.391996962155914, 0.02869983288053483], + [-1.30357573e-07, 1.74828946e-03], + [8.21332876e-01, 4.54582590e-02], + [1.64266588e00, 2.53356134e-01], + [2.46399889e00, 9.00558131e-01], + [3.28533190e00, 1.62029335e00], + [4.10666490e00, 1.99160666e00], + [4.92799791e00, 2.59777233e00], + [5.74933092e00, 4.50470780e-01], + [6.57066392e00, 2.93647488e-01], + [7.39199693e00, 2.86997789e-02], ] np.testing.assert_allclose(