scnc
/
phono3py
mirror of https://github.com/phonopy/phono3py.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #201 from phonopy/r0-ave-default 

Set fc3-r0-average default
This commit is contained in:
Atsushi Togo 2023-12-25 18:04:08 +09:00 committed by GitHub
parent 7f55f9bee9 37eb8bb74c
commit d946de7e01
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
15 changed files with 3669 additions and 2346 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
.github/workflows/phono3py-pytest-conda-mkl-v2.yml vendored Normal file
View File

@ -0,0 +1,51 @@
name: Pytest with mkl and --v2 option
on:
pull_request:
branches: [ develop ]
push:
branches-ignore:
- publish-gh-pages
- master
- rc
jobs:
build-linux:
runs-on: ubuntu-latest
defaults:
run:
shell: bash -l {0}
strategy:
matrix:
python-version: ["3.11"]
steps:
- uses: actions/checkout@v3
# Use conda-incubator/setup-miniconda for precise control of conda infrastructure
- uses: conda-incubator/setup-miniconda@v2
with:
miniforge-version: latest
- name: Install dependent packages
run: |
conda activate test
conda install --yes python=${{ matrix.python-version }}
#conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py scipy pytest codecov pytest-cov spglib alm cmake c-compiler
conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest codecov pytest-cov spglib alm cmake c-compiler
- name: Install phonopy develop branch
run: |
conda activate test
git clone --depth 1 https://github.com/phonopy/phonopy.git
cd phonopy
pip install -e . -vvv
cd ..
- name: Install phono3py
run: |
conda activate test
pip install -e . -vvv
- name: Run pytest
run: |
pytest --v2 -v --cov=./ --cov-report=xml test
- name: Upload coverage to Codecov
uses: codecov/codecov-action@v3
with:
verbose: true

7
c/interaction.c
View File

@ -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) {

127
c/real_to_reciprocal.c
View File

@ -43,8 +43,15 @@
#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,
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);
@ -70,18 +77,27 @@ 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;
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,
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) {
@ -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,50 +184,44 @@ 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.
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, i, j, k,
2);
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 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] =
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[(j * 3 + m) * num_band *
num_band +
(i * 3 + l) * num_band + k * 3 +
n],
fc3_reciprocal[(i * 3 + l) * num_band * num_band +
(j * 3 + m) * num_band + k * 3 + n],
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, i, j, k,
3);
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 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] =
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[(k * 3 + n) * num_band *
num_band +
(j * 3 + m) * num_band + i * 3 +
l],
fc3_reciprocal[(i * 3 + l) * num_band * num_band +
(j * 3 + m) * num_band + k * 3 + n],
fc3_rec);
}
}
}
}
}
}
static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,

6
phono3py/api_phono3py.py
View File

@ -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

64
phono3py/cui/load.py
View File

@ -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,
)

50
test/conductivity/test_kappa_LBTE.py
View File

@ -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)

31
test/conductivity/test_kappa_LBTE_Wigner.py
View File

@ -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)
'''

348
test/conductivity/test_kappa_RTA.py
View File

@ -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,106 @@ 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]), "],")
if si_pbesol._make_r0_average:
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,
[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.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,
[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.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.sum(gN_ref, axis=1), gN[0, 0].sum(axis=1), atol=0.05)
np.testing.assert_allclose(np.sum(gU_ref, axis=1), gU[0, 0].sum(axis=1), atol=0.05)
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 +324,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,

365
test/conductivity/test_kappa_RTA_Wigner.py
View File

@ -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,124 @@ def test_kappa_RTA_si_N_U(si_pbesol):
)
gN, gU = ph3.thermal_conductivity.get_gamma_N_U()
# 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.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,
[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.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,
[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.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.sum(gN_ref, axis=1), gN[0, 0].sum(axis=1), atol=0.05)
np.testing.assert_allclose(np.sum(gU_ref, axis=1), gU[0, 0].sum(axis=1), atol=0.05)
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(

60
test/conftest.py
View File

@ -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,
)

394
test/other/test_kaccum.py
View File

@ -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, 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(

2494
test/phonon3/test_imag_self_energy.py
View File

File diff suppressed because it is too large Load Diff

204
test/phonon3/test_interaction.py
View File

@ -11,7 +11,12 @@ from phonopy.structure.cells import get_smallest_vectors
from phono3py import Phono3py
from phono3py.phonon3.interaction import Interaction
itr_RTA_Si = [
@pytest.mark.parametrize("lang", ["C", "Python"])
def test_interaction_RTA_si(si_pbesol: Phono3py, lang: Literal["C", "Python"]):
"""Test interaction_strength of Si."""
if si_pbesol._make_r0_average:
ref_itr_RTA_Si = [
4.522052e-08,
4.896362e-08,
4.614211e-08,
@ -22,9 +27,35 @@ itr_RTA_Si = [
4.645423e-08,
4.659572e-08,
4.730222e-08,
]
]
else:
ref_itr_RTA_Si = [
4.522052e-08,
4.896362e-08,
4.614211e-08,
4.744361e-08,
4.832248e-08,
4.698535e-08,
4.597876e-08,
4.645423e-08,
4.659572e-08,
4.730222e-08,
]
itr_RTA_AlN = [
itr = _get_irt(si_pbesol, [4, 4, 4])
itr.set_grid_point(1)
itr.run(lang=lang)
# _show(itr)
# (10, 6, 6, 6)
np.testing.assert_allclose(
itr.interaction_strength.sum(axis=(1, 2, 3)), ref_itr_RTA_Si, rtol=0, atol=1e-10
)
def test_interaction_RTA_AlN(aln_lda: Phono3py):
"""Test interaction_strength of AlN."""
if aln_lda._make_r0_average:
ref_itr_RTA_AlN = [
7.456796e-08,
7.242121e-08,
7.068141e-08,
@ -89,109 +120,82 @@ itr_RTA_AlN = [
6.710084e-08,
6.698233e-08,
6.694871e-08,
]
]
else:
ref_itr_RTA_AlN = [
7.456796e-08,
7.242121e-08,
7.068141e-08,
7.059521e-08,
7.289497e-08,
7.127172e-08,
7.082734e-08,
7.394367e-08,
7.084351e-08,
7.083299e-08,
7.085792e-08,
7.124150e-08,
7.048386e-08,
7.062840e-08,
7.036795e-08,
7.043995e-08,
7.366440e-08,
7.136803e-08,
6.988469e-08,
6.989518e-08,
7.179516e-08,
7.038043e-08,
7.011416e-08,
7.278196e-08,
6.999028e-08,
7.009615e-08,
7.018236e-08,
7.025054e-08,
6.977425e-08,
6.993095e-08,
6.962119e-08,
6.964423e-08,
7.121739e-08,
6.939940e-08,
6.834705e-08,
6.847351e-08,
6.977063e-08,
6.872065e-08,
6.863218e-08,
7.055696e-08,
6.836064e-08,
6.854052e-08,
6.864199e-08,
6.849059e-08,
6.826958e-08,
6.837379e-08,
6.808307e-08,
6.804480e-08,
6.961289e-08,
6.816170e-08,
6.730028e-08,
6.746055e-08,
6.851460e-08,
6.764892e-08,
6.754060e-08,
6.913662e-08,
6.729303e-08,
6.736722e-08,
6.734663e-08,
6.743441e-08,
6.713107e-08,
6.710084e-08,
6.698233e-08,
6.694871e-08,
]
itr_RTA_AlN_r0_ave = [
7.451662e-08,
7.248965e-08,
7.068341e-08,
7.038257e-08,
7.291756e-08,
7.130737e-08,
7.073777e-08,
7.391677e-08,
7.086006e-08,
7.077947e-08,
7.084617e-08,
7.128499e-08,
7.036749e-08,
7.057057e-08,
7.032233e-08,
7.041152e-08,
7.366076e-08,
7.149180e-08,
6.992398e-08,
6.972641e-08,
7.187472e-08,
7.046525e-08,
7.006265e-08,
7.280888e-08,
7.005107e-08,
7.008412e-08,
7.019342e-08,
7.034898e-08,
6.969589e-08,
6.990595e-08,
6.961230e-08,
6.966506e-08,
7.121203e-08,
6.954065e-08,
6.839930e-08,
6.831861e-08,
6.986872e-08,
6.882433e-08,
6.859035e-08,
7.059587e-08,
6.844120e-08,
6.854220e-08,
6.865246e-08,
6.861039e-08,
6.819986e-08,
6.835737e-08,
6.808818e-08,
6.808952e-08,
6.954409e-08,
6.824539e-08,
6.732596e-08,
6.726087e-08,
6.854755e-08,
6.770693e-08,
6.745257e-08,
6.910600e-08,
6.733516e-08,
6.731725e-08,
6.731323e-08,
6.749548e-08,
6.700802e-08,
6.703569e-08,
6.694826e-08,
6.694365e-08,
]
@pytest.mark.parametrize("lang", ["C", "Python"])
def test_interaction_RTA_si(si_pbesol: Phono3py, lang: Literal["C", "Python"]):
"""Test interaction_strength of Si."""
itr = _get_irt(si_pbesol, [4, 4, 4])
itr.set_grid_point(1)
itr.run(lang=lang)
# _show(itr)
# (10, 6, 6, 6)
np.testing.assert_allclose(
itr.interaction_strength.sum(axis=(1, 2, 3)), itr_RTA_Si, rtol=0, atol=1e-10
)
def test_interaction_RTA_AlN(aln_lda: Phono3py):
"""Test interaction_strength of AlN."""
itr = _get_irt(aln_lda, [7, 7, 7])
itr.set_grid_point(1)
itr.run()
# _show(itr)
np.testing.assert_allclose(
itr.interaction_strength.sum(axis=(1, 2, 3)), itr_RTA_AlN, rtol=0, atol=1e-10
)
def test_interaction_RTA_AlN_r0_ave(aln_lda: Phono3py):
"""Test interaction_strength of AlN."""
itr = _get_irt(aln_lda, [7, 7, 7], make_r0_average=True)
itr.set_grid_point(1)
itr.run()
# _show(itr)
_show(itr)
np.testing.assert_allclose(
itr.interaction_strength.sum(axis=(1, 2, 3)),
itr_RTA_AlN_r0_ave,
ref_itr_RTA_AlN,
rtol=0,
atol=1e-10,
)

493
test/phonon3/test_real_self_energy.py
View File

@ -4,12 +4,128 @@ import numpy as np
from phono3py import Phono3py
from phono3py.phonon3.real_self_energy import ImagToReal
si_pbesol_Delta = [
[-0.0057666, -0.0057666, -0.01639729, -0.14809965, -0.15091765, -0.15091765],
[-0.02078728, -0.02102094, -0.06573269, -0.11432603, -0.1366966, -0.14371315],
]
si_pbesol_Delta_fps = [
def test_real_self_energy_with_band_indices(si_pbesol: Phono3py):
"""Real part of self energy spectrum of Si.
* at frequencies of band indices.
"""
if si_pbesol._make_r0_average:
ref_Delta = [
[
-0.0057326,
-0.0057326,
-0.01633157,
-0.14797815,
-0.15086179,
-0.15086179,
],
[
-0.02080177,
-0.02104276,
-0.06575259,
-0.11439863,
-0.13666415,
-0.14372152,
],
]
else:
ref_Delta = [
[
-0.0057666,
-0.0057666,
-0.01639729,
-0.14809965,
-0.15091765,
-0.15091765,
],
[
-0.02078728,
-0.02102094,
-0.06573269,
-0.11432603,
-0.1366966,
-0.14371315,
],
]
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
_, delta = si_pbesol.run_real_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
],
write_hdf5=False,
frequency_points_at_bands=True,
)
# print(delta[0, 0, :])
np.testing.assert_allclose(ref_Delta, delta[0, 0, :], atol=0.01)
def test_real_self_energy_with_frequency_points(si_pbesol: Phono3py):
"""Real part of self energy spectrum of Si.
* specified frquency points
"""
if si_pbesol._make_r0_average:
ref_Delta_fps = [
[
-0.00573260,
-0.00589756,
-0.00839449,
-0.00966402,
-0.00573260,
-0.00589756,
-0.00839449,
-0.00966402,
-0.01486569,
-0.01633157,
-0.01993995,
-0.02073417,
-0.15512961,
-0.14758526,
-0.14797816,
-0.14222444,
-0.15683457,
-0.15691468,
-0.15974324,
-0.15086181,
-0.15683457,
-0.15691468,
-0.15974324,
-0.15086181,
],
[
-0.01990876,
-0.02077841,
-0.01798683,
-0.01933049,
-0.02158770,
-0.02193100,
-0.02186381,
-0.01877096,
-0.05743552,
-0.05244391,
-0.06241623,
-0.05643933,
-0.13068294,
-0.11925247,
-0.13470590,
-0.13105466,
-0.15192994,
-0.14204618,
-0.14364976,
-0.14173656,
-0.14762541,
-0.13904288,
-0.14284911,
-0.14112034,
],
]
else:
ref_Delta_fps = [
[
-0.00576660,
-0.00594616,
@ -62,10 +178,219 @@ si_pbesol_Delta_fps = [
-0.14275290,
-0.14100562,
],
]
]
# imag-self-energy Si-PBEsol 50x50x50 gp=5, bi=4, 101 points, 300K
im_part = [
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
frequency_points = [1.469947, 3.085309, 14.997187, 15.129080]
fps, delta = si_pbesol.run_real_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
],
frequency_points=frequency_points,
write_hdf5=False,
frequency_points_at_bands=False,
)
# for b in delta[0, 0, 0]:
# print("".join(f"{s:.8f}, " % s for s in b))
# for b in delta[0, 0, 1]:
# print("".join(f"{s:.8f}, " % s for s in b))
np.testing.assert_allclose(frequency_points, fps, atol=1e-5)
np.testing.assert_allclose(ref_Delta_fps[0], delta[0, 0, 0].ravel(), atol=0.01)
np.testing.assert_allclose(ref_Delta_fps[1], delta[0, 0, 1].ravel(), atol=0.01)
def test_real_self_energy_nacl_nac_npoints(nacl_pbe: Phono3py):
"""Real part of self energy spectrum of NaCl.
* at 10 frequency points sampled uniformly.
* at q->0
"""
ref_freq_points_nacl_nac = [
0.0,
1.63223063,
3.26446125,
4.89669188,
6.5289225,
8.16115313,
9.79338375,
11.42561438,
13.057845,
14.69007563,
]
if nacl_pbe._make_r0_average:
ref_delta_nacl_nac = [
[
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
],
[
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
],
[
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
],
[
0.00000000,
-0.43955299,
-0.24659961,
-0.45227602,
-0.80296691,
-0.35433156,
0.54514663,
0.29417714,
0.16346531,
0.11216578,
],
[
0.00000000,
-0.43955299,
-0.24659961,
-0.45227602,
-0.80296691,
-0.35433156,
0.54514663,
0.29417714,
0.16346531,
0.11216578,
],
[
0.00000000,
-0.27249573,
-0.15287654,
-0.28038322,
-0.49778992,
-0.21966368,
0.33795727,
0.18237167,
0.10133841,
0.06953586,
],
]
else:
ref_delta_nacl_nac = [
[
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
],
[
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
],
[
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
],
[
0.00000000,
-0.43835860,
-0.24732790,
-0.45125795,
-0.80781808,
-0.36217377,
0.54848873,
0.29335955,
0.16329652,
0.11208513,
],
[
0.00000000,
-0.43835860,
-0.24732790,
-0.45125795,
-0.80781808,
-0.36217377,
0.54848873,
0.29335955,
0.16329652,
0.11208513,
],
[
0.00000000,
-0.27175528,
-0.15332803,
-0.27975208,
-0.50079735,
-0.22452538,
0.34002916,
0.18186482,
0.10123376,
0.06948586,
],
]
nacl_pbe.mesh_numbers = [9, 9, 9]
nacl_pbe.init_phph_interaction(nac_q_direction=[1, 0, 0])
fps, delta = nacl_pbe.run_real_self_energy(
[nacl_pbe.grid.gp_Gamma], [300], num_frequency_points=10
)
# for line in delta[0, 0, 0]:
# print("[", ",".join([f"{val:.8f}" for val in line]), "],")
np.testing.assert_allclose(ref_freq_points_nacl_nac, fps.ravel(), rtol=0, atol=1e-5)
np.testing.assert_allclose(ref_delta_nacl_nac, delta[0, 0, 0], rtol=0, atol=1e-2)
def test_ImagToReal():
"""Test ImagToReal class (Kramers–Kronig relation)."""
# imag-self-energy Si-PBEsol 50x50x50 gp=5, bi=4, 101 points, 300K
im_part = [
[0.0000000, 0.0000000, -0.0000000, -0.1750223, 0.1534611, -0.1794583],
[0.3069222, 0.0180686, 0.3069222, -0.1692049, 0.4603833, -0.1636592],
[0.6138444, 0.0158916, 0.6138444, -0.1621245, 0.7673055, -0.1625514],
@ -167,158 +492,8 @@ im_part = [
[30.0783775, 0.1397870, 30.0783775, 0.5578788, 30.2318386, 0.5591731],
[30.3852997, 0.0468188, 30.3852997, 0.4941608, 30.5387608, 0.4706733],
[30.6922219, 0.0000000, 30.6922219, 0.4173657, 30.8456830, 0.3821416],
]
delta_nacl_nac = [
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
-0.43835832,
-0.43835832,
-0.27175512,
0.00000000,
0.00000000,
0.00000000,
-0.24732785,
-0.24732785,
-0.15332798,
0.00000000,
0.00000000,
0.00000000,
-0.45125778,
-0.45125778,
-0.27975198,
0.00000000,
0.00000000,
0.00000000,
-0.80781739,
-0.80781739,
-0.50079699,
0.00000000,
0.00000000,
0.00000000,
-0.36217245,
-0.36217245,
-0.22452456,
0.00000000,
0.00000000,
0.00000000,
0.54848940,
0.54848940,
0.34002959,
0.00000000,
0.00000000,
0.00000000,
0.29335952,
0.29335952,
0.18186480,
0.00000000,
0.00000000,
0.00000000,
0.16329650,
0.16329650,
0.10123375,
0.00000000,
0.00000000,
0.00000000,
0.11208511,
0.11208511,
0.06948585,
]
freq_points_nacl_nac = [
0.0,
1.63223063,
3.26446125,
4.89669188,
6.5289225,
8.16115313,
9.79338375,
11.42561438,
13.057845,
14.69007563,
]
]
def test_real_self_energy_with_band_indices(si_pbesol):
"""Real part of self energy spectrum of Si.
* at frequencies of band indices.
"""
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
_, delta = si_pbesol.run_real_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
],
write_hdf5=False,
frequency_points_at_bands=True,
)
np.testing.assert_allclose(si_pbesol_Delta, delta[0, 0, :], atol=0.01)
def test_real_self_energy_with_frequency_points(si_pbesol: Phono3py):
"""Real part of self energy spectrum of Si.
* specified frquency points
"""
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
frequency_points = [1.469947, 3.085309, 14.997187, 15.129080]
fps, delta = si_pbesol.run_real_self_energy(
si_pbesol.grid.grg2bzg[[1, 103]],
[
300,
],
frequency_points=frequency_points,
write_hdf5=False,
frequency_points_at_bands=False,
)
np.testing.assert_allclose(frequency_points, fps, atol=1e-5)
np.testing.assert_allclose(
si_pbesol_Delta_fps[0], delta[0, 0, 0].ravel(), atol=0.01
)
np.testing.assert_allclose(
si_pbesol_Delta_fps[1], delta[0, 0, 1].ravel(), atol=0.01
)
# for b in delta[0, 0, 0]:
# print("".join("%.8f, " % s for s in b))
# for b in delta[0, 1, 0]:
# print("".join("%.8f, " % s for s in b))
def test_real_self_energy_nacl_nac_npoints(nacl_pbe: Phono3py):
"""Real part of self energy spectrum of NaCl.
* at 10 frequency points sampled uniformly.
* at q->0
"""
nacl_pbe.mesh_numbers = [9, 9, 9]
nacl_pbe.init_phph_interaction(nac_q_direction=[1, 0, 0])
fps, delta = nacl_pbe.run_real_self_energy(
[nacl_pbe.grid.gp_Gamma], [300], num_frequency_points=10
)
# for line in np.swapaxes(delta, -1, -2).ravel().reshape(-1, 6):
# print(("%10.8f, " * 6) % tuple(line))
# print(fps.ravel())
np.testing.assert_allclose(freq_points_nacl_nac, fps.ravel(), rtol=0, atol=1e-5)
np.testing.assert_allclose(
delta_nacl_nac, np.swapaxes(delta, -1, -2).ravel(), rtol=0, atol=1e-2
)
def test_ImagToReal():
"""Test ImagToReal class (Kramers–Kronig relation)."""
vals = np.array(im_part)
i2r = ImagToReal(vals[:, 1], vals[:, 0])
i2r.run()

913
test/phonon3/test_spectral_function.py
View File

@ -1,289 +1,619 @@
"""Test spectral_function.py."""
import numpy as np
import pytest
from phono3py import Phono3py
from phono3py.phonon3.spectral_function import SpectralFunction
shifts = [
-0.0049592,
-0.0049592,
-0.0120983,
-0.1226471,
-0.1214069,
-0.1214069,
-0.0051678,
-0.0051678,
-0.0128471,
-0.1224616,
-0.1200362,
-0.1200362,
-0.0055308,
-0.0055308,
-0.0122157,
-0.1093754,
-0.1077399,
-0.1077399,
-0.0037992,
-0.0037992,
-0.0089979,
-0.0955525,
-0.0958995,
-0.0958995,
-0.0034397,
-0.0034397,
-0.0107575,
-0.1068741,
-0.1067815,
-0.1067815,
-0.0017800,
-0.0017800,
-0.0102865,
-0.1348585,
-0.1275650,
-0.1275650,
0.0006728,
0.0006728,
-0.0065349,
-0.2011702,
-0.2015991,
-0.2015991,
0.0021133,
0.0021133,
0.0020353,
-0.0740009,
-0.0833644,
-0.0833644,
0.0037739,
0.0037739,
0.0121357,
0.1597195,
0.1585307,
0.1585307,
0.0026257,
0.0026257,
0.0103523,
0.1626420,
0.1634832,
0.1634832,
-0.0189694,
-0.0188985,
-0.0415773,
-0.0955391,
-0.1180182,
-0.1126508,
-0.0194533,
-0.0191057,
-0.0420358,
-0.0913521,
-0.1140995,
-0.1075009,
-0.0233933,
-0.0219600,
-0.0466734,
-0.0865867,
-0.1086070,
-0.1014454,
-0.0140271,
-0.0150165,
-0.0344515,
-0.0755416,
-0.1018518,
-0.0951606,
-0.0058780,
-0.0089457,
-0.0256867,
-0.0775726,
-0.1070427,
-0.1018654,
-0.0069737,
-0.0092857,
-0.0333909,
-0.1014042,
-0.1320678,
-0.1288315,
-0.0030075,
-0.0060858,
-0.0245855,
-0.1186313,
-0.1963719,
-0.1857004,
0.0058243,
0.0030539,
-0.0049966,
-0.0583228,
-0.0921850,
-0.0893692,
0.0141517,
0.0149365,
0.0312156,
0.0898626,
0.1454759,
0.1347802,
0.0110954,
0.0137260,
0.0427527,
0.1280421,
0.1715647,
0.1648037,
]
spec_funcs = [
-0.0000000,
-0.0000000,
0.0000000,
-0.0000000,
-0.0000000,
-0.0000000,
0.0000165,
0.0000165,
0.0022357,
0.0001249,
0.0001318,
0.0001318,
0.0000027,
0.0000027,
0.0000592,
0.0005056,
0.0004722,
0.0004722,
0.0000016,
0.0000016,
0.0000118,
0.0008887,
0.0008281,
0.0008281,
0.0000001,
0.0000001,
0.0000007,
0.0062806,
0.0053809,
0.0053809,
0.0000003,
0.0000003,
0.0000026,
0.0025466,
0.0035400,
0.0035400,
0.0000001,
0.0000001,
0.0000010,
0.0004469,
0.0002613,
0.0002613,
0.0000001,
0.0000001,
0.0000011,
0.0011874,
0.0012518,
0.0012518,
0.0000000,
0.0000000,
0.0000003,
0.0002964,
0.0003119,
0.0003119,
0.0000000,
0.0000000,
0.0000000,
0.0000000,
0.0000000,
0.0000000,
0.0000000,
-0.0000000,
0.0000000,
0.0000000,
0.0000000,
0.0000000,
0.0032185,
0.0005560,
0.0001434,
0.0001868,
0.0001602,
0.0001670,
0.0000726,
0.0002220,
0.0016122,
0.0004059,
0.0003461,
0.0003024,
0.0000878,
0.0001280,
0.0033720,
0.0027162,
0.0013271,
0.0011427,
0.0000049,
0.0000094,
0.0000892,
0.0137571,
0.0405108,
0.0173674,
0.0000052,
0.0000084,
0.0000567,
0.0005807,
0.0013438,
0.0014025,
0.0000021,
0.0000034,
0.0000405,
0.0003176,
0.0003022,
0.0003523,
0.0000019,
0.0000034,
0.0000235,
0.0003419,
0.0008701,
0.0008436,
0.0000004,
0.0000010,
0.0000162,
0.0001759,
0.0002910,
0.0002978,
0.0000000,
0.0000000,
0.0000000,
0.0000000,
0.0000000,
0.0000000,
]
def test_SpectralFunction(si_pbesol: Phono3py):
@pytest.mark.parametrize("check_bands", [False, True])
def test_SpectralFunction(si_pbesol: Phono3py, check_bands: bool):
"""Spectral function of Si."""
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
sf = SpectralFunction(
si_pbesol.phph_interaction,
si_pbesol.grid.grg2bzg[[1, 103]],
temperatures=[
300,
if si_pbesol._make_r0_average:
ref_shifts = [
[
[
-0.00491549,
-0.00513379,
-0.00549296,
-0.00377646,
-0.00342792,
-0.00178271,
0.00063823,
0.00208652,
0.00376214,
0.00262103,
],
num_frequency_points=10,
log_level=1,
)
sf.run()
[
-0.00491549,
-0.00513379,
-0.00549296,
-0.00377646,
-0.00342792,
-0.00178271,
0.00063823,
0.00208652,
0.00376214,
0.00262103,
],
[
-0.01200352,
-0.01272740,
-0.01213832,
-0.00902171,
-0.01079296,
-0.01033632,
-0.00661888,
0.00196845,
0.01212704,
0.01037479,
],
[
-0.12224881,
-0.12186476,
-0.10950402,
-0.09619987,
-0.10705652,
-0.13503653,
-0.20145401,
-0.07468261,
0.15954608,
0.16311580,
],
[
-0.12117175,
-0.11966157,
-0.10786945,
-0.09635481,
-0.10686396,
-0.12762500,
-0.20163789,
-0.08374266,
0.15834565,
0.16375939,
],
[
-0.12117175,
-0.11966157,
-0.10786945,
-0.09635481,
-0.10686396,
-0.12762500,
-0.20163789,
-0.08374266,
0.15834565,
0.16375939,
],
],
[
[
-0.01895098,
-0.01944799,
-0.02343492,
-0.01391542,
-0.00562335,
-0.00675556,
-0.00299433,
0.00568717,
0.01407410,
0.01103152,
],
[
-0.01883020,
-0.01904414,
-0.02189420,
-0.01497364,
-0.00887001,
-0.00907876,
-0.00600166,
0.00297018,
0.01485332,
0.01364674,
],
[
-0.04158644,
-0.04204489,
-0.04660642,
-0.03449069,
-0.02581911,
-0.03329617,
-0.02458896,
-0.00519089,
0.03119712,
0.04283032,
],
[
-0.09542060,
-0.09131037,
-0.08642350,
-0.07536168,
-0.07744657,
-0.10100751,
-0.11855897,
-0.05851935,
0.08976609,
0.12789279,
],
[
-0.11827973,
-0.11431900,
-0.10875016,
-0.10196855,
-0.10719767,
-0.13236285,
-0.19702133,
-0.09252739,
0.14584477,
0.17195801,
],
[
-0.11277741,
-0.10757643,
-0.10150954,
-0.09524697,
-0.10202478,
-0.12884057,
-0.18556233,
-0.08932224,
0.13482432,
0.16483631,
],
],
]
ref_spec_funcs = [
[
[
0.00000000,
0.00001578,
0.00000272,
0.00000161,
0.00000013,
0.00000032,
0.00000006,
0.00000006,
0.00000001,
0.00000000,
],
[
0.00000000,
0.00001578,
0.00000272,
0.00000161,
0.00000013,
0.00000032,
0.00000006,
0.00000006,
0.00000001,
0.00000000,
],
[
0.00000000,
0.00223977,
0.00005791,
0.00001160,
0.00000066,
0.00000259,
0.00000098,
0.00000112,
0.00000029,
0.00000000,
],
[
0.00000000,
0.00012503,
0.00049470,
0.00089500,
0.00632551,
0.00256627,
0.00044659,
0.00118653,
0.00029733,
0.00000000,
],
[
0.00000000,
0.00013193,
0.00044796,
0.00083448,
0.00540572,
0.00352615,
0.00026154,
0.00125032,
0.00031326,
0.00000000,
],
[
0.00000000,
0.00013193,
0.00044796,
0.00083448,
0.00540572,
0.00352615,
0.00026154,
0.00125032,
0.00031326,
0.00000000,
],
],
[
[
-0.00000000,
0.00320096,
0.00007269,
0.00008693,
0.00000484,
0.00000526,
0.00000203,
0.00000189,
0.00000039,
0.00000000,
],
[
0.00000000,
0.00055197,
0.00022117,
0.00012818,
0.00000952,
0.00000843,
0.00000335,
0.00000334,
0.00000105,
0.00000000,
],
[
0.00000000,
0.00014311,
0.00162453,
0.00335197,
0.00009023,
0.00005696,
0.00004058,
0.00002353,
0.00001610,
0.00000000,
],
[
0.00000000,
0.00018525,
0.00041136,
0.00272557,
0.01393003,
0.00057594,
0.00031515,
0.00034235,
0.00017612,
0.00000000,
],
[
0.00000000,
0.00016156,
0.00034722,
0.00132145,
0.04066875,
0.00134953,
0.00030062,
0.00087157,
0.00029123,
0.00000000,
],
[
-0.00000000,
0.00016722,
0.00030198,
0.00114459,
0.01745510,
0.00141671,
0.00035262,
0.00084328,
0.00029789,
0.00000000,
],
],
]
else:
ref_shifts = [
[
[
-0.00496121,
-0.00517009,
-0.00553453,
-0.00380260,
-0.00343693,
-0.00177772,
0.00067068,
0.00210892,
0.00377279,
0.00262923,
],
[
-0.00496121,
-0.00517009,
-0.00553453,
-0.00380260,
-0.00343693,
-0.00177772,
0.00067068,
0.00210892,
0.00377279,
0.00262923,
],
[
-0.01202453,
-0.01275071,
-0.01219115,
-0.00903023,
-0.01073474,
-0.01031603,
-0.00662244,
0.00195097,
0.01211624,
0.01038809,
],
[
-0.12223984,
-0.12186467,
-0.10941590,
-0.09614195,
-0.10716126,
-0.13519843,
-0.20168604,
-0.07470229,
0.15969340,
0.16315933,
],
[
-0.12114484,
-0.11964957,
-0.10778748,
-0.09630527,
-0.10697739,
-0.12776042,
-0.20180332,
-0.08373516,
0.15847215,
0.16377770,
],
[
-0.12114484,
-0.11964957,
-0.10778748,
-0.09630527,
-0.10697739,
-0.12776042,
-0.20180332,
-0.08373516,
0.15847215,
0.16377770,
],
],
[
[
-0.01893616,
-0.01942744,
-0.02339230,
-0.01392117,
-0.00568286,
-0.00676300,
-0.00298881,
0.00566943,
0.01407270,
0.01103835,
],
[
-0.01882494,
-0.01903756,
-0.02188396,
-0.01496531,
-0.00891242,
-0.00912826,
-0.00600410,
0.00298462,
0.01486886,
0.01365297,
],
[
-0.04155678,
-0.04201981,
-0.04661298,
-0.03446840,
-0.02577765,
-0.03332493,
-0.02460421,
-0.00520459,
0.03117184,
0.04283480,
],
[
-0.09551912,
-0.09141204,
-0.08650838,
-0.07531933,
-0.07736040,
-0.10097208,
-0.11850788,
-0.05857319,
0.08971321,
0.12793090,
],
[
-0.11821481,
-0.11425389,
-0.10865996,
-0.10189830,
-0.10716084,
-0.13231357,
-0.19690540,
-0.09252776,
0.14571718,
0.17189918,
],
[
-0.11276994,
-0.10757084,
-0.10142181,
-0.09519851,
-0.10205844,
-0.12882962,
-0.18549798,
-0.08931099,
0.13476362,
0.16481222,
],
],
]
ref_spec_funcs = [
[
[
0.00000000,
0.00001654,
0.00000271,
0.00000163,
0.00000013,
0.00000032,
0.00000006,
0.00000006,
0.00000001,
0.00000000,
],
[
0.00000000,
0.00001654,
0.00000271,
0.00000163,
0.00000013,
0.00000032,
0.00000006,
0.00000006,
0.00000001,
0.00000000,
],
[
0.00000000,
0.00223584,
0.00005772,
0.00001181,
0.00000067,
0.00000258,
0.00000098,
0.00000112,
0.00000029,
0.00000000,
],
[
0.00000000,
0.00012499,
0.00049591,
0.00088903,
0.00629350,
0.00256216,
0.00044658,
0.00118784,
0.00029727,
0.00000000,
],
[
0.00000000,
0.00013178,
0.00044892,
0.00082915,
0.00537908,
0.00352736,
0.00026142,
0.00125148,
0.00031313,
0.00000000,
],
[
0.00000000,
0.00013178,
0.00044892,
0.00082915,
0.00537908,
0.00352736,
0.00026142,
0.00125148,
0.00031313,
0.00000000,
],
],
[
[
-0.00000000,
0.00320916,
0.00007278,
0.00008659,
0.00000484,
0.00000528,
0.00000202,
0.00000189,
0.00000039,
0.00000000,
],
[
0.00000000,
0.00055185,
0.00022088,
0.00012806,
0.00000942,
0.00000847,
0.00000336,
0.00000334,
0.00000105,
0.00000000,
],
[
0.00000000,
0.00014260,
0.00161534,
0.00335775,
0.00008970,
0.00005681,
0.00004062,
0.00002350,
0.00001611,
0.00000000,
],
[
0.00000000,
0.00018553,
0.00041303,
0.00273911,
0.01389025,
0.00057626,
0.00031509,
0.00034209,
0.00017625,
0.00000000,
],
[
0.00000000,
0.00016153,
0.00034741,
0.00131818,
0.04062696,
0.00134833,
0.00030074,
0.00087080,
0.00029123,
0.00000000,
],
[
-0.00000000,
0.00016730,
0.00030334,
0.00113844,
0.01748363,
0.00141678,
0.00035292,
0.00084285,
0.00029790,
0.00000000,
],
],
]
# for line in np.swapaxes(sf.spectral_functions, -2, -1).reshape(-1, 6):
# print(("%.7f, " * 6) % tuple(line))
# raise
np.testing.assert_allclose(
shifts, np.swapaxes(sf.shifts, -2, -1).ravel(), atol=1e-2
)
np.testing.assert_allclose(
spec_funcs,
np.swapaxes(sf.spectral_functions, -2, -1).ravel(),
atol=1e-2,
rtol=1e-2,
)
def test_SpectralFunction_band_indices(si_pbesol: Phono3py):
"""Spectral function of Si."""
si_pbesol.mesh_numbers = [9, 9, 9]
if check_bands:
si_pbesol.band_indices = [[4, 5]]
si_pbesol.init_phph_interaction()
sf = SpectralFunction(
@ -297,18 +627,35 @@ def test_SpectralFunction_band_indices(si_pbesol: Phono3py):
)
sf.run()
# for line in np.swapaxes(sf.spectral_functions, -2, -1).reshape(-1, 6):
# print(("%.7f, " * 6) % tuple(line))
# raise
if check_bands:
np.testing.assert_allclose(
np.reshape(shifts, (-1, 6))[:, [4, 5]],
np.swapaxes(sf.shifts, -2, -1).reshape(-1, 2),
np.array(ref_shifts)[:, [4, 5], :],
sf.shifts[0, 0],
atol=1e-2,
)
np.testing.assert_allclose(
np.reshape(spec_funcs, (-1, 6))[:, [4, 5]],
np.swapaxes(sf.spectral_functions, -2, -1).reshape(-1, 2),
np.array(ref_spec_funcs)[:, [4, 5], :],
sf.spectral_functions[0, 0],
atol=1e-2,
rtol=1e-2,
)
else:
# for line in sf.shifts[0, 0, 0]:
# print("[", ",".join([f"{val:.8f}" for val in line]), "],")
# print("")
# for line in sf.shifts[0, 0, 1]:
# print("[", ",".join([f"{val:.8f}" for val in line]), "],")
for line in sf.spectral_functions[0, 0, 0]:
print("[", ",".join([f"{val:.8f}" for val in line]), "],")
print("")
for line in sf.spectral_functions[0, 0, 1]:
print("[", ",".join([f"{val:.8f}" for val in line]), "],")
np.testing.assert_allclose(ref_shifts, sf.shifts[0, 0], atol=1e-2)
np.testing.assert_allclose(
ref_spec_funcs, sf.spectral_functions[0, 0], atol=1e-2, rtol=1e-2
)
if check_bands:
si_pbesol.band_indices = None