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

Merge branch 'develop' into publish-gh-pages

This commit is contained in:
Atsushi Togo 2023-12-25 22:43:56 +09:00
parent 259f3e87b5 8a0d9ebf0e
commit 8287ddb6f8
42 changed files with 4268 additions and 2934 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

6
.pre-commit-config.yaml
View File

@ -19,7 +19,7 @@ repos:
- "--ignore=E203,W503"
- repo: https://github.com/psf/black
rev: 23.10.1
rev: 23.12.0
hooks:
- id: black
args:
@ -31,13 +31,13 @@ repos:
- id: pydocstyle
- repo: https://github.com/pycqa/isort
rev: 5.12.0
rev: 5.13.2
hooks:
- id: isort
name: isort (python)
- repo: https://github.com/Takishima/cmake-pre-commit-hooks
rev: v1.9.4
rev: v1.9.5
hooks:
- id: clang-format
args:

28
c/_phono3py.c
View File

@ -275,6 +275,7 @@ static PyObject *py_get_interaction(PyObject *self, PyObject *args) {
PyArrayObject *py_p2s_map;
PyArrayObject *py_s2p_map;
PyArrayObject *py_band_indices;
PyArrayObject *py_all_shortest;
double cutoff_frequency;
long symmetrize_fc3_q;
long make_r0_average;
@ -293,6 +294,7 @@ static PyObject *py_get_interaction(PyObject *self, PyObject *args) {
double(*svecs)[3];
long(*multi)[2];
double *masses;
char *all_shortest;
long *p2s;
long *s2p;
long *band_indices;
@ -300,12 +302,12 @@ static PyObject *py_get_interaction(PyObject *self, PyObject *args) {
long i;
long is_compact_fc3;
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOlldl", &py_fc3_normal_squared,
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOllOdl", &py_fc3_normal_squared,
&py_g_zero, &py_frequencies, &py_eigenvectors,
&py_triplets, &py_bz_grid_addresses, &py_D_diag,
&py_Q, &py_fc3, &py_svecs, &py_multi, &py_masses,
&py_p2s_map, &py_s2p_map, &py_band_indices,
&symmetrize_fc3_q, &make_r0_average,
&symmetrize_fc3_q, &make_r0_average, &py_all_shortest,
&cutoff_frequency, &openmp_per_triplets)) {
return NULL;
}
@ -336,12 +338,13 @@ static PyObject *py_get_interaction(PyObject *self, PyObject *args) {
p2s = (long *)PyArray_DATA(py_p2s_map);
s2p = (long *)PyArray_DATA(py_s2p_map);
band_indices = (long *)PyArray_DATA(py_band_indices);
all_shortest = (char *)PyArray_DATA(py_all_shortest);
ph3py_get_interaction(fc3_normal_squared, g_zero, freqs, eigvecs, triplets,
num_triplets, bz_grid_addresses, D_diag, Q, fc3,
is_compact_fc3, svecs, multi_dims, multi, masses, p2s,
s2p, band_indices, symmetrize_fc3_q, make_r0_average,
cutoff_frequency, openmp_per_triplets);
all_shortest, cutoff_frequency, openmp_per_triplets);
free(fc3_normal_squared);
fc3_normal_squared = NULL;
@ -370,6 +373,7 @@ static PyObject *py_get_pp_collision(PyObject *self, PyObject *args) {
PyArrayObject *py_s2p_map;
PyArrayObject *py_band_indices;
PyArrayObject *py_temperatures;
PyArrayObject *py_all_shortest;
double cutoff_frequency;
long is_NU;
long symmetrize_fc3_q;
@ -396,18 +400,19 @@ static PyObject *py_get_pp_collision(PyObject *self, PyObject *args) {
long *s2p;
Larray *band_indices;
Darray *temperatures;
char *all_shortest;
long multi_dims[2];
long i;
long is_compact_fc3;
if (!PyArg_ParseTuple(
args, "OOOOOOOOlOOOOOOOOOOllldl", &py_gamma,
args, "OOOOOOOOlOOOOOOOOOOlllOdl", &py_gamma,
&py_relative_grid_address, &py_frequencies, &py_eigenvectors,
&py_triplets, &py_triplet_weights, &py_bz_grid_addresses,
&py_bz_map, &bz_grid_type, &py_D_diag, &py_Q, &py_fc3, &py_svecs,
&py_multi, &py_masses, &py_p2s_map, &py_s2p_map, &py_band_indices,
&py_temperatures, &is_NU, &symmetrize_fc3_q, &make_r0_average,
&cutoff_frequency, &openmp_per_triplets)) {
&py_all_shortest, &cutoff_frequency, &openmp_per_triplets)) {
return NULL;
}
@ -439,13 +444,14 @@ static PyObject *py_get_pp_collision(PyObject *self, PyObject *args) {
s2p = (long *)PyArray_DATA(py_s2p_map);
band_indices = convert_to_larray(py_band_indices);
temperatures = convert_to_darray(py_temperatures);
all_shortest = (char *)PyArray_DATA(py_all_shortest);
ph3py_get_pp_collision(
gamma, relative_grid_address, frequencies, eigenvectors, triplets,
num_triplets, triplet_weights, bz_grid_addresses, bz_map, bz_grid_type,
D_diag, Q, fc3, is_compact_fc3, svecs, multi_dims, multi, masses, p2s,
s2p, band_indices, temperatures, is_NU, symmetrize_fc3_q,
make_r0_average, cutoff_frequency, openmp_per_triplets);
make_r0_average, all_shortest, cutoff_frequency, openmp_per_triplets);
free(band_indices);
band_indices = NULL;
@ -473,6 +479,7 @@ static PyObject *py_get_pp_collision_with_sigma(PyObject *self,
PyArrayObject *py_s2p_map;
PyArrayObject *py_band_indices;
PyArrayObject *py_temperatures;
PyArrayObject *py_all_shortest;
long is_NU;
long symmetrize_fc3_q;
double sigma;
@ -498,17 +505,19 @@ static PyObject *py_get_pp_collision_with_sigma(PyObject *self,
long *s2p;
Larray *band_indices;
Darray *temperatures;
char *all_shortest;
long multi_dims[2];
long i;
long is_compact_fc3;
if (!PyArg_ParseTuple(
args, "OddOOOOOOOOOOOOOOOllldl", &py_gamma, &sigma, &sigma_cutoff,
args, "OddOOOOOOOOOOOOOOOlllOdl", &py_gamma, &sigma, &sigma_cutoff,
&py_frequencies, &py_eigenvectors, &py_triplets,
&py_triplet_weights, &py_bz_grid_addresses, &py_D_diag, &py_Q,
&py_fc3, &py_svecs, &py_multi, &py_masses, &py_p2s_map, &py_s2p_map,
&py_band_indices, &py_temperatures, &is_NU, &symmetrize_fc3_q,
&make_r0_average, &cutoff_frequency, &openmp_per_triplets)) {
&make_r0_average, &py_all_shortest, &cutoff_frequency,
&openmp_per_triplets)) {
return NULL;
}
@ -537,13 +546,14 @@ static PyObject *py_get_pp_collision_with_sigma(PyObject *self,
s2p = (long *)PyArray_DATA(py_s2p_map);
band_indices = convert_to_larray(py_band_indices);
temperatures = convert_to_darray(py_temperatures);
all_shortest = (char *)PyArray_DATA(py_all_shortest);
ph3py_get_pp_collision_with_sigma(
gamma, sigma, sigma_cutoff, frequencies, eigenvectors, triplets,
num_triplets, triplet_weights, bz_grid_addresses, D_diag, Q, fc3,
is_compact_fc3, svecs, multi_dims, multi, masses, p2s, s2p,
band_indices, temperatures, is_NU, symmetrize_fc3_q, make_r0_average,
cutoff_frequency, openmp_per_triplets);
all_shortest, cutoff_frequency, openmp_per_triplets);
free(band_indices);
band_indices = NULL;

122
c/interaction.c
View File

@ -46,42 +46,37 @@
static const long index_exchange[6][3] = {{0, 1, 2}, {2, 0, 1}, {1, 2, 0},
{2, 1, 0}, {0, 2, 1}, {1, 0, 2}};
static void real_to_normal(
double *fc3_normal_squared, const long (*g_pos)[4], const long num_g_pos,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *eigvecs0,
const lapack_complex_double *eigvecs1,
const lapack_complex_double *eigvecs2, const double *fc3,
const long is_compact_fc3, const double q_vecs[3][3], /* q0, q1, q2 */
const double (*svecs)[3], const long multi_dims[2],
const long (*multiplicity)[2], const double *masses, const long *p2s_map,
const long *s2p_map, const long *band_indices, const long num_band,
const double cutoff_frequency, const long triplet_index,
const long num_triplets, const long make_r0_average,
const long openmp_at_bands);
static void real_to_normal(double *fc3_normal_squared, const long (*g_pos)[4],
const long num_g_pos, const double *freqs0,
const double *freqs1, const double *freqs2,
const lapack_complex_double *eigvecs0,
const lapack_complex_double *eigvecs1,
const lapack_complex_double *eigvecs2,
const double *fc3, const long is_compact_fc3,
const double q_vecs[3][3], /* q0, q1, q2 */
const AtomTriplets *atom_triplets,
const double *masses, const long *band_indices,
const long num_band, const double cutoff_frequency,
const long triplet_index, const long num_triplets,
const long openmp_at_bands);
static void real_to_normal_sym_q(
double *fc3_normal_squared, const long (*g_pos)[4], const long num_g_pos,
double *const freqs[3], lapack_complex_double *const eigvecs[3],
const double *fc3, const long is_compact_fc3,
const double q_vecs[3][3], /* q0, q1, q2 */
const double (*svecs)[3], const long multi_dims[2],
const long (*multiplicity)[2], const double *masses, const long *p2s_map,
const long *s2p_map, const long *band_indices, const long num_band0,
const long num_band, const double cutoff_frequency,
const long triplet_index, const long num_triplets,
const long make_r0_average, const long openmp_at_bands);
const AtomTriplets *atom_triplets, const double *masses,
const long *band_indices, const long num_band0, const long num_band,
const double cutoff_frequency, const long triplet_index,
const long num_triplets, const long openmp_at_bands);
/* fc3_normal_squared[num_triplets, num_band0, num_band, num_band] */
void itr_get_interaction(
Darray *fc3_normal_squared, const char *g_zero, const Darray *frequencies,
const lapack_complex_double *eigenvectors, const long (*triplets)[3],
const long num_triplets, const ConstBZGrid *bzgrid, const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency, const long make_r0_average,
const long openmp_per_triplets) {
const long is_compact_fc3, const AtomTriplets *atom_triplets,
const double *masses, const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency, const long openmp_per_triplets) {
long(*g_pos)[4];
long i;
long num_band, num_band0, num_band_prod, num_g_pos;
@ -94,7 +89,7 @@ void itr_get_interaction(
#ifdef _OPENMP
#pragma omp parallel for schedule(guided) private( \
num_g_pos, g_pos) if (openmp_per_triplets)
num_g_pos, g_pos) if (openmp_per_triplets)
#endif
for (i = 0; i < num_triplets; i++) {
g_pos = (long(*)[4])malloc(sizeof(long[4]) * num_band_prod);
@ -104,9 +99,8 @@ void itr_get_interaction(
itr_get_interaction_at_triplet(
fc3_normal_squared->data + i * num_band_prod, num_band0, num_band,
g_pos, num_g_pos, frequencies->data, eigenvectors, triplets[i],
bzgrid, fc3, is_compact_fc3, svecs, multi_dims, multiplicity,
masses, p2s_map, s2p_map, band_indices, symmetrize_fc3_q,
cutoff_frequency, i, num_triplets, make_r0_average,
bzgrid, fc3, is_compact_fc3, atom_triplets, masses, band_indices,
symmetrize_fc3_q, cutoff_frequency, i, num_triplets,
1 - openmp_per_triplets);
free(g_pos);
@ -119,13 +113,12 @@ void itr_get_interaction_at_triplet(
const long (*g_pos)[4], const long num_g_pos, const double *frequencies,
const lapack_complex_double *eigenvectors, const long triplet[3],
const ConstBZGrid *bzgrid, const double *fc3, const long is_compact_fc3,
const double (*svecs)[3], const long multi_dims[2],
const long (*multiplicity)[2], const double *masses, const long *p2s_map,
const long *s2p_map, const long *band_indices, const long symmetrize_fc3_q,
const AtomTriplets *atom_triplets, const double *masses,
const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency,
const long triplet_index, /* only for print */
const long num_triplets, /* only for print */
const long make_r0_average, const long openmp_at_bands) {
const long openmp_at_bands) {
long j, k;
double *freqs[3];
lapack_complex_double *eigvecs[3];
@ -155,9 +148,8 @@ void itr_get_interaction_at_triplet(
real_to_normal_sym_q(
fc3_normal_squared, g_pos, num_g_pos, freqs, eigvecs, fc3,
is_compact_fc3, q_vecs, /* q0, q1, q2 */
svecs, multi_dims, multiplicity, masses, p2s_map, s2p_map,
band_indices, num_band0, num_band, cutoff_frequency, triplet_index,
num_triplets, make_r0_average, openmp_at_bands);
atom_triplets, masses, band_indices, num_band0, num_band,
cutoff_frequency, triplet_index, num_triplets, openmp_at_bands);
for (j = 0; j < 3; j++) {
free(freqs[j]);
freqs[j] = NULL;
@ -173,26 +165,25 @@ void itr_get_interaction_at_triplet(
eigenvectors + triplet[1] * num_band * num_band,
eigenvectors + triplet[2] * num_band * num_band, fc3,
is_compact_fc3, q_vecs, /* q0, q1, q2 */
svecs, multi_dims, multiplicity, masses, p2s_map,
s2p_map, band_indices, num_band, cutoff_frequency,
triplet_index, num_triplets, make_r0_average,
atom_triplets, masses, band_indices, num_band,
cutoff_frequency, triplet_index, num_triplets,
openmp_at_bands);
}
}
static void real_to_normal(
double *fc3_normal_squared, const long (*g_pos)[4], const long num_g_pos,
const double *freqs0, const double *freqs1, const double *freqs2,
const lapack_complex_double *eigvecs0,
const lapack_complex_double *eigvecs1,
const lapack_complex_double *eigvecs2, const double *fc3,
const long is_compact_fc3, const double q_vecs[3][3], /* q0, q1, q2 */
const double (*svecs)[3], const long multi_dims[2],
const long (*multiplicity)[2], const double *masses, const long *p2s_map,
const long *s2p_map, const long *band_indices, const long num_band,
const double cutoff_frequency, const long triplet_index,
const long num_triplets, const long make_r0_average,
const long openmp_at_bands) {
static void real_to_normal(double *fc3_normal_squared, const long (*g_pos)[4],
const long num_g_pos, const double *freqs0,
const double *freqs1, const double *freqs2,
const lapack_complex_double *eigvecs0,
const lapack_complex_double *eigvecs1,
const lapack_complex_double *eigvecs2,
const double *fc3, const long is_compact_fc3,
const double q_vecs[3][3], /* q0, q1, q2 */
const AtomTriplets *atom_triplets,
const double *masses, const long *band_indices,
const long num_band, const double cutoff_frequency,
const long triplet_index, const long num_triplets,
const long openmp_at_bands) {
lapack_complex_double *fc3_reciprocal;
lapack_complex_double comp_zero;
long i;
@ -203,16 +194,8 @@ static void real_to_normal(
for (i = 0; i < num_band * num_band * num_band; i++) {
fc3_reciprocal[i] = comp_zero;
}
r2r_real_to_reciprocal(fc3_reciprocal, q_vecs, fc3, is_compact_fc3, svecs,
multi_dims, multiplicity, p2s_map, s2p_map,
make_r0_average, openmp_at_bands);
if (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);
}
}
r2r_real_to_reciprocal(fc3_reciprocal, q_vecs, fc3, is_compact_fc3,
atom_triplets, openmp_at_bands);
#ifdef MEASURE_R2N
if (openmp_at_bands && num_triplets > 0) {
@ -234,12 +217,10 @@ static void real_to_normal_sym_q(
double *const freqs[3], lapack_complex_double *const eigvecs[3],
const double *fc3, const long is_compact_fc3,
const double q_vecs[3][3], /* q0, q1, q2 */
const double (*svecs)[3], const long multi_dims[2],
const long (*multiplicity)[2], const double *masses, const long *p2s_map,
const long *s2p_map, const long *band_indices, const long num_band0,
const long num_band, const double cutoff_frequency,
const long triplet_index, const long num_triplets,
const long make_r0_average, const long openmp_at_bands) {
const AtomTriplets *atom_triplets, const double *masses,
const long *band_indices, const long num_band0, const long num_band,
const double cutoff_frequency, const long triplet_index,
const long num_triplets, const long openmp_at_bands) {
long i, j, k, l;
long band_ex[3];
double q_vecs_ex[3][3];
@ -264,9 +245,8 @@ static void real_to_normal_sym_q(
freqs[index_exchange[i][2]], eigvecs[index_exchange[i][0]],
eigvecs[index_exchange[i][1]], eigvecs[index_exchange[i][2]], fc3,
is_compact_fc3, q_vecs_ex, /* q0, q1, q2 */
svecs, multi_dims, multiplicity, masses, p2s_map, s2p_map,
band_indices, num_band, cutoff_frequency, triplet_index,
num_triplets, make_r0_average, openmp_at_bands);
atom_triplets, masses, band_indices, num_band, cutoff_frequency,
triplet_index, num_triplets, openmp_at_bands);
for (j = 0; j < num_band0; j++) {
for (k = 0; k < num_band; k++) {
for (l = 0; l < num_band; l++) {

17
c/interaction.h
View File

@ -38,28 +38,25 @@
#include "bzgrid.h"
#include "lapack_wrapper.h"
#include "phonoc_array.h"
#include "real_to_reciprocal.h"
void itr_get_interaction(
Darray *fc3_normal_squared, const char *g_zero, const Darray *frequencies,
const lapack_complex_double *eigenvectors, const long (*triplets)[3],
const long num_triplets, const ConstBZGrid *bzgrid, const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency, const long make_r0_average,
const long openmp_per_triplets);
const long is_compact_fc3, const AtomTriplets *atom_triplets,
const double *masses, const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency, const long openmp_per_triplets);
void itr_get_interaction_at_triplet(
double *fc3_normal_squared, const long num_band0, const long num_band,
const long (*g_pos)[4], const long num_g_pos, const double *frequencies,
const lapack_complex_double *eigenvectors, const long triplet[3],
const ConstBZGrid *bzgrid, const double *fc3, const long is_compact_fc3,
const double (*svecs)[3], const long multi_dims[2],
const long (*multiplicity)[2], const double *masses, const long *p2s_map,
const long *s2p_map, const long *band_indices, const long symmetrize_fc3_q,
const AtomTriplets *atom_triplets, const double *masses,
const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency,
const long triplet_index, /* only for print */
const long num_triplets, /* only for print */
const long make_r0_average, const long openmp_at_bands);
const long openmp_at_bands);
#endif

94
c/phono3py.c
View File

@ -49,6 +49,7 @@
#include "phonoc_array.h"
#include "pp_collision.h"
#include "real_self_energy.h"
#include "real_to_reciprocal.h"
#include "tetrahedron_method.h"
#include "triplet.h"
#include "triplet_iw.h"
@ -66,9 +67,10 @@ long ph3py_get_interaction(
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const long *band_indices, const long symmetrize_fc3_q,
const long make_r0_average, const double cutoff_frequency,
const long openmp_per_triplets) {
const long make_r0_average, const char *all_shortest,
const double cutoff_frequency, const long openmp_per_triplets) {
ConstBZGrid *bzgrid;
AtomTriplets *atom_triplets;
long i, j;
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) {
@ -85,12 +87,30 @@ long ph3py_get_interaction(
}
}
if ((atom_triplets = (AtomTriplets *)malloc(sizeof(AtomTriplets))) ==
NULL) {
warning_print("Memory could not be allocated.");
return 0;
}
atom_triplets->svecs = svecs;
atom_triplets->multi_dims[0] = multi_dims[0];
atom_triplets->multi_dims[1] = multi_dims[1];
atom_triplets->multiplicity = multiplicity;
atom_triplets->p2s_map = p2s_map;
atom_triplets->s2p_map = s2p_map;
atom_triplets->make_r0_average = make_r0_average;
atom_triplets->all_shortest = all_shortest;
itr_get_interaction(fc3_normal_squared, g_zero, frequencies,
(lapack_complex_double *)eigenvectors, triplets,
num_triplets, bzgrid, fc3, is_compact_fc3, svecs,
multi_dims, multiplicity, masses, p2s_map, s2p_map,
band_indices, symmetrize_fc3_q, cutoff_frequency,
make_r0_average, openmp_per_triplets);
num_triplets, bzgrid, fc3, is_compact_fc3,
atom_triplets, masses, band_indices, symmetrize_fc3_q,
cutoff_frequency, openmp_per_triplets);
free(atom_triplets);
atom_triplets = NULL;
free(bzgrid);
bzgrid = NULL;
@ -109,9 +129,11 @@ long ph3py_get_pp_collision(
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const Larray *band_indices, const Darray *temperatures, const long is_NU,
const long symmetrize_fc3_q, const double make_r0_average,
const double cutoff_frequency, const long openmp_per_triplets) {
const long symmetrize_fc3_q, const long make_r0_average,
const char *all_shortest, const double cutoff_frequency,
const long openmp_per_triplets) {
ConstBZGrid *bzgrid;
AtomTriplets *atom_triplets;
long i, j;
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) {
@ -130,13 +152,30 @@ long ph3py_get_pp_collision(
}
}
if ((atom_triplets = (AtomTriplets *)malloc(sizeof(AtomTriplets))) ==
NULL) {
warning_print("Memory could not be allocated.");
return 0;
}
atom_triplets->svecs = svecs;
atom_triplets->multi_dims[0] = multi_dims[0];
atom_triplets->multi_dims[1] = multi_dims[1];
atom_triplets->multiplicity = multiplicity;
atom_triplets->p2s_map = p2s_map;
atom_triplets->s2p_map = s2p_map;
atom_triplets->make_r0_average = make_r0_average;
atom_triplets->all_shortest = all_shortest;
ppc_get_pp_collision(imag_self_energy, relative_grid_address, frequencies,
(lapack_complex_double *)eigenvectors, triplets,
num_triplets, triplet_weights, bzgrid, fc3,
is_compact_fc3, svecs, multi_dims, multiplicity,
masses, p2s_map, s2p_map, band_indices, temperatures,
is_NU, symmetrize_fc3_q, cutoff_frequency,
make_r0_average, openmp_per_triplets);
is_compact_fc3, atom_triplets, masses, band_indices,
temperatures, is_NU, symmetrize_fc3_q,
cutoff_frequency, openmp_per_triplets);
free(atom_triplets);
atom_triplets = NULL;
free(bzgrid);
bzgrid = NULL;
@ -154,9 +193,11 @@ long ph3py_get_pp_collision_with_sigma(
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const Larray *band_indices, const Darray *temperatures, const long is_NU,
const long symmetrize_fc3_q, const double make_r0_average,
const double cutoff_frequency, const long openmp_per_triplets) {
const long symmetrize_fc3_q, const long make_r0_average,
const char *all_shortest, const double cutoff_frequency,
const long openmp_per_triplets) {
ConstBZGrid *bzgrid;
AtomTriplets *atom_triplets;
long i, j;
if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) {
@ -173,14 +214,31 @@ long ph3py_get_pp_collision_with_sigma(
}
}
if ((atom_triplets = (AtomTriplets *)malloc(sizeof(AtomTriplets))) ==
NULL) {
warning_print("Memory could not be allocated.");
return 0;
}
atom_triplets->svecs = svecs;
atom_triplets->multi_dims[0] = multi_dims[0];
atom_triplets->multi_dims[1] = multi_dims[1];
atom_triplets->multiplicity = multiplicity;
atom_triplets->p2s_map = p2s_map;
atom_triplets->s2p_map = s2p_map;
atom_triplets->make_r0_average = make_r0_average;
atom_triplets->all_shortest = all_shortest;
ppc_get_pp_collision_with_sigma(
imag_self_energy, sigma, sigma_cutoff, frequencies,
(lapack_complex_double *)eigenvectors, triplets, num_triplets,
triplet_weights, bzgrid, fc3, is_compact_fc3, svecs, multi_dims,
multiplicity, masses, p2s_map, s2p_map, band_indices, temperatures,
is_NU, symmetrize_fc3_q, cutoff_frequency, make_r0_average,
triplet_weights, bzgrid, fc3, is_compact_fc3, atom_triplets, masses,
band_indices, temperatures, is_NU, symmetrize_fc3_q, cutoff_frequency,
openmp_per_triplets);
free(atom_triplets);
atom_triplets = NULL;
free(bzgrid);
bzgrid = NULL;
@ -600,7 +658,7 @@ void ph3py_expand_collision_matrix(double *collision_matrix,
j * num_column * num_column);
#ifdef _OPENMP
#pragma omp parallel for private(ir_gp, adrs_shift_plus, colmat_copy, l, gp_r, \
m, n, p)
m, n, p)
#endif
for (k = 0; k < num_ir_gp; k++) {
ir_gp = ir_grid_points[k];

11
c/phono3py.h
View File

@ -51,7 +51,8 @@ long ph3py_get_interaction(
const long multi_dims[2], const long (*multi)[2], const double *masses,
const long *p2s_map, const long *s2p_map, const long *band_indices,
const long symmetrize_fc3_q, const long make_r0_average,
const double cutoff_frequency, const long openmp_per_triplets);
const char *all_shortest, const double cutoff_frequency,
const long openmp_per_triplets);
long ph3py_get_pp_collision(
double *imag_self_energy,
const long relative_grid_address[24][4][3], /* thm */
@ -64,8 +65,8 @@ long ph3py_get_pp_collision(
const long multi_dims[2], const long (*multi)[2], const double *masses,
const long *p2s_map, const long *s2p_map, const Larray *band_indices,
const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q,
const double make_r0_average, const double cutoff_frequency,
const long openmp_per_triplets);
const long make_r0_average, const char *all_shortest,
const double cutoff_frequency, const long openmp_per_triplets);
long ph3py_get_pp_collision_with_sigma(
double *imag_self_energy, const double sigma, const double sigma_cutoff,
const double *frequencies, const _lapack_complex_double *eigenvectors,
@ -76,8 +77,8 @@ long ph3py_get_pp_collision_with_sigma(
const long multi_dims[2], const long (*multi)[2], const double *masses,
const long *p2s_map, const long *s2p_map, const Larray *band_indices,
const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q,
const double make_r0_average, const double cutoff_frequency,
const long openmp_per_triplets);
const long make_r0_average, const char *all_shortest,
const double cutoff_frequency, const long openmp_per_triplets);
void ph3py_get_imag_self_energy_at_bands_with_g(
double *imag_self_energy, const Darray *fc3_normal_squared,
const double *frequencies, const long (*triplets)[3],

63
c/pp_collision.c
View File

@ -42,6 +42,7 @@
#include "lapack_wrapper.h"
#include "phonoc_array.h"
#include "phonoc_utils.h"
#include "real_to_reciprocal.h"
#include "triplet.h"
#include "triplet_iw.h"
@ -51,12 +52,9 @@ static void get_collision(
const char *g_zero, const double *frequencies,
const lapack_complex_double *eigenvectors, const long triplet[3],
const long triplet_weight, const ConstBZGrid *bzgrid, const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency, const long make_r0_average,
const long openmp_per_triplets);
const long is_compact_fc3, const AtomTriplets *atom_triplets,
const double *masses, const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency, const long openmp_per_triplets);
static void finalize_ise(double *imag_self_energy, const double *ise,
const long (*bz_grid_address)[3],
const long (*triplets)[3], const long num_triplets,
@ -69,12 +67,10 @@ void ppc_get_pp_collision(
const double *frequencies, const lapack_complex_double *eigenvectors,
const long (*triplets)[3], const long num_triplets,
const long *triplet_weights, const ConstBZGrid *bzgrid, const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const Larray *band_indices, const Darray *temperatures, const long is_NU,
const long symmetrize_fc3_q, const double cutoff_frequency,
const long make_r0_average, const long openmp_per_triplets) {
const long is_compact_fc3, const AtomTriplets *atom_triplets,
const double *masses, const Larray *band_indices,
const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q,
const double cutoff_frequency, const long openmp_per_triplets) {
long i;
long num_band, num_band0, num_band_prod, num_temps;
double *ise, *freqs_at_gp, *g;
@ -87,7 +83,7 @@ void ppc_get_pp_collision(
g_zero = NULL;
num_band0 = band_indices->dims[0];
num_band = multi_dims[1] * 3;
num_band = atom_triplets->multi_dims[1] * 3;
num_band_prod = num_band0 * num_band * num_band;
num_temps = temperatures->dims[0];
ise =
@ -103,7 +99,7 @@ void ppc_get_pp_collision(
#ifdef _OPENMP
#pragma omp parallel for schedule(guided) private( \
g, g_zero) if (openmp_per_triplets)
g, g_zero) if (openmp_per_triplets)
#endif
for (i = 0; i < num_triplets; i++) {
g = (double *)malloc(sizeof(double) * 2 * num_band_prod);
@ -118,9 +114,8 @@ void ppc_get_pp_collision(
get_collision(ise + i * num_temps * num_band0, num_band0, num_band,
num_temps, temperatures->data, g, g_zero, frequencies,
eigenvectors, triplets[i], triplet_weights[i], bzgrid,
fc3, is_compact_fc3, svecs, multi_dims, multiplicity,
masses, p2s_map, s2p_map, band_indices->data,
symmetrize_fc3_q, cutoff_frequency, make_r0_average,
fc3, is_compact_fc3, atom_triplets, masses,
band_indices->data, symmetrize_fc3_q, cutoff_frequency,
openmp_per_triplets);
free(g_zero);
@ -143,12 +138,10 @@ void ppc_get_pp_collision_with_sigma(
const double *frequencies, const lapack_complex_double *eigenvectors,
const long (*triplets)[3], const long num_triplets,
const long *triplet_weights, const ConstBZGrid *bzgrid, const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const Larray *band_indices, const Darray *temperatures, const long is_NU,
const long symmetrize_fc3_q, const double cutoff_frequency,
const long make_r0_average, const long openmp_per_triplets) {
const long is_compact_fc3, const AtomTriplets *atom_triplets,
const double *masses, const Larray *band_indices,
const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q,
const double cutoff_frequency, const long openmp_per_triplets) {
long i;
long num_band, num_band0, num_band_prod, num_temps;
long const_adrs_shift;
@ -162,7 +155,7 @@ void ppc_get_pp_collision_with_sigma(
g_zero = NULL;
num_band0 = band_indices->dims[0];
num_band = multi_dims[1] * 3;
num_band = atom_triplets->multi_dims[1] * 3;
num_band_prod = num_band0 * num_band * num_band;
num_temps = temperatures->dims[0];
const_adrs_shift = num_band_prod;
@ -179,7 +172,7 @@ void ppc_get_pp_collision_with_sigma(
#ifdef _OPENMP
#pragma omp parallel for schedule(guided) private( \
g, g_zero) if (openmp_per_triplets)
g, g_zero) if (openmp_per_triplets)
#endif
for (i = 0; i < num_triplets; i++) {
g = (double *)malloc(sizeof(double) * 2 * num_band_prod);
@ -191,9 +184,8 @@ void ppc_get_pp_collision_with_sigma(
get_collision(ise + i * num_temps * num_band0, num_band0, num_band,
num_temps, temperatures->data, g, g_zero, frequencies,
eigenvectors, triplets[i], triplet_weights[i], bzgrid,
fc3, is_compact_fc3, svecs, multi_dims, multiplicity,
masses, p2s_map, s2p_map, band_indices->data,
symmetrize_fc3_q, cutoff_frequency, make_r0_average,
fc3, is_compact_fc3, atom_triplets, masses,
band_indices->data, symmetrize_fc3_q, cutoff_frequency,
openmp_per_triplets);
free(g_zero);
@ -217,12 +209,9 @@ static void get_collision(
const char *g_zero, const double *frequencies,
const lapack_complex_double *eigenvectors, const long triplet[3],
const long triplet_weight, const ConstBZGrid *bzgrid, const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency, const long make_r0_average,
const long openmp_per_triplets) {
const long is_compact_fc3, const AtomTriplets *atom_triplets,
const double *masses, const long *band_indices, const long symmetrize_fc3_q,
const double cutoff_frequency, const long openmp_per_triplets) {
long i;
long num_band_prod, num_g_pos;
double *fc3_normal_squared;
@ -243,9 +232,9 @@ static void get_collision(
itr_get_interaction_at_triplet(
fc3_normal_squared, num_band0, num_band, g_pos, num_g_pos, frequencies,
eigenvectors, triplet, bzgrid, fc3, is_compact_fc3, svecs, multi_dims,
multiplicity, masses, p2s_map, s2p_map, band_indices, symmetrize_fc3_q,
cutoff_frequency, 0, 0, make_r0_average, 1 - openmp_per_triplets);
eigenvectors, triplet, bzgrid, fc3, is_compact_fc3, atom_triplets,
masses, band_indices, symmetrize_fc3_q, cutoff_frequency, 0, 0,
1 - openmp_per_triplets);
ise_imag_self_energy_at_triplet(
ise, num_band0, num_band, fc3_normal_squared, frequencies, triplet,

21
c/pp_collision.h
View File

@ -38,6 +38,7 @@
#include "bzgrid.h"
#include "lapack_wrapper.h"
#include "phonoc_array.h"
#include "real_to_reciprocal.h"
void ppc_get_pp_collision(
double *imag_self_energy,
@ -45,22 +46,18 @@ void ppc_get_pp_collision(
const double *frequencies, const lapack_complex_double *eigenvectors,
const long (*triplets)[3], const long num_triplets,
const long *triplet_weights, const ConstBZGrid *bzgrid, const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const Larray *band_indices, const Darray *temperatures, const long is_NU,
const long symmetrize_fc3_q, const double cutoff_frequency,
const long make_r0_average, const long openmp_per_triplets);
const long is_compact_fc3, const AtomTriplets *atom_triplets,
const double *masses, const Larray *band_indices,
const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q,
const double cutoff_frequency, const long openmp_per_triplets);
void ppc_get_pp_collision_with_sigma(
double *imag_self_energy, const double sigma, const double sigma_cutoff,
const double *frequencies, const lapack_complex_double *eigenvectors,
const long (*triplets)[3], const long num_triplets,
const long *triplet_weights, const ConstBZGrid *bzgrid, const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const double *masses, const long *p2s_map, const long *s2p_map,
const Larray *band_indices, const Darray *temperatures, const long is_NU,
const long symmetrize_fc3_q, const double cutoff_frequency,
const long make_r0_average, const long openmp_per_triplets);
const long is_compact_fc3, const AtomTriplets *atom_triplets,
const double *masses, const Larray *band_indices,
const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q,
const double cutoff_frequency, const long openmp_per_triplets);
#endif

378
c/real_to_reciprocal.c
View File

@ -43,161 +43,194 @@
#include "phonoc_array.h"
#include "phonoc_const.h"
static void real_to_reciprocal(
lapack_complex_double *fc3_reciprocal, const double q_vecs[3][3],
const double *fc3, const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const long *p2s_map, const long *s2p_map, const long make_r0_average,
const long openmp_at_bands);
static void real_to_reciprocal_elements(
lapack_complex_double *fc3_rec_elem, const double q1[3], const double q2[3],
const double *fc3, const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2], const long *p2s,
const long *s2p, const long pi0, const long pi1, const long pi2);
static void real_to_reciprocal_legacy(lapack_complex_double *fc3_reciprocal,
const double q_vecs[3][3],
const double *fc3,
const long is_compact_fc3,
const AtomTriplets *atom_triplets,
const long openmp_at_bands);
static void real_to_reciprocal_r0_average(lapack_complex_double *fc3_reciprocal,
const double q_vecs[3][3],
const double *fc3,
const long is_compact_fc3,
const AtomTriplets *atom_triplets,
const long openmp_at_bands);
static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
const double q1[3], const double q2[3],
const double *fc3,
const long is_compact_fc3,
const AtomTriplets *atom_triplets,
const long pi0, const long pi1,
const long pi2, const long leg_index);
static lapack_complex_double get_phase_factor(const double q[3],
const double (*svecs)[3],
const long multi[2]);
static lapack_complex_double get_pre_phase_factor(const long i_patom,
const double q_vecs[3][3],
const double (*svecs)[3],
const long multi_dims[2],
const long (*multiplicity)[2],
const long *p2s_map);
static lapack_complex_double get_pre_phase_factor(
const long i_patom, const double q_vecs[3][3],
const AtomTriplets *atom_triplets);
static lapack_complex_double sum_lapack_complex_double(lapack_complex_double a,
lapack_complex_double b);
/* fc3_reciprocal[num_patom, num_patom, num_patom, 3, 3, 3] */
void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
const double q_vecs[3][3], const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2],
const long (*multiplicity)[2], const long *p2s_map,
const long *s2p_map, const long make_r0_average,
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, svecs,
multi_dims, multiplicity, p2s_map, s2p_map,
make_r0_average, openmp_at_bands);
long i, num_band;
if (atom_triplets->make_r0_average) {
real_to_reciprocal_r0_average(fc3_reciprocal, q_vecs, fc3,
is_compact_fc3, atom_triplets,
openmp_at_bands);
num_band = atom_triplets->multi_dims[1] * 3;
for (i = 0; i < num_band * num_band * num_band; i++) {
fc3_reciprocal[i] = lapack_make_complex_double(
lapack_complex_double_real(fc3_reciprocal[i]) / 3,
lapack_complex_double_imag(fc3_reciprocal[i]) / 3);
}
} else {
real_to_reciprocal_legacy(fc3_reciprocal, q_vecs, fc3, is_compact_fc3,
atom_triplets, openmp_at_bands);
}
}
static void real_to_reciprocal(
lapack_complex_double *fc3_reciprocal, const double q_vecs[3][3],
const double *fc3, const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2],
const long *p2s_map, const long *s2p_map, const long make_r0_average,
const long openmp_at_bands) {
long i, j, k, l, m, n, jk, ik, ji;
static void real_to_reciprocal_legacy(lapack_complex_double *fc3_reciprocal,
const double q_vecs[3][3],
const double *fc3,
const long is_compact_fc3,
const AtomTriplets *atom_triplets,
const long openmp_at_bands) {
long i, j, k, l, m, n, ijk;
long num_patom, num_band;
lapack_complex_double pre_phase_factor, fc3_rec_elem[27], fc3_rec;
num_patom = multi_dims[1];
num_patom = atom_triplets->multi_dims[1];
num_band = num_patom * 3;
for (i = 0; i < num_patom; i++) {
pre_phase_factor = get_pre_phase_factor(i, q_vecs, svecs, multi_dims,
multiplicity, p2s_map);
#ifdef _OPENMP
#pragma omp parallel for private(j, k, l, m, n, fc3_rec_elem, \
fc3_rec) if (openmp_at_bands)
#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 (jk = 0; jk < num_patom * num_patom; jk++) {
j = jk / num_patom;
k = jk % num_patom;
real_to_reciprocal_elements(
fc3_rec_elem, q_vecs[1], q_vecs[2], fc3, is_compact_fc3, svecs,
multi_dims, multiplicity, p2s_map, s2p_map, i, j, k);
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);
}
}
}
}
}
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;
if (make_r0_average) {
for (j = 0; j < num_patom; j++) {
pre_phase_factor = get_pre_phase_factor(
j, q_vecs, svecs, multi_dims, multiplicity, p2s_map);
#ifdef _OPENMP
#pragma omp parallel for private(i, k, l, m, n, fc3_rec_elem, \
fc3_rec) if (openmp_at_bands)
#endif
for (ik = 0; ik < num_patom * num_patom; ik++) {
i = ik / num_patom;
k = ik % num_patom;
real_to_reciprocal_elements(
fc3_rec_elem, q_vecs[0], q_vecs[2], fc3, is_compact_fc3,
svecs, multi_dims, multiplicity, p2s_map, s2p_map, j, i, k);
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
fc3_rec = phonoc_complex_prod(
fc3_rec_elem[m * 9 + l * 3 + n],
pre_phase_factor);
fc3_reciprocal[(i * 3 + l) * num_band * num_band +
(j * 3 + m) * num_band + k * 3 + n] =
sum_lapack_complex_double(
fc3_reciprocal[(i * 3 + l) * num_band *
num_band +
(j * 3 + m) * num_band +
k * 3 + n],
fc3_rec);
}
}
}
}
}
pre_phase_factor = get_pre_phase_factor(i, q_vecs, atom_triplets);
for (k = 0; k < num_patom; k++) {
pre_phase_factor = get_pre_phase_factor(
k, q_vecs, svecs, multi_dims, multiplicity, p2s_map);
#ifdef _OPENMP
#pragma omp parallel for private(j, i, l, m, n, fc3_rec_elem, \
fc3_rec) if (openmp_at_bands)
#endif
for (ji = 0; ji < num_patom * num_patom; ji++) {
j = ji / num_patom;
i = ji % num_patom;
real_to_reciprocal_elements(
fc3_rec_elem, q_vecs[1], q_vecs[0], fc3, is_compact_fc3,
svecs, multi_dims, multiplicity, p2s_map, s2p_map, k, j, i);
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
fc3_rec = phonoc_complex_prod(
fc3_rec_elem[n * 9 + m * 3 + l],
pre_phase_factor);
real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[2], fc3,
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] =
sum_lapack_complex_double(
fc3_reciprocal[(i * 3 + l) * num_band *
num_band +
(j * 3 + m) * num_band +
k * 3 + n],
fc3_rec);
}
}
(j * 3 + m) * num_band + k * 3 + n],
fc3_rec);
}
}
}
}
}
static void real_to_reciprocal_elements(
lapack_complex_double *fc3_rec_elem, const double q1[3], const double q2[3],
const double *fc3, const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2], const long (*multiplicity)[2], const long *p2s,
const long *s2p, const long pi0, const long pi1, const long pi2) {
// 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++) {
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);
}
}
}
// fc3_rec is stored in a way swapping jm <-> il.
pre_phase_factor = get_pre_phase_factor(j, q_vecs, atom_triplets);
real_to_reciprocal_elements(fc3_rec_elem, q_vecs[0], q_vecs[2], fc3,
is_compact_fc3, atom_triplets, j, i, k, 2);
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
fc3_rec = phonoc_complex_prod(
fc3_rec_elem[m * 9 + l * 3 + n], pre_phase_factor);
fc3_reciprocal[(i * 3 + l) * num_band * num_band +
(j * 3 + m) * num_band + k * 3 + n] =
sum_lapack_complex_double(
fc3_reciprocal[(i * 3 + l) * num_band * num_band +
(j * 3 + m) * num_band + k * 3 + n],
fc3_rec);
}
}
}
// fc3_rec is stored in a way swapping kn <-> il.
pre_phase_factor = get_pre_phase_factor(k, q_vecs, atom_triplets);
real_to_reciprocal_elements(fc3_rec_elem, q_vecs[1], q_vecs[0], fc3,
is_compact_fc3, atom_triplets, k, j, i, 3);
for (l = 0; l < 3; l++) {
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
fc3_rec = phonoc_complex_prod(
fc3_rec_elem[n * 9 + m * 3 + l], pre_phase_factor);
fc3_reciprocal[(i * 3 + l) * num_band * num_band +
(j * 3 + m) * num_band + k * 3 + n] =
sum_lapack_complex_double(
fc3_reciprocal[(i * 3 + l) * num_band * num_band +
(j * 3 + m) * num_band + k * 3 + n],
fc3_rec);
}
}
}
}
}
static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
const double q1[3], const double q2[3],
const double *fc3,
const long is_compact_fc3,
const AtomTriplets *atom_triplets,
const long pi0, const long pi1,
const long pi2, const long leg_index) {
long i, j, k, l;
long num_satom, adrs_shift, adrs_vec1, adrs_vec2;
lapack_complex_double phase_factor, phase_factor1, phase_factor2;
@ -208,36 +241,60 @@ static void real_to_reciprocal_elements(
fc3_rec_imag[i] = 0;
}
num_satom = multi_dims[0];
num_satom = atom_triplets->multi_dims[0];
if (is_compact_fc3) {
i = pi0;
} else {
i = p2s[pi0];
i = atom_triplets->p2s_map[pi0];
}
for (j = 0; j < num_satom; j++) {
if (s2p[j] != p2s[pi1]) {
if (atom_triplets->s2p_map[j] != atom_triplets->p2s_map[pi1]) {
continue;
}
adrs_vec1 = j * multi_dims[1] + pi0;
phase_factor1 = get_phase_factor(q1, svecs, multiplicity[adrs_vec1]);
adrs_vec1 = j * atom_triplets->multi_dims[1] + pi0;
phase_factor1 = get_phase_factor(
q1, atom_triplets->svecs, atom_triplets->multiplicity[adrs_vec1]);
for (k = 0; k < num_satom; k++) {
if (s2p[k] != p2s[pi2]) {
if (atom_triplets->s2p_map[k] != atom_triplets->p2s_map[pi2]) {
continue;
}
adrs_vec2 = k * multi_dims[1] + pi0;
if (leg_index > 1) {
if (atom_triplets->all_shortest[pi0 * num_satom * num_satom +
j * num_satom + k]) {
continue;
}
}
adrs_vec2 = k * atom_triplets->multi_dims[1] + pi0;
phase_factor2 =
get_phase_factor(q2, svecs, multiplicity[adrs_vec2]);
get_phase_factor(q2, atom_triplets->svecs,
atom_triplets->multiplicity[adrs_vec2]);
adrs_shift =
i * 27 * num_satom * num_satom + j * 27 * num_satom + k * 27;
phase_factor = phonoc_complex_prod(phase_factor1, phase_factor2);
for (l = 0; l < 27; l++) {
fc3_rec_real[l] += lapack_complex_double_real(phase_factor) *
fc3[adrs_shift + l];
fc3_rec_imag[l] += lapack_complex_double_imag(phase_factor) *
fc3[adrs_shift + l];
if ((leg_index == 1) &&
(atom_triplets->all_shortest[pi0 * num_satom * num_satom +
j * num_satom + k])) {
for (l = 0; l < 27; l++) {
fc3_rec_real[l] +=
lapack_complex_double_real(phase_factor) *
fc3[adrs_shift + l] * 3;
fc3_rec_imag[l] +=
lapack_complex_double_imag(phase_factor) *
fc3[adrs_shift + l] * 3;
}
} else {
for (l = 0; l < 27; l++) {
fc3_rec_real[l] +=
lapack_complex_double_real(phase_factor) *
fc3[adrs_shift + l];
fc3_rec_imag[l] +=
lapack_complex_double_imag(phase_factor) *
fc3[adrs_shift + l];
}
}
}
}
@ -248,33 +305,26 @@ static void real_to_reciprocal_elements(
}
}
static lapack_complex_double get_pre_phase_factor(const long i_patom,
const double q_vecs[3][3],
const double (*svecs)[3],
const long multi_dims[2],
const long (*multiplicity)[2],
const long *p2s_map) {
long i, j, svecs_adrs;
double pre_phase, sum_real, sum_imag;
// This function doesn't need to think about position +
// lattice-translation because q+q'+q''=G.
static lapack_complex_double get_pre_phase_factor(
const long i_patom, const double q_vecs[3][3],
const AtomTriplets *atom_triplets) {
long j, svecs_adrs;
double pre_phase;
lapack_complex_double pre_phase_factor;
svecs_adrs = p2s_map[i_patom] * multi_dims[1];
sum_real = 0;
sum_imag = 0;
for (i = 0; i < multiplicity[svecs_adrs][0]; i++) {
pre_phase = 0;
for (j = 0; j < 3; j++) {
pre_phase += svecs[multiplicity[svecs_adrs][1] + i][j] *
(q_vecs[0][j] + q_vecs[1][j] + q_vecs[2][j]);
}
pre_phase *= M_2PI;
sum_real += cos(pre_phase);
sum_imag += sin(pre_phase);
svecs_adrs = atom_triplets->p2s_map[i_patom] * atom_triplets->multi_dims[1];
pre_phase = 0;
for (j = 0; j < 3; j++) {
pre_phase +=
atom_triplets
->svecs[atom_triplets->multiplicity[svecs_adrs][1]][j] *
(q_vecs[0][j] + q_vecs[1][j] + q_vecs[2][j]);
}
sum_real /= multiplicity[svecs_adrs][0];
sum_imag /= multiplicity[svecs_adrs][0];
pre_phase_factor = lapack_make_complex_double(sum_real, sum_imag);
pre_phase *= M_2PI;
pre_phase_factor =
lapack_make_complex_double(cos(pre_phase), sin(pre_phase));
return pre_phase_factor;
}

16
c/real_to_reciprocal.h
View File

@ -38,11 +38,19 @@
#include "lapack_wrapper.h"
#include "phonoc_array.h"
typedef struct {
const double (*svecs)[3];
long multi_dims[2];
const long (*multiplicity)[2];
const long *p2s_map;
const long *s2p_map;
long make_r0_average;
const char *all_shortest;
} AtomTriplets;
void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
const double q_vecs[3][3], const double *fc3,
const long is_compact_fc3, const double (*svecs)[3],
const long multi_dims[2],
const long (*multiplicity)[2], const long *p2s_map,
const long *s2p_map, const long make_r0_average,
const long is_compact_fc3,
const AtomTriplets *atom_triplets,
const long openmp_at_bands);
#endif

11
doc/changelog.md
View File

@ -2,6 +2,17 @@
# Change Log
## Dec-25-2023: Version 2.9.0
- Pre-release of version 3.0.
- `--v3` option enables phono3py version 3 behaviour. In phono3py-v3, it is
planned to replace $\sum_{l'l''}\Phi_{\alpha\beta\gamma}(0\kappa, l'\kappa',
l''\kappa'') \cdots$ in Eq.(41) of
<https://journals.jps.jp/doi/10.7566/JPSJ.92.012001> by
$[\sum_{l'l''}\Phi_{\alpha\beta\gamma}(0\kappa, l'\kappa', l''\kappa'') \cdots + \sum_{ll''}\Phi_{\alpha\beta\gamma}(l\kappa, 0\kappa', l''\kappa'') \cdots + \sum_{ll'}\Phi_{\alpha\beta\gamma}(l\kappa, l'\kappa', 0\kappa'') \cdots] / 3$
for better treatment of lattice sum in supercell although this requires more
computational demand. In phono3py-v3, `--v2` option will be prepared.
## Dec-4-2023: Version 2.8.0
- Maintenance release

4
doc/conf.py
View File

@ -57,9 +57,9 @@ copyright = "2015, Atsushi Togo"
# built documents.
#
# The short X.Y version.
version = "2.8"
version = "2.9"
# The full version, including alpha/beta/rc tags.
release = "2.8.0"
release = "2.9.0"
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.

11
phono3py/api_isotope.py
View File

@ -33,8 +33,6 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import warnings
import numpy as np
from phonopy.units import VaspToTHz
@ -53,7 +51,6 @@ class Phono3pyIsotope:
sigmas=None,
frequency_factor_to_THz=VaspToTHz,
use_grg=False,
store_dense_gp_map=True,
symprec=1e-5,
cutoff_frequency=None,
lapack_zheev_uplo="L",
@ -66,13 +63,6 @@ class Phono3pyIsotope:
else:
self._sigmas = sigmas
if not store_dense_gp_map:
warnings.warn(
"Phono3pyIsotope init parameter of store_dense_gp_map is deprecated. "
"This will be always set True.",
DeprecationWarning,
)
self._iso = Isotope(
mesh,
primitive,
@ -80,7 +70,6 @@ class Phono3pyIsotope:
band_indices=band_indices,
frequency_factor_to_THz=frequency_factor_to_THz,
use_grg=use_grg,
store_dense_gp_map=store_dense_gp_map,
symprec=symprec,
cutoff_frequency=cutoff_frequency,
lapack_zheev_uplo=lapack_zheev_uplo,

14
phono3py/api_jointdos.py
View File

@ -33,8 +33,6 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
import warnings
import numpy as np
from phonopy.structure.cells import Primitive, Supercell
from phonopy.structure.symmetry import Symmetry
@ -72,7 +70,6 @@ class Phono3pyJointDos:
SNF_coordinates="reciprocal",
is_mesh_symmetry=True,
is_symmetry=True,
store_dense_gp_map=True,
symprec=1e-5,
output_filename=None,
log_level=0,
@ -94,14 +91,6 @@ class Phono3pyJointDos:
self._is_mesh_symmetry = is_mesh_symmetry
self._is_symmetry = is_symmetry
if not store_dense_gp_map:
warnings.warn(
"Phono3pyJointDos init parameter of store_dense_gp_map is deprecated. "
"This will be set always True.",
DeprecationWarning,
)
self._store_dense_gp_map = store_dense_gp_map
self._use_grg = use_grg
self._SNF_coordinates = SNF_coordinates
self._symprec = symprec
@ -165,7 +154,7 @@ class Phono3pyJointDos:
use_grg=self._use_grg,
force_SNF=False,
SNF_coordinates=self._SNF_coordinates,
store_dense_gp_map=self._store_dense_gp_map,
store_dense_gp_map=True,
)
def initialize(self, mesh_numbers):
@ -180,7 +169,6 @@ class Phono3pyJointDos:
frequency_factor_to_THz=self._frequency_factor_to_THz,
frequency_scale_factor=self._frequency_scale_factor,
is_mesh_symmetry=self._is_mesh_symmetry,
store_dense_gp_map=self._store_dense_gp_map,
symprec=self._symprec,
filename=self._filename,
log_level=self._log_level,

95
phono3py/api_phono3py.py
View File

@ -34,7 +34,8 @@
# POSSIBILITY OF SUCH DAMAGE.
import warnings
from typing import Optional
from collections.abc import Sequence
from typing import Optional, Union
import numpy as np
from phonopy.exception import ForceCalculatorRequiredError
@ -152,13 +153,10 @@ class Phono3py:
is_mesh_symmetry=True,
use_grg=False,
SNF_coordinates="reciprocal",
symmetrize_fc3q=None,
store_dense_gp_map=True,
store_dense_svecs=True,
make_r0_average: bool = False,
symprec=1e-5,
calculator: Optional[str] = None,
log_level=0,
lapack_zheev_uplo=None,
):
"""Init method.
@ -212,14 +210,11 @@ class Phono3py:
`reciprocal` or `direct`. Space of coordinates to generate grid
generating matrix either in direct or reciprocal space. The default
is `reciprocal`.
symmetrize_fc3q : Deprecated.
See Phono3py.init_phph_interaction().
store_dense_gp_map : bool, optional, Deprecated.
Use dense format of BZ grid system. Default is True.
store_dense_svecs : bool, optional, Deprecated.
Shortest vectors are stored in the dense array format. This is
expected to be always True. Setting False is for rough compatibility
with v1.x. Default is True.
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.
calculator : str, optional.
@ -227,8 +222,6 @@ class Phono3py:
of physical units. Default is None, which is equivalent to "vasp".
log_level : int, optional
Verbosity control. Default is 0. This can be 0, 1, or 2.
lapack_zheev_uplo : Deprecated.
See Phono3py.init_phph_interaction().
"""
self._symprec = symprec
@ -238,21 +231,7 @@ class Phono3py:
self._use_grg = use_grg
self._SNF_coordinates = SNF_coordinates
if not store_dense_gp_map:
warnings.warn(
"Phono3py init parameter of store_dense_gp_map is deprecated. "
"This will be set always True.",
DeprecationWarning,
)
self._store_dense_gp_map = store_dense_gp_map
if not store_dense_svecs:
warnings.warn(
"Phono3py init parameter of store_dense_svecs is deprecated. "
"This will be set always True.",
DeprecationWarning,
)
self._store_dense_svecs = store_dense_svecs
self._make_r0_average = make_r0_average
self._cutoff_frequency = cutoff_frequency
self._calculator: Optional[str] = calculator
@ -360,26 +339,6 @@ class Phono3py:
)
self.sigma_cutoff = sigma_cutoff
if symmetrize_fc3q is not None:
warnings.warn(
"Phono3py init parameter of symmetrize_fc3q is deprecated. "
"Set this at Phono3py.init_phph_interaction().",
DeprecationWarning,
)
self._symmetrize_fc3q = symmetrize_fc3q
else:
self._symmetrize_fc3q = None
if lapack_zheev_uplo is not None:
warnings.warn(
"Phono3py init parameter of lapack_zheev_uplo is deprecated. "
"Set this at Phono3py.init_phph_interaction().",
DeprecationWarning,
)
self._lapack_zheev_uplo = lapack_zheev_uplo
else:
self._lapack_zheev_uplo = None
@property
def version(self):
"""Return phono3py release version number.
@ -1041,7 +1000,7 @@ class Phono3py:
return self._bz_grid.D_diag
@mesh_numbers.setter
def mesh_numbers(self, mesh_numbers):
def mesh_numbers(self, mesh_numbers: Union[int, float, Sequence, np.ndarray]):
self._set_mesh_numbers(mesh_numbers)
@property
@ -1330,8 +1289,7 @@ class Phono3py:
nac_q_direction=None,
constant_averaged_interaction=None,
frequency_scale_factor=None,
symmetrize_fc3q=False,
make_r0_average=False,
symmetrize_fc3q: bool = False,
lapack_zheev_uplo="L",
openmp_per_triplets=None,
):
@ -1366,10 +1324,6 @@ class Phono3py:
symmetrize_fc3q : bool, optional
fc3 in phonon space is symmetrized by permutation symmetry.
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.
lapack_zheev_uplo : str, optional
'L' or 'U'. Default is 'L'. This is passed to LAPACK zheev
used for phonon solver.
@ -1388,16 +1342,6 @@ class Phono3py:
msg = "Phono3py.fc2 of instance is not found."
raise RuntimeError(msg)
if self._symmetrize_fc3q is None:
_symmetrize_fc3q = symmetrize_fc3q
else:
_symmetrize_fc3q = self._symmetrize_fc3q
if self._lapack_zheev_uplo is None:
_lapack_zheev_uplo = lapack_zheev_uplo
else:
_lapack_zheev_uplo = self._lapack_zheev_uplo
self._interaction = Interaction(
self._primitive,
self._bz_grid,
@ -1409,9 +1353,9 @@ class Phono3py:
frequency_scale_factor=frequency_scale_factor,
cutoff_frequency=self._cutoff_frequency,
is_mesh_symmetry=self._is_mesh_symmetry,
symmetrize_fc3q=_symmetrize_fc3q,
make_r0_average=make_r0_average,
lapack_zheev_uplo=_lapack_zheev_uplo,
symmetrize_fc3q=symmetrize_fc3q,
make_r0_average=self._make_r0_average,
lapack_zheev_uplo=lapack_zheev_uplo,
openmp_per_triplets=openmp_per_triplets,
)
self._interaction.nac_q_direction = nac_q_direction
@ -2544,9 +2488,7 @@ class Phono3py:
else:
t_mat = np.dot(inv_supercell_matrix, primitive_matrix)
return get_primitive(
supercell, t_mat, self._symprec, store_dense_svecs=self._store_dense_svecs
)
return get_primitive(supercell, t_mat, self._symprec, store_dense_svecs=True)
def _determine_primitive_matrix(self, primitive_matrix):
pmat = get_primitive_matrix(primitive_matrix, symprec=self._symprec)
@ -2555,7 +2497,10 @@ class Phono3py:
else:
return pmat
def _set_mesh_numbers(self, mesh):
def _set_mesh_numbers(
self,
mesh: Union[int, float, Sequence, np.ndarray],
):
# initialization related to mesh
self._interaction = None
@ -2567,7 +2512,7 @@ class Phono3py:
use_grg=self._use_grg,
force_SNF=False,
SNF_coordinates=self._SNF_coordinates,
store_dense_gp_map=self._store_dense_gp_map,
store_dense_gp_map=True,
)
def _init_dynamical_matrix(self):

4
phono3py/conductivity/direct_solution.py
View File

@ -504,9 +504,7 @@ class ConductivityLBTEBase(ConductivityBase):
self._collision.run_interaction(is_full_pp=self._is_full_pp)
if self._is_full_pp and j == 0:
self._averaged_pp_interaction[
i_gp
] = self._pp.get_averaged_interaction()
self._averaged_pp_interaction[i_gp] = self._pp.averaged_interaction
for k, t in enumerate(self._temperatures):
self._collision.temperature = t

16
phono3py/conductivity/rta.py
View File

@ -115,7 +115,7 @@ class ConductivityRTABase(ConductivityBase):
log_level=log_level,
)
self._use_const_ave_pp = self._pp.get_constant_averaged_interaction()
self._use_const_ave_pp = self._pp.constant_averaged_interaction
self._read_pp = read_pp
self._store_pp = store_pp
self._pp_filename = pp_filename
@ -262,10 +262,10 @@ class ConductivityRTABase(ConductivityBase):
if self._log_level:
print(
"Constant ph-ph interaction of %6.3e is used."
% self._pp.get_constant_averaged_interaction()
% self._pp.constant_averaged_interaction
)
self._collision.run_interaction()
self._averaged_pp_interaction[i] = self._pp.get_averaged_interaction()
self._averaged_pp_interaction[i] = self._pp.averaged_interaction
elif j != 0 and (self._is_full_pp or self._sigma_cutoff is None):
if self._log_level:
print("Existing ph-ph interaction is used.")
@ -274,9 +274,7 @@ class ConductivityRTABase(ConductivityBase):
print("Calculating ph-ph interaction...")
self._collision.run_interaction(is_full_pp=self._is_full_pp)
if self._is_full_pp:
self._averaged_pp_interaction[
i
] = self._pp.get_averaged_interaction()
self._averaged_pp_interaction[i] = self._pp.averaged_interaction
# Number of triplets depends on q-point.
# So this is allocated each time.
@ -304,7 +302,7 @@ class ConductivityRTABase(ConductivityBase):
k
] = self._collision.get_detailed_imag_self_energy()
def _set_gamma_at_sigmas_lowmem(self, i, make_r0_average=False):
def _set_gamma_at_sigmas_lowmem(self, i):
"""Calculate gamma without storing ph-ph interaction strength.
`svecs` and `multi` below must not be simply replaced by
@ -386,6 +384,7 @@ class ConductivityRTABase(ConductivityBase):
self._is_N_U * 1,
self._pp.symmetrize_fc3q * 1,
self._pp.make_r0_average * 1,
self._pp.all_shortest,
self._pp.cutoff_frequency,
openmp_per_triplets * 1,
)
@ -416,11 +415,12 @@ class ConductivityRTABase(ConductivityBase):
self._is_N_U * 1,
self._pp.symmetrize_fc3q * 1,
self._pp.make_r0_average * 1,
self._pp.all_shortest,
self._pp.cutoff_frequency,
openmp_per_triplets * 1,
)
col_unit_conv = self._collision.unit_conversion_factor
pp_unit_conv = self._pp.get_unit_conversion_factor()
pp_unit_conv = self._pp.unit_conversion_factor
if self._is_N_U:
col = collisions.sum(axis=0)
col_N = collisions[0]

74
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,35 +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,
store_dense_gp_map=True,
store_dense_svecs=True,
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.
@ -229,12 +231,11 @@ def load(
cells. Default is False.
use_grg : bool, optional
Use generalized regular grid when True. Default is False.
store_dense_gp_map : bool, optional, Deprecated
Use new format of BZ grid system. Default is True.
store_dense_svecs : bool, optional, Deprecated
Shortest vectors are stored in the dense array format. This is expected
to be always True. Setting False is for rough compatibility with v1.x.
Default is True.
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
@ -305,8 +306,7 @@ def load(
is_symmetry=is_symmetry,
is_mesh_symmetry=is_mesh_symmetry,
use_grg=use_grg,
store_dense_gp_map=store_dense_gp_map,
store_dense_svecs=store_dense_svecs,
make_r0_average=make_r0_average,
calculator=calculator,
log_level=log_level,
)

16
phono3py/cui/phono3py_argparse.py
View File

@ -251,13 +251,13 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
default=None,
help="Supercell dimension for extra fc2",
)
parser.add_argument(
"--emulate-v2",
dest="emulate_v2",
action="store_true",
default=False,
help="Emulate v1.x grid system and shortest vectors.",
)
# parser.add_argument(
# "--emulate-v2",
# dest="emulate_v2",
# action="store_true",
# default=False,
# help="Emulate v2.x behaviour.",
# )
parser.add_argument(
"--factor",
dest="frequency_conversion_factor",
@ -281,7 +281,7 @@ def get_parser(fc_symmetry=False, is_nac=False, load_phono3py_yaml=False):
help="Read third order force constants",
)
parser.add_argument(
"--fc3-r0-average",
"--v3",
dest="is_fc3_r0_average",
action="store_true",
default=False,

4
phono3py/cui/phono3py_script.py
View File

@ -741,6 +741,7 @@ def init_phono3py(
is_symmetry=settings.is_symmetry,
is_mesh_symmetry=settings.is_mesh_symmetry,
use_grg=settings.use_grg,
make_r0_average=settings.is_fc3_r0_average,
symprec=symprec,
calculator=interface_mode,
log_level=log_level,
@ -1036,7 +1037,6 @@ def init_phph_interaction(
constant_averaged_interaction=ave_pp,
frequency_scale_factor=updated_settings["frequency_scale_factor"],
symmetrize_fc3q=settings.is_symmetrize_fc3_q,
make_r0_average=settings.is_fc3_r0_average,
lapack_zheev_uplo=settings.lapack_zheev_uplo,
)
@ -1446,7 +1446,7 @@ def main(**argparse_control):
is_reducible_collision_matrix=settings.is_reducible_collision_matrix, # noqa E501
is_kappa_star=settings.is_kappa_star,
gv_delta_q=settings.group_velocity_delta_q,
is_full_pp=settings.is_full_pp,
is_full_pp=(settings.is_full_pp or settings.is_symmetrize_fc3_q),
pinv_cutoff=settings.pinv_cutoff,
pinv_solver=settings.pinv_solver,
pinv_method=settings.pinv_method,

2
phono3py/cui/settings.py
View File

@ -832,7 +832,7 @@ class Phono3pyConfParser(ConfParser):
if "cutoff_pair_distance" in params:
self._settings.set_cutoff_pair_distance(params["cutoff_pair_distance"])
# Emulate v1.x grid system and shortest vectors.
# Emulate v2.x behaviour
if "emulate_v2" in params:
self._settings.set_emulate_v2(params["emulate_v2"])

3
phono3py/other/isotope.py
View File

@ -100,7 +100,6 @@ class Isotope:
bz_grid=None,
frequency_factor_to_THz=VaspToTHz,
use_grg=False,
store_dense_gp_map=True,
symprec=1e-5,
cutoff_frequency=None,
lapack_zheev_uplo="L",
@ -146,7 +145,7 @@ class Isotope:
lattice=self._primitive.cell,
symmetry_dataset=primitive_symmetry.dataset,
use_grg=use_grg,
store_dense_gp_map=store_dense_gp_map,
store_dense_gp_map=True,
)
def set_grid_point(self, grid_point):

6
phono3py/other/tetrahedron_method.py
View File

@ -36,8 +36,10 @@
import numpy as np
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phono3py.phonon.grid import BZGrid
def get_unique_grid_points(grid_points, bz_grid):
def get_unique_grid_points(grid_points, bz_grid: BZGrid):
"""Collect grid points on tetrahedron vertices around input grid points.
Find grid points of 24 tetrahedra around each grid point and
@ -88,7 +90,7 @@ def get_unique_grid_points(grid_points, bz_grid):
def get_integration_weights(
sampling_points,
grid_values,
bz_grid,
bz_grid: BZGrid,
grid_points=None,
bzgp2irgp_map=None,
function="I",

14
phono3py/phonon3/gruneisen.py
View File

@ -37,7 +37,8 @@ import sys
import numpy as np
from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
from phonopy.structure.cells import sparse_to_dense_svecs
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import Primitive
from phonopy.structure.grid_points import get_qpoints
from phonopy.units import VaspToTHz
@ -126,8 +127,8 @@ class Gruneisen:
self,
fc2,
fc3,
supercell,
primitive,
supercell: PhonopyAtoms,
primitive: Primitive,
nac_params=None,
nac_q_direction=None,
ion_clamped=False,
@ -151,12 +152,7 @@ class Gruneisen:
)
self._nac_q_direction = nac_q_direction
svecs, multi = self._pcell.get_smallest_vectors()
if self._pcell.store_dense_svecs:
self._svecs = svecs
self._multi = multi
else:
self._svecs, self._multi = sparse_to_dense_svecs(svecs, multi)
self._svecs, self._multi = self._pcell.get_smallest_vectors()
if self._ion_clamped:
num_atom_prim = len(self._pcell)

186
phono3py/phonon3/interaction.py
View File

@ -32,16 +32,15 @@
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
from __future__ import annotations
import warnings
from collections.abc import Sequence
from typing import Literal, Optional, Union
import numpy as np
from phonopy.harmonic.dynamical_matrix import get_dynamical_matrix
from phonopy.structure.cells import (
Primitive,
compute_all_sg_permutations,
sparse_to_dense_svecs,
)
from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, get_dynamical_matrix
from phonopy.structure.cells import Primitive, compute_all_sg_permutations
from phonopy.structure.symmetry import Symmetry
from phonopy.units import AMU, EV, Angstrom, Hbar, THz, VaspToTHz
@ -94,18 +93,18 @@ class Interaction:
primitive: Primitive,
bz_grid: BZGrid,
primitive_symmetry: Symmetry,
fc3=None,
band_indices=None,
constant_averaged_interaction=None,
frequency_factor_to_THz=VaspToTHz,
frequency_scale_factor=None,
unit_conversion=None,
is_mesh_symmetry=True,
symmetrize_fc3q=False,
make_r0_average=False,
cutoff_frequency=None,
lapack_zheev_uplo="L",
openmp_per_triplets=None,
fc3: Optional[np.ndarray] = None,
band_indices: Optional[Union[np.ndarray, Sequence]] = None,
constant_averaged_interaction: Optional[float] = None,
frequency_factor_to_THz: float = VaspToTHz,
frequency_scale_factor: Optional[float] = None,
unit_conversion: Optional[float] = None,
is_mesh_symmetry: bool = True,
symmetrize_fc3q: bool = False,
make_r0_average: bool = False,
cutoff_frequency: Optional[float] = None,
lapack_zheev_uplo: Literal["L", "U"] = "L",
openmp_per_triplets: Optional[bool] = None,
):
"""Init method."""
self._primitive = primitive
@ -168,18 +167,19 @@ class Interaction:
self._band_index_count = 0
svecs, multi = self._primitive.get_smallest_vectors()
if self._primitive.store_dense_svecs:
self._svecs = svecs
self._multi = multi
else:
self._svecs, self._multi = sparse_to_dense_svecs(svecs, multi)
self._svecs, self._multi = self._primitive.get_smallest_vectors()
self._masses = np.array(self._primitive.masses, dtype="double")
self._p2s = np.array(self._primitive.p2s_map, dtype="int_")
self._s2p = np.array(self._primitive.s2p_map, dtype="int_")
n_satom, n_patom, _ = self._multi.shape
self._all_shortest = np.zeros(
(n_patom, n_satom, n_satom), dtype="byte", order="C"
)
self._get_all_shortest()
def run(self, lang="C", g_zero=None):
def run(
self, lang: Literal["C", "Python"] = "C", g_zero: Optional[np.ndarray] = None
):
"""Run ph-ph interaction calculation."""
if (self._phonon_done == 0).any():
self.run_phonon_solver()
@ -203,7 +203,7 @@ class Interaction:
)
@property
def interaction_strength(self):
def interaction_strength(self) -> np.ndarray:
"""Return ph-ph interaction strength.
Returns
@ -225,7 +225,7 @@ class Interaction:
return self.interaction_strength
@property
def mesh_numbers(self):
def mesh_numbers(self) -> np.ndarray:
"""Return mesh numbers.
Returns
@ -246,12 +246,12 @@ class Interaction:
return self.mesh_numbers
@property
def is_mesh_symmetry(self):
def is_mesh_symmetry(self) -> bool:
"""Whether symmetry of grid is utilized or not."""
return self._is_mesh_symmetry
@property
def fc3(self):
def fc3(self) -> np.ndarray:
"""Return fc3."""
return self._fc3
@ -264,7 +264,7 @@ class Interaction:
return self.fc3
@property
def dynamical_matrix(self):
def dynamical_matrix(self) -> Optional[DynamicalMatrix]:
"""Return DynamicalMatrix class instance."""
return self._dm
@ -278,7 +278,7 @@ class Interaction:
return self.dynamical_matrix
@property
def primitive(self):
def primitive(self) -> Primitive:
"""Return Primitive class instance."""
return self._primitive
@ -292,11 +292,13 @@ class Interaction:
return self.primitive
@property
def primitive_symmetry(self):
def primitive_symmetry(self) -> Symmetry:
"""Return Symmetry class instance of primitive cell."""
return self._primitive_symmetry
def get_triplets_at_q(self):
def get_triplets_at_q(
self,
) -> tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray):
"""Return grid point triplets information.
triplets_at_q is in BZ-grid.
@ -313,12 +315,12 @@ class Interaction:
)
@property
def bz_grid(self):
def bz_grid(self) -> BZGrid:
"""Return BZGrid class instance."""
return self._bz_grid
@property
def band_indices(self):
def band_indices(self) -> np.ndarray:
"""Return band indices.
Returns
@ -339,12 +341,12 @@ class Interaction:
return self.band_indices
@property
def nac_params(self):
def nac_params(self) -> dict:
"""Return NAC params."""
return self._nac_params
@property
def nac_q_direction(self):
def nac_q_direction(self) -> Optional[np.ndarray]:
"""Return q-direction used for NAC at q->0.
Direction of q-vector watching from Gamma point used for
@ -383,7 +385,7 @@ class Interaction:
self.nac_q_direction = nac_q_direction
@property
def zero_value_positions(self):
def zero_value_positions(self) -> Optional[np.ndarray]:
"""Return zero ph-ph interaction elements information.
Returns
@ -402,7 +404,7 @@ class Interaction:
)
return self.zero_value_positions
def get_phonons(self):
def get_phonons(self) -> tuple(np.ndarray, np.ndarray, np.ndarray):
"""Return phonons on grid.
Returns
@ -423,7 +425,7 @@ class Interaction:
return self._frequencies, self._eigenvectors, self._phonon_done
@property
def frequency_factor_to_THz(self):
def frequency_factor_to_THz(self) -> float:
"""Return phonon frequency conversion factor to THz."""
return self._frequency_factor_to_THz
@ -437,7 +439,7 @@ class Interaction:
return self.frequency_factor_to_THz
@property
def lapack_zheev_uplo(self):
def lapack_zheev_uplo(self) -> Literal["L", "U"]:
"""Return U or L for lapack zheev solver."""
return self._lapack_zheev_uplo
@ -451,7 +453,7 @@ class Interaction:
return self.lapack_zheev_uplo
@property
def cutoff_frequency(self):
def cutoff_frequency(self) -> float:
"""Return cutoff phonon frequency to judge imaginary phonon."""
return self._cutoff_frequency
@ -465,17 +467,17 @@ class Interaction:
return self.cutoff_frequency
@property
def openmp_per_triplets(self):
def openmp_per_triplets(self) -> bool:
"""Return whether OpenMP distribution over triplets or bands."""
return self._openmp_per_triplets
@property
def symmetrize_fc3q(self):
def symmetrize_fc3q(self) -> bool:
"""Return boolean of symmetrize_fc3q."""
return self._symmetrize_fc3q
@property
def make_r0_average(self):
def make_r0_average(self) -> bool:
"""Return boolean of make_r0_average.
This flag is used to activate averaging of fc3 transformation
@ -486,7 +488,20 @@ class Interaction:
"""
return self._make_r0_average
def get_averaged_interaction(self):
@property
def all_shortest(self) -> np.ndarray:
"""Return boolean of make_r0_average.
This flag is used to activate averaging of fc3 transformation
from real space to reciprocal space around three atoms. With False,
it is done at the first atom. With True, it is done at three atoms
and averaged.
"""
return self._all_shortest
@property
def averaged_interaction(self) -> np.ndarray:
"""Return sum over phonon triplets of interaction strength.
See Eq.(21) of PRB 91, 094306 (2015)
@ -498,18 +513,49 @@ class Interaction:
v_sum = np.dot(w, v.sum(axis=2).sum(axis=2))
return v_sum / np.prod(v.shape[2:])
def get_primitive_and_supercell_correspondence(self):
def get_averaged_interaction(self):
"""Return sum over phonon triplets of interaction strength."""
warnings.warn(
"Use attribute, Interaction.averaged_interaction "
"instead of Interaction.get_averaged_interaction().",
DeprecationWarning,
)
return self.averaged_interaction
def get_primitive_and_supercell_correspondence(
self,
) -> tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray):
"""Return atomic pair information."""
return (self._svecs, self._multi, self._p2s, self._s2p, self._masses)
def get_unit_conversion_factor(self):
@property
def unit_conversion_factor(self) -> float:
"""Return unit conversion factor."""
return self._unit_conversion
def get_constant_averaged_interaction(self):
def get_unit_conversion_factor(self):
"""Return unit conversion factor."""
warnings.warn(
"Use attribute, Interaction.unit_conversion_factor "
"instead of Interaction.get_unit_conversion_factor().",
DeprecationWarning,
)
return self.unit_conversion_factor
@property
def constant_averaged_interaction(self) -> float:
"""Return constant averaged interaction."""
return self._constant_averaged_interaction
def get_constant_averaged_interaction(self):
"""Return constant averaged interaction."""
warnings.warn(
"Use attribute, Interaction.constant_averaged_interaction "
"instead of Interaction.get_constant_averaged_interaction().",
DeprecationWarning,
)
return self.constant_averaged_interaction
def set_interaction_strength(self, pp_strength, g_zero=None):
"""Set interaction strength."""
self._interaction_strength = pp_strength
@ -868,7 +914,7 @@ class Interaction:
else:
self._band_indices = np.array(band_indices, dtype="int_")
def _run_c(self, g_zero, make_r0_average=False):
def _run_c(self, g_zero):
import phono3py._phono3py as phono3c
num_band = len(self._primitive) * 3
@ -907,6 +953,7 @@ class Interaction:
self._band_indices,
self._symmetrize_fc3q * 1,
self._make_r0_average * 1,
self._all_shortest,
self._cutoff_frequency,
openmp_per_triplets * 1,
)
@ -984,6 +1031,43 @@ class Interaction:
self._eigenvectors_at_gamma = self._eigenvectors[gp_Gamma].copy()
self._phonon_done[gp_Gamma] = 0
def _get_all_shortest(self):
"""Return array indicating distances among three atoms are all shortest.
multi.shape = (n_satom, n_patom)
svecs : distance with respect to primitive cell basis
perms.shape = (n_pure_trans, n_satom)
"""
svecs = self._svecs
multi = self._multi
n_satom, n_patom, _ = multi.shape
perms = self._primitive.atomic_permutations
s2pp_map = [self._primitive.p2p_map[i] for i in self._s2p]
lattice = self._primitive.cell
for i_patom, j_atom in np.ndindex((n_patom, n_satom)):
if multi[j_atom, i_patom, 0] > 1:
continue
j_patom = s2pp_map[j_atom]
i_perm = np.where(perms[:, j_atom] == self._p2s[j_patom])[0]
assert len(i_perm) == 1
for k_atom in range(n_satom):
if multi[k_atom, i_patom, 0] > 1:
continue
k_atom_mapped = perms[i_perm[0], k_atom]
if multi[k_atom_mapped, j_patom, 0] > 1:
continue
vec_jk = (
svecs[multi[k_atom, i_patom, 1]] - svecs[multi[j_atom, i_patom, 1]]
)
d_jk = np.linalg.norm(vec_jk @ lattice)
d_jk_mapped = np.linalg.norm(
svecs[multi[k_atom_mapped, j_patom, 1]] @ lattice
)
if abs(d_jk_mapped - d_jk) < self._symprec:
self._all_shortest[i_patom, j_atom, k_atom] = 1
def all_bands_exist(interaction: Interaction):
"""Return if all bands are selected or not."""

2
phono3py/phonon3/joint_dos.py
View File

@ -68,7 +68,6 @@ class JointDos:
frequency_factor_to_THz=VaspToTHz,
frequency_scale_factor=1.0,
is_mesh_symmetry=True,
store_dense_gp_map=False,
symprec=1e-5,
filename=None,
log_level=False,
@ -93,7 +92,6 @@ class JointDos:
self._frequency_factor_to_THz = frequency_factor_to_THz
self._frequency_scale_factor = frequency_scale_factor
self._is_mesh_symmetry = is_mesh_symmetry
self._store_dense_gp_map = store_dense_gp_map
self._symprec = symprec
self._filename = filename
self._log_level = log_level

8
phono3py/phonon3/real_to_reciprocal.py
View File

@ -34,7 +34,6 @@
# POSSIBILITY OF SUCH DAMAGE.
import numpy as np
from phonopy.structure.cells import sparse_to_dense_svecs
class RealToReciprocal:
@ -50,12 +49,7 @@ class RealToReciprocal:
self._p2s_map = primitive.p2s_map
self._s2p_map = primitive.s2p_map
# Reduce supercell atom index to primitive index
svecs, multi = self._primitive.get_smallest_vectors()
if self._primitive.store_dense_svecs:
self._svecs = svecs
self._multi = multi
else:
self._svecs, self._multi = sparse_to_dense_svecs(svecs, multi)
self._svecs, self._multi = self._primitive.get_smallest_vectors()
self._fc3_reciprocal = None
def run(self, triplet):

7
phono3py/phonon3/triplets.py
View File

@ -70,8 +70,6 @@ def get_triplets_at_q(
Inversion symemtry is added if it doesn't exist. Default is True.
swappable : bool, optional
q1 and q2 among (q0, q1, q2) can be swapped. Deafult is True.
store_dense_gp_map : bool, optional
See the detail in the docstring of ``_relocate_BZ_grid_address``.
Returns
-------
@ -343,7 +341,10 @@ def _get_BZ_triplets_at_q(bz_grid_index, bz_grid: BZGrid, map_triplets):
def _set_triplets_integration_weights_c(
g, g_zero, pp: Union["Interaction", "JointDos"], frequency_points
g: np.ndarray,
g_zero: np.ndarray,
pp: Union["Interaction", "JointDos"],
frequency_points,
):
import phono3py._phono3py as phono3c

2
phono3py/version.py
View File

@ -33,4 +33,4 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
__version__ = "2.8.0"
__version__ = "2.9.0"

47
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,26 +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]
aln_lda.init_phph_interaction(make_r0_average=True)
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,
]
@ -88,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)
'''

359
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,182 +214,119 @@ def test_kappa_RTA_si_N_U(si_pbesol):
is_N_U=is_N_U,
)
gN, gU = ph3.thermal_conductivity.get_gamma_N_U()
# for line in gN.reshape(-1, 6):
# print("[", ",".join([f"{val:.8f}" for val in line]), "],")
# for line in gU.reshape(-1, 6):
# print("[", ",".join([f"{val:.8f}" for val in line]), "],")
gN_ref = [
0.00000000,
0.00000000,
0.00000000,
0.07402084,
0.07402084,
0.07402084,
0.00078535,
0.00078535,
0.00917995,
0.02178049,
0.04470075,
0.04470075,
0.00173337,
0.00173337,
0.01240191,
0.00198981,
0.03165195,
0.03165195,
0.00224713,
0.00224713,
0.00860026,
0.03083611,
0.03083611,
0.02142118,
0.00277534,
0.00330170,
0.02727451,
0.00356415,
0.01847744,
0.01320643,
0.00155072,
0.00365611,
0.01641919,
0.00650083,
0.02576069,
0.01161589,
0.00411969,
0.00411969,
0.00168211,
0.00168211,
0.01560092,
0.01560092,
0.00620091,
0.00620091,
0.03764912,
0.03764912,
0.02668523,
0.02668523,
]
gU_ref = [
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00015178,
0.00015178,
0.00076936,
0.00727539,
0.00113112,
0.00113112,
0.00022696,
0.00022696,
0.00072558,
0.00000108,
0.00021968,
0.00021968,
0.00079397,
0.00079397,
0.00111068,
0.00424761,
0.00424761,
0.00697760,
0.00221593,
0.00259510,
0.01996296,
0.00498962,
0.01258375,
0.00513825,
0.00148802,
0.00161955,
0.01589219,
0.00646134,
0.00577275,
0.00849711,
0.00313208,
0.00313208,
0.00036610,
0.00036610,
0.01135335,
0.01135335,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
]
if si_pbesol._make_r0_average:
gN_ref = [
[0.00000000, 0.00000000, 0.00000000, 0.07898606, 0.07898606, 0.07898606],
[0.00079647, 0.00079647, 0.00913611, 0.01911102, 0.04553001, 0.04553001],
[0.00173868, 0.00173868, 0.01404937, 0.00201732, 0.03354033, 0.03354033],
[0.00223616, 0.00223616, 0.01039331, 0.02860916, 0.02860916, 0.01987485],
[0.00291788, 0.00356241, 0.02858543, 0.00367742, 0.02065990, 0.01533763],
[0.00146333, 0.00343175, 0.01596851, 0.00626596, 0.02431620, 0.01091592],
[0.00396766, 0.00396766, 0.00159161, 0.00159161, 0.01479018, 0.01479018],
[0.00682740, 0.00682740, 0.03983399, 0.03983399, 0.02728522, 0.02728522],
]
gU_ref = [
[0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
[0.00015184, 0.00015184, 0.00075965, 0.00736940, 0.00114177, 0.00114177],
[0.00022400, 0.00022400, 0.00072237, 0.00000112, 0.00022016, 0.00022016],
[0.00079188, 0.00079188, 0.00106579, 0.00418717, 0.00418717, 0.00712761],
[0.00219252, 0.00262840, 0.01927670, 0.00491388, 0.01254730, 0.00519414],
[0.00146999, 0.00168024, 0.01596274, 0.00641979, 0.00597353, 0.00859841],
[0.00307881, 0.00307881, 0.00036554, 0.00036554, 0.01176737, 0.01176737],
[0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
]
else:
gN_ref = [
[0.00000000, 0.00000000, 0.00000000, 0.07832198, 0.07832198, 0.07832198],
[0.00079578, 0.00079578, 0.00909025, 0.01917012, 0.04557656, 0.04557656],
[0.00176235, 0.00176235, 0.01414436, 0.00204092, 0.03361112, 0.03361112],
[0.00221919, 0.00221919, 0.01020133, 0.02889554, 0.02889554, 0.01995543],
[0.00292189, 0.00356099, 0.02855954, 0.00370530, 0.02071850, 0.01533334],
[0.00147656, 0.00342335, 0.01589430, 0.00630792, 0.02427768, 0.01099287],
[0.00400675, 0.00400675, 0.00162186, 0.00162186, 0.01478489, 0.01478489],
[0.00676576, 0.00676576, 0.03984290, 0.03984290, 0.02715102, 0.02715102],
]
gU_ref = [
[0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
[0.00015178, 0.00015178, 0.00076936, 0.00727539, 0.00113112, 0.00113112],
[0.00022696, 0.00022696, 0.00072558, 0.00000108, 0.00021968, 0.00021968],
[0.00079397, 0.00079397, 0.00111068, 0.00424761, 0.00424761, 0.00697760],
[0.00219456, 0.00261878, 0.01928629, 0.00490046, 0.01249235, 0.00517685],
[0.00149539, 0.00161230, 0.01594274, 0.00653088, 0.00593572, 0.00849890],
[0.00311169, 0.00311169, 0.00036610, 0.00036610, 0.01171667, 0.01171667],
[0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
]
# print(np.sum(gN), np.sum(gU))
np.testing.assert_allclose(
np.sum([gN_ref, gU_ref], axis=0), gN.ravel() + gU.ravel(), atol=1e-2
)
np.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
def _get_kappa(
ph3,
ph3: Phono3py,
mesh,
is_isotope=False,
is_full_pp=False,
openmp_per_triplets=None,
make_r0_average=False,
):
ph3.mesh_numbers = mesh
ph3.init_phph_interaction(
make_r0_average=make_r0_average, openmp_per_triplets=openmp_per_triplets
)
ph3.init_phph_interaction(openmp_per_triplets=openmp_per_triplets)
ph3.run_thermal_conductivity(
temperatures=[
300,

373
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()
gN_ref = [
0.00000000,
0.00000000,
0.00000000,
0.07402084,
0.07402084,
0.07402084,
0.00078535,
0.00078535,
0.00917995,
0.02178049,
0.04470075,
0.04470075,
0.00173337,
0.00173337,
0.01240191,
0.00198981,
0.03165195,
0.03165195,
0.00224713,
0.00224713,
0.00860026,
0.03083611,
0.03083611,
0.02142118,
0.00277534,
0.00330170,
0.02727451,
0.00356415,
0.01847744,
0.01320643,
0.00155072,
0.00365611,
0.01641919,
0.00650083,
0.02576069,
0.01161589,
0.00411969,
0.00411969,
0.00168211,
0.00168211,
0.01560092,
0.01560092,
0.00620091,
0.00620091,
0.03764912,
0.03764912,
0.02668523,
0.02668523,
]
gU_ref = [
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00015178,
0.00015178,
0.00076936,
0.00727539,
0.00113112,
0.00113112,
0.00022696,
0.00022696,
0.00072558,
0.00000108,
0.00021968,
0.00021968,
0.00079397,
0.00079397,
0.00111068,
0.00424761,
0.00424761,
0.00697760,
0.00221593,
0.00259510,
0.01996296,
0.00498962,
0.01258375,
0.00513825,
0.00148802,
0.00161955,
0.01589219,
0.00646134,
0.00577275,
0.00849711,
0.00313208,
0.00313208,
0.00036610,
0.00036610,
0.01135335,
0.01135335,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
0.00000000,
]
# for line in gU.reshape(-1, 6):
# print("[", ",".join([f"{val:.8f}" for val in line]), "],")
if si_pbesol._make_r0_average:
gN_ref = [
[0.0, 0.0, 0.0, 0.07898606, 0.07898606, 0.07898606],
[0.00079647, 0.00079647, 0.00913611, 0.01911102, 0.04553001, 0.04553001],
[0.00173868, 0.00173868, 0.01404937, 0.00201732, 0.03354033, 0.03354033],
[0.00223616, 0.00223616, 0.01039331, 0.02860916, 0.02860916, 0.01987485],
[0.00291788, 0.00356241, 0.02858543, 0.00367742, 0.0206599, 0.01533763],
[0.00146333, 0.00343175, 0.01596851, 0.00626596, 0.0243162, 0.01091592],
[0.00396766, 0.00396766, 0.00159161, 0.00159161, 0.01479018, 0.01479018],
[0.0068274, 0.0068274, 0.03983399, 0.03983399, 0.02728522, 0.02728522],
]
gU_ref = [
[0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
[0.00015184, 0.00015184, 0.00075965, 0.00736940, 0.00114177, 0.00114177],
[0.00022400, 0.00022400, 0.00072237, 0.00000112, 0.00022016, 0.00022016],
[0.00079188, 0.00079188, 0.00106579, 0.00418717, 0.00418717, 0.00712761],
[0.00219252, 0.00262840, 0.01927670, 0.00491388, 0.01254730, 0.00519414],
[0.00146999, 0.00168024, 0.01596274, 0.00641979, 0.00597353, 0.00859841],
[0.00307881, 0.00307881, 0.00036554, 0.00036554, 0.01176737, 0.01176737],
[0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
]
else:
gN_ref = [
[0.0, 0.0, 0.0, 0.07832198, 0.07832198, 0.07832198],
[0.00079578, 0.00079578, 0.00909025, 0.01917012, 0.04557656, 0.04557656],
[0.00176235, 0.00176235, 0.01414436, 0.00204092, 0.03361112, 0.03361112],
[0.00221919, 0.00221919, 0.01020133, 0.02889554, 0.02889554, 0.01995543],
[0.00292189, 0.00356099, 0.02855954, 0.0037053, 0.0207185, 0.01533334],
[0.00147656, 0.00342335, 0.0158943, 0.00630792, 0.02427768, 0.01099287],
[0.00400675, 0.00400675, 0.00162186, 0.00162186, 0.01478489, 0.01478489],
[0.00676576, 0.00676576, 0.0398429, 0.0398429, 0.02715102, 0.02715102],
]
gU_ref = [
[0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
[0.00015178, 0.00015178, 0.00076936, 0.00727539, 0.00113112, 0.00113112],
[0.00022696, 0.00022696, 0.00072558, 0.00000108, 0.00021968, 0.00021968],
[0.00079397, 0.00079397, 0.00111068, 0.00424761, 0.00424761, 0.00697760],
[0.00219456, 0.00261878, 0.01928629, 0.00490046, 0.01249235, 0.00517685],
[0.00149539, 0.00161230, 0.01594274, 0.00653088, 0.00593572, 0.00849890],
[0.00311169, 0.00311169, 0.00036610, 0.00036610, 0.01171667, 0.01171667],
[0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000],
]
# print(np.sum(gN), np.sum(gU))
np.testing.assert_allclose(
np.sum([gN_ref, gU_ref], axis=0), gN.ravel() + gU.ravel(), atol=1e-2
)
np.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(

111
test/conftest.py
View File

@ -12,16 +12,15 @@ import phono3py
from phono3py import Phono3py
cwd = Path(__file__).parent
store_dense_gp_map = True
def pytest_addoption(parser):
"""Activate v1 emulation with --v1 option."""
"""Activate v2 emulation with --v2 option."""
parser.addoption(
"--v1",
action="store_false",
default=True,
help="Run with phono3py v1.x emulation.",
"--v2",
action="store_true",
default=False,
help="Run with phono3py v2.x emulation.",
)
@ -59,12 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
forces_fc3_filename=forces_fc3_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
make_r0_average=not enable_v2,
log_level=1,
)
@ -80,13 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
forces_fc3_filename=forces_fc3_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
use_grg=True,
make_r0_average=not enable_v2,
log_level=1,
)
@ -101,14 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
forces_fc3_filename=forces_fc3_filename,
is_symmetry=False,
produce_fc=False,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
make_r0_average=not enable_v2,
log_level=1,
)
@ -123,14 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
forces_fc3_filename=forces_fc3_filename,
is_mesh_symmetry=False,
produce_fc=False,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
make_r0_average=not enable_v2,
log_level=1,
)
@ -145,13 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
forces_fc3_filename=forces_fc3_filename,
is_compact_fc=True,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
make_r0_average=not enable_v2,
log_level=1,
)
@ -165,11 +159,10 @@ def si_pbesol_111(request) -> Phono3py:
"""
yaml_filename = cwd / "phono3py_params_Si111.yaml"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
make_r0_average=not enable_v2,
log_level=1,
)
@ -186,12 +179,11 @@ def si_pbesol_111_alm(request) -> Phono3py:
pytest.importorskip("alm")
yaml_filename = cwd / "phono3py_params_Si111.yaml"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
fc_calculator="alm",
make_r0_average=not enable_v2,
log_level=1,
)
@ -218,11 +210,10 @@ def si_pbesol_111_222_fd(request) -> Phono3py:
"""
yaml_filename = cwd / "phono3py_params_Si-111-222.yaml"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
make_r0_average=not enable_v2,
log_level=1,
)
@ -239,12 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
fc_calculator="alm",
make_r0_average=not enable_v2,
log_level=1,
)
@ -261,12 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
fc_calculator="alm|",
make_r0_average=not enable_v2,
log_level=1,
)
@ -283,12 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
fc_calculator="|alm",
make_r0_average=not enable_v2,
log_level=1,
)
@ -306,13 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
fc_calculator="alm",
fc_calculator_options="cutoff = 3",
make_r0_average=not enable_v2,
log_level=1,
)
@ -330,13 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
fc_calculator="alm",
fc_calculator_options="cutoff = 3|",
make_r0_average=not enable_v2,
log_level=1,
)
@ -354,13 +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("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
fc_calculator="alm",
fc_calculator_options="|cutoff = 3",
make_r0_average=not enable_v2,
log_level=1,
)
@ -374,11 +359,10 @@ def nacl_pbe(request) -> Phono3py:
"""
yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
make_r0_average=not enable_v2,
log_level=1,
)
@ -392,12 +376,11 @@ def nacl_pbe_compact_fc(request) -> Phono3py:
"""
yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
is_compact_fc=True,
make_r0_average=not enable_v2,
log_level=1,
)
@ -410,12 +393,11 @@ def nacl_pbe_cutoff_fc3(request) -> Phono3py:
"""
yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
ph3 = phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
produce_fc=False,
make_r0_average=not enable_v2,
log_level=1,
)
forces = ph3.forces
@ -444,12 +426,11 @@ def nacl_pbe_cutoff_fc3_all_forces(request) -> Phono3py:
"""
yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
ph3 = phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
produce_fc=False,
make_r0_average=not enable_v2,
log_level=1,
)
forces = ph3.forces
@ -469,12 +450,11 @@ def nacl_pbe_cutoff_fc3_compact_fc(request) -> Phono3py:
"""
yaml_filename = cwd / "phono3py_params_NaCl222.yaml.xz"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
ph3 = phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
produce_fc=False,
make_r0_average=not enable_v2,
log_level=1,
)
forces = ph3.forces
@ -493,11 +473,10 @@ def aln_lda(request) -> Phono3py:
"""
yaml_filename = cwd / "phono3py_params_AlN332.yaml.xz"
enable_v2 = request.config.getoption("--v1")
enable_v2 = request.config.getoption("--v2")
return phono3py.load(
yaml_filename,
store_dense_gp_map=enable_v2,
store_dense_svecs=enable_v2,
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(

2506
test/phonon3/test_imag_self_energy.py
View File

File diff suppressed because it is too large Load Diff

386
test/phonon3/test_interaction.py
View File

@ -1,191 +1,201 @@
"""Test Interaction class."""
from __future__ import annotations
from collections.abc import Sequence
from typing import Literal, Optional, Union
import numpy as np
import pytest
from phonopy.structure.cells import get_smallest_vectors
from phono3py import Phono3py
from phono3py.phonon3.interaction import Interaction
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 = [
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", "Py"])
def test_interaction_RTA_si(si_pbesol, lang):
@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,
4.744361e-08,
4.832248e-08,
4.698535e-08,
4.597876e-08,
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 = _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
itr.interaction_strength.sum(axis=(1, 2, 3)), ref_itr_RTA_Si, rtol=0, atol=1e-10
)
def test_interaction_RTA_AlN(aln_lda):
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,
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,
]
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 = _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):
"""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,
)
@ -271,7 +281,7 @@ def test_interaction_run_phonon_solver_at_gamma_NaCl(nacl_pbe: Phono3py):
)
def test_phonon_solver_expand_RTA_si(si_pbesol):
def test_phonon_solver_expand_RTA_si(si_pbesol: Phono3py):
"""Test phonon solver with eigenvector rotation of Si.
Eigenvectors can be different but frequencies must be almost the same.
@ -286,19 +296,53 @@ def test_phonon_solver_expand_RTA_si(si_pbesol):
np.testing.assert_allclose(freqs, freqs_expanded, rtol=0, atol=1e-6)
def test_get_all_shortest(aln_lda: Phono3py):
"""Test Interaction._get_all_shortest."""
ph3 = aln_lda
ph3.mesh_numbers = 30
itr = Interaction(
ph3.primitive,
ph3.grid,
ph3.primitive_symmetry,
cutoff_frequency=1e-5,
)
s_svecs, s_multi = get_smallest_vectors(
ph3.supercell.cell,
ph3.supercell.scaled_positions,
ph3.supercell.scaled_positions,
store_dense_svecs=True,
)
s_lattice = ph3.supercell.cell
p_lattice = itr.primitive.cell
shortests = itr._all_shortest
svecs, multi, _, _, _ = itr.get_primitive_and_supercell_correspondence()
n_satom, n_patom, _ = multi.shape
for i, j, k in np.ndindex((n_patom, n_satom, n_satom)):
is_found = 0
if multi[j, i, 0] == 1 and multi[k, i, 0] == 1 and s_multi[j, k, 0] == 1:
d_jk_shortest = np.linalg.norm(s_svecs[s_multi[j, k, 1]] @ s_lattice)
vec_ij = svecs[multi[j, i, 1]]
vec_ik = svecs[multi[k, i, 1]]
vec_jk = vec_ik - vec_ij
d_jk = np.linalg.norm(vec_jk @ p_lattice)
if abs(d_jk - d_jk_shortest) < ph3.symmetry.tolerance:
is_found = 1
assert shortests[i, j, k] == is_found
def _get_irt(
ph3: Phono3py,
mesh,
nac_params=None,
solve_dynamical_matrices=True,
make_r0_average=False,
mesh: Union[int, float, Sequence, np.ndarray],
nac_params: Optional[dict] = None,
solve_dynamical_matrices: bool = True,
make_r0_average: bool = False,
):
ph3.mesh_numbers = mesh
itr = Interaction(
ph3.primitive,
ph3.grid,
ph3.primitive_symmetry,
ph3.fc3,
fc3=ph3.fc3,
make_r0_average=make_r0_average,
cutoff_frequency=1e-4,
)
@ -320,7 +364,7 @@ def _get_irt(
return itr
def _show(itr):
def _show(itr: Interaction):
itr_vals = itr.interaction_strength.sum(axis=(1, 2, 3))
for i, v in enumerate(itr_vals):
print("%e, " % v, end="")

86
test/phonon3/test_joint_dos.py
View File

@ -205,14 +205,8 @@ nacl_jdos_gamma_at_300K = [
]
@pytest.mark.parametrize("gp,store_dense_gp_map", [(105, True)])
def test_jdos_si(si_pbesol: Phono3py, gp: int, store_dense_gp_map: bool):
"""Test joint-DOS by Si.
store_dense_gp_map=False : 103
store_dense_gp_map=True : 105
"""
def test_jdos_si(si_pbesol: Phono3py):
"""Test joint-DOS by Si."""
si_pbesol.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
si_pbesol.phonon_supercell,
@ -220,10 +214,9 @@ def test_jdos_si(si_pbesol: Phono3py, gp: int, store_dense_gp_map: bool):
si_pbesol.fc2,
mesh=si_pbesol.mesh_numbers,
num_frequency_points=10,
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([gp])
jdos.run([105])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(si_freq_points, jdos.frequency_points, atol=1e-5)
@ -233,8 +226,7 @@ def test_jdos_si(si_pbesol: Phono3py, gp: int, store_dense_gp_map: bool):
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(105, True)])
def test_jdso_si_nomeshsym(si_pbesol: Phono3py, gp: int, store_dense_gp_map: bool):
def test_jdso_si_nomeshsym(si_pbesol: Phono3py):
"""Test joint-DOS without considering mesh symmetry by Si."""
si_pbesol.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
@ -244,10 +236,9 @@ def test_jdso_si_nomeshsym(si_pbesol: Phono3py, gp: int, store_dense_gp_map: boo
mesh=si_pbesol.mesh_numbers,
num_frequency_points=10,
is_mesh_symmetry=False,
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([gp])
jdos.run([105])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(si_freq_points, jdos.frequency_points, atol=1e-5)
# print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
@ -256,8 +247,7 @@ def test_jdso_si_nomeshsym(si_pbesol: Phono3py, gp: int, store_dense_gp_map: boo
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(105, True)])
def test_jdos_nacl(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
def test_jdos_nacl(nacl_pbe: Phono3py):
"""Test joint-DOS by NaCl."""
nacl_pbe.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
@ -267,10 +257,9 @@ def test_jdos_nacl(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
mesh=nacl_pbe.mesh_numbers,
nac_params=nacl_pbe.nac_params,
num_frequency_points=10,
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([gp])
jdos.run([105])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(nacl_freq_points, jdos.frequency_points, atol=1e-5)
# print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
@ -279,8 +268,7 @@ def test_jdos_nacl(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(0, True)])
def test_jdos_nacl_gamma(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
def test_jdos_nacl_gamma(nacl_pbe: Phono3py):
"""Test joint-DOS at Gamma-point by NaCl."""
nacl_pbe.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
@ -291,10 +279,9 @@ def test_jdos_nacl_gamma(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
nac_params=nacl_pbe.nac_params,
nac_q_direction=[1, 0, 0],
num_frequency_points=10,
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([gp])
jdos.run([0])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(nacl_freq_points_gamma, jdos.frequency_points, atol=1e-5)
# print(", ".join(["%.7f" % jd for jd in jdos.joint_dos.ravel()]))
@ -303,8 +290,7 @@ def test_jdos_nacl_gamma(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(105, True)])
def test_jdos_nacl_at_300K(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
def test_jdos_nacl_at_300K(nacl_pbe: Phono3py):
"""Test joint-DOS at 300K by NaCl."""
nacl_pbe.mesh_numbers = [9, 9, 9]
jdos = Phono3pyJointDos(
@ -317,10 +303,9 @@ def test_jdos_nacl_at_300K(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool
temperatures=[
300,
],
store_dense_gp_map=store_dense_gp_map,
log_level=1,
)
jdos.run([gp])
jdos.run([105])
# print(", ".join(["%.7f" % fp for fp in jdos.frequency_points]))
np.testing.assert_allclose(
nacl_freq_points_at_300K, jdos.frequency_points, atol=1e-5
@ -363,29 +348,22 @@ def test_jdos_nacl_nac_gamma_at_300K_npoints(nacl_pbe: Phono3py):
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(0, False), (0, True)])
def test_jdos_nac_direction_phonon_NaCl(
nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool
):
def test_jdos_nac_direction_phonon_NaCl(nacl_pbe: Phono3py):
"""Test JDOS of NaCl with nac_q_direction."""
jdos = _get_jdos(
nacl_pbe,
[7, 7, 7],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.nac_q_direction = [1, 0, 0]
jdos.set_grid_point(gp)
jdos.set_grid_point(0)
frequencies, _, _ = jdos.get_phonons()
np.testing.assert_allclose(
frequencies[0], [0, 0, 0, 4.59488262, 4.59488262, 7.41183870], rtol=0, atol=1e-6
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(0, False), (0, True)])
def test_jdos_nac_direction_phonon_NaCl_second_error(
nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool
):
def test_jdos_nac_direction_phonon_NaCl_second_error(nacl_pbe: Phono3py):
"""Test JDOS of NaCl with nac_q_direction.
Second setting non-gamma grid point must raise exception.
@ -395,18 +373,14 @@ def test_jdos_nac_direction_phonon_NaCl_second_error(
nacl_pbe,
[7, 7, 7],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.nac_q_direction = [1, 0, 0]
jdos.set_grid_point(gp)
jdos.set_grid_point(0)
with pytest.raises(RuntimeError):
jdos.set_grid_point(1)
@pytest.mark.parametrize("gp1,gp2,store_dense_gp_map", [(0, 1, False), (0, 1, True)])
def test_jdos_nac_direction_phonon_NaCl_second_no_error(
nacl_pbe: Phono3py, gp1: int, gp2: int, store_dense_gp_map: bool
):
def test_jdos_nac_direction_phonon_NaCl_second_no_error(nacl_pbe: Phono3py):
"""Test JDOS of NaCl with nac_q_direction.
Second setting non-gamma grid point should not raise exception because
@ -418,28 +392,25 @@ def test_jdos_nac_direction_phonon_NaCl_second_no_error(
nacl_pbe,
[7, 7, 7],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.nac_q_direction = [1, 0, 0]
jdos.set_grid_point(gp1)
jdos.set_grid_point(0)
jdos.nac_q_direction = None
jdos.set_grid_point(gp2)
jdos.set_grid_point(1)
frequencies, _, _ = jdos.get_phonons()
np.testing.assert_allclose(
frequencies[0], [0, 0, 0, 4.59488262, 4.59488262, 4.59488262], rtol=0, atol=1e-6
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdos_nac_NaCl_300K_C(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
def test_jdos_nac_NaCl_300K_C(nacl_pbe: Phono3py):
"""Test running JDOS of NaCl in C mode."""
jdos = _get_jdos(
nacl_pbe,
[9, 9, 9],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.set_grid_point(gp)
jdos.set_grid_point(105)
jdos.frequency_points = nacl_freq_points_at_300K
jdos.temperature = 300
jdos.run_phonon_solver()
@ -450,16 +421,14 @@ def test_jdos_nac_NaCl_300K_C(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: b
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdos_nac_NaCl_300K_Py(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
def test_jdos_nac_NaCl_300K_Py(nacl_pbe: Phono3py):
"""Test running JDOS of NaCl in Py (JDOS) mode."""
jdos = _get_jdos(
nacl_pbe,
[9, 9, 9],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.set_grid_point(gp)
jdos.set_grid_point(105)
jdos.frequency_points = nacl_freq_points_at_300K
jdos.temperature = 300
jdos.run_phonon_solver()
@ -470,16 +439,14 @@ def test_jdos_nac_NaCl_300K_Py(nacl_pbe: Phono3py, gp: int, store_dense_gp_map:
)
@pytest.mark.parametrize("gp,store_dense_gp_map", [(103, False), (105, True)])
def test_jdos_nac_NaCl_300K_PyPy(nacl_pbe: Phono3py, gp: int, store_dense_gp_map: bool):
def test_jdos_nac_NaCl_300K_PyPy(nacl_pbe: Phono3py):
"""Test running JDOS of NaCl in Py (JDOS) and Py (tetrahedron) mode."""
jdos = _get_jdos(
nacl_pbe,
[9, 9, 9],
nac_params=nacl_pbe.nac_params,
store_dense_gp_map=store_dense_gp_map,
)
jdos.set_grid_point(gp)
jdos.set_grid_point(105)
jdos.frequency_points = nacl_freq_points_at_300K
jdos.temperature = 300
jdos.run_phonon_solver()
@ -490,12 +457,12 @@ def test_jdos_nac_NaCl_300K_PyPy(nacl_pbe: Phono3py, gp: int, store_dense_gp_map
)
def _get_jdos(ph3: Phono3py, mesh, nac_params=None, store_dense_gp_map=False):
def _get_jdos(ph3: Phono3py, mesh, nac_params=None):
bz_grid = BZGrid(
mesh,
lattice=ph3.primitive.cell,
symmetry_dataset=ph3.primitive_symmetry.dataset,
store_dense_gp_map=store_dense_gp_map,
store_dense_gp_map=True,
)
jdos = JointDos(
ph3.primitive,
@ -503,7 +470,6 @@ def _get_jdos(ph3: Phono3py, mesh, nac_params=None, store_dense_gp_map=False):
bz_grid,
ph3.fc2,
nac_params=nac_params,
store_dense_gp_map=store_dense_gp_map,
cutoff_frequency=1e-4,
)
return jdos

693
test/phonon3/test_real_self_energy.py
View File

@ -4,252 +4,51 @@ 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 = [
[
-0.00576660,
-0.00594616,
-0.00840087,
-0.00960344,
-0.00576660,
-0.00594616,
-0.00840087,
-0.00960344,
-0.01493508,
-0.01639729,
-0.01997820,
-0.02070427,
-0.15511645,
-0.14747203,
-0.14809966,
-0.14230763,
-0.15674925,
-0.15684992,
-0.15983868,
-0.15091767,
-0.15674925,
-0.15684992,
-0.15983868,
-0.15091767,
],
[
-0.01990306,
-0.02077094,
-0.01798066,
-0.01935581,
-0.02158076,
-0.02190634,
-0.02195633,
-0.01882258,
-0.05740055,
-0.05240406,
-0.06252644,
-0.05651015,
-0.13072273,
-0.11929265,
-0.13472599,
-0.13105120,
-0.15191900,
-0.14202698,
-0.14371246,
-0.14168892,
-0.14760248,
-0.13907618,
-0.14275290,
-0.14100562,
],
]
# 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],
[0.9207667, 0.0238034, 0.9207667, -0.1556618, 1.0742278, -0.1484696],
[1.2276889, 0.0131519, 1.2276889, -0.1515719, 1.3811500, -0.1525501],
[1.5346111, 0.0124129, 1.5346111, -0.1541357, 1.6880722, -0.1566124],
[1.8415333, 0.0165506, 1.8415333, -0.1556995, 1.9949944, -0.1567434],
[2.1484555, 0.0196960, 2.1484555, -0.1556492, 2.3019166, -0.1563791],
[2.4553778, 0.0225276, 2.4553778, -0.1553117, 2.6088389, -0.1562531],
[2.7623000, 0.0269193, 2.7623000, -0.1536529, 2.9157611, -0.1520174],
[3.0692222, 0.0247502, 3.0692222, -0.1530976, 3.2226833, -0.1545476],
[3.3761444, 0.0253261, 3.3761444, -0.1565178, 3.5296055, -0.1617560],
[3.6830666, 0.0360216, 3.6830666, -0.1586552, 3.8365277, -0.1618262],
[3.9899888, 0.0510642, 3.9899888, -0.1523272, 4.1434499, -0.1463620],
[4.2969111, 0.0494221, 4.2969111, -0.1428688, 4.4503722, -0.1388442],
[4.6038333, 0.0474328, 4.6038333, -0.1361604, 4.7572944, -0.1307086],
[4.9107555, 0.0370451, 4.9107555, -0.1345154, 5.0642166, -0.1359581],
[5.2176777, 0.0375397, 5.2176777, -0.1361887, 5.3711388, -0.1356714],
[5.5245999, 0.0345936, 5.5245999, -0.1377925, 5.6780610, -0.1390885],
[5.8315222, 0.0327132, 5.8315222, -0.1425343, 5.9849833, -0.1472336],
[6.1384444, 0.0374592, 6.1384444, -0.1489307, 6.2919055, -0.1562694],
[6.4453666, 0.0555985, 6.4453666, -0.1481088, 6.5988277, -0.1452657],
[6.7522888, 0.0586765, 6.7522888, -0.1409779, 6.9057499, -0.1375563],
[7.0592110, 0.0560470, 7.0592110, -0.1371741, 7.2126721, -0.1365273],
[7.3661333, 0.0550400, 7.3661333, -0.1372613, 7.5195944, -0.1387434],
[7.6730555, 0.0565263, 7.6730555, -0.1399341, 7.8265166, -0.1440837],
[7.9799777, 0.0624044, 7.9799777, -0.1457614, 8.1334388, -0.1582010],
[8.2868999, 0.1009214, 8.2868999, -0.1369623, 8.4403610, -0.1261905],
[8.5938221, 0.1090238, 8.5938221, -0.1087372, 8.7472832, -0.0870371],
[8.9007444, 0.0807253, 8.9007444, -0.0903838, 9.0542055, -0.0835210],
[9.2076666, 0.0658383, 9.2076666, -0.0875470, 9.3611277, -0.0853012],
[9.5145888, 0.0592222, 9.5145888, -0.0867791, 9.6680499, -0.0838069],
[9.8215110, 0.0521071, 9.8215110, -0.0853895, 9.9749721, -0.0820186],
[10.1284332, 0.0419588, 10.1284332, -0.0862870, 10.2818943, -0.0860513],
[10.4353554, 0.0364097, 10.4353554, -0.0892629, 10.5888165, -0.0892228],
[10.7422777, 0.0316193, 10.7422777, -0.0922975, 10.8957388, -0.0926032],
[11.0491999, 0.0274798, 11.0491999, -0.0955924, 11.2026610, -0.0960711],
[11.3561221, 0.0233126, 11.3561221, -0.0994131, 11.5095832, -0.1004271],
[11.6630443, 0.0193464, 11.6630443, -0.1043388, 11.8165054, -0.1065527],
[11.9699665, 0.0174292, 11.9699665, -0.1099907, 12.1234276, -0.1125414],
[12.2768888, 0.0164239, 12.2768888, -0.1157951, 12.4303499, -0.1188108],
[12.5838110, 0.0165308, 12.5838110, -0.1218120, 12.7372721, -0.1253937],
[12.8907332, 0.0184723, 12.8907332, -0.1275275, 13.0441943, -0.1309579],
[13.1976554, 0.0211215, 13.1976554, -0.1325454, 13.3511165, -0.1358727],
[13.5045776, 0.0244407, 13.5045776, -0.1369047, 13.6580387, -0.1400194],
[13.8114999, 0.0282908, 13.8114999, -0.1404205, 13.9649610, -0.1429165],
[14.1184221, 0.0313175, 14.1184221, -0.1435582, 14.2718832, -0.1462968],
[14.4253443, 0.0348581, 14.4253443, -0.1467576, 14.5788054, -0.1495649],
[14.7322665, 0.0386150, 14.7322665, -0.1500845, 14.8857276, -0.1533339],
[15.0391887, 0.0424503, 15.0391887, -0.1550110, 15.1926498, -0.1625259],
[15.3461109, 0.0605747, 15.3461109, -0.1550584, 15.4995720, -0.1547951],
[15.6530332, 0.0727408, 15.6530332, -0.1443151, 15.8064943, -0.1346934],
[15.9599554, 0.0625829, 15.9599554, -0.1346277, 16.1134165, -0.1305515],
[16.2668776, 0.0547742, 16.2668776, -0.1324223, 16.4203387, -0.1301107],
[16.5737998, 0.0459500, 16.5737998, -0.1348752, 16.7272609, -0.1365713],
[16.8807220, 0.0432888, 16.8807220, -0.1397155, 17.0341831, -0.1410867],
[17.1876443, 0.0398939, 17.1876443, -0.1450244, 17.3411054, -0.1471968],
[17.4945665, 0.0351088, 17.4945665, -0.1540409, 17.6480276, -0.1617873],
[17.8014887, 0.0424833, 17.8014887, -0.1635035, 17.9549498, -0.1699703],
[18.1084109, 0.0556675, 18.1084109, -0.1649880, 18.2618720, -0.1626782],
[18.4153331, 0.0532241, 18.4153331, -0.1631579, 18.5687942, -0.1618524],
[18.7222554, 0.0458460, 18.7222554, -0.1672846, 18.8757165, -0.1709147],
[19.0291776, 0.0443150, 19.0291776, -0.1758322, 19.1826387, -0.1808661],
[19.3360998, 0.0454289, 19.3360998, -0.1853104, 19.4895609, -0.1911626],
[19.6430220, 0.0483353, 19.6430220, -0.1952856, 19.7964831, -0.2019949],
[19.9499442, 0.0534799, 19.9499442, -0.2053378, 20.1034053, -0.2123304],
[20.2568665, 0.0599861, 20.2568665, -0.2151359, 20.4103276, -0.2225551],
[20.5637887, 0.0685856, 20.5637887, -0.2240731, 20.7172498, -0.2305568],
[20.8707109, 0.0755373, 20.8707109, -0.2331629, 21.0241720, -0.2419628],
[21.1776331, 0.0906438, 21.1776331, -0.2384003, 21.3310942, -0.2380604],
[21.4845553, 0.0827018, 21.4845553, -0.2476712, 21.6380164, -0.2593767],
[21.7914775, 0.0907297, 21.7914775, -0.2660044, 21.9449386, -0.2793854],
[22.0983998, 0.1025946, 22.0983998, -0.2859215, 22.2518609, -0.3026711],
[22.4053220, 0.1210133, 22.4053220, -0.3084041, 22.5587831, -0.3297156],
[22.7122442, 0.1517773, 22.7122442, -0.3307280, 22.8657053, -0.3535969],
[23.0191664, 0.1918115, 23.0191664, -0.3505455, 23.1726275, -0.3769971],
[23.3260886, 0.2484881, 23.3260886, -0.3659809, 23.4795497, -0.3980745],
[23.6330109, 0.3449959, 23.6330109, -0.3587633, 23.7864720, -0.3785512],
[23.9399331, 0.4932745, 23.9399331, -0.2683146, 24.0933942, -0.1740930],
[24.2468553, 0.4000672, 24.2468553, -0.1714061, 24.4003164, -0.1388929],
[24.5537775, 0.3476085, 24.5537775, -0.1478397, 24.7072386, -0.1341295],
[24.8606997, 0.3031722, 24.8606997, -0.1588841, 25.0141608, -0.1796324],
[25.1676220, 0.3332258, 25.1676220, -0.1716723, 25.3210831, -0.1741573],
[25.4745442, 0.3514490, 25.4745442, -0.1621265, 25.6280053, -0.1553650],
[25.7814664, 0.3514801, 25.7814664, -0.1503569, 25.9349275, -0.1458311],
[26.0883886, 0.3490583, 26.0883886, -0.1437619, 26.2418497, -0.1387991],
[26.3953108, 0.3346525, 26.3953108, -0.1471754, 26.5487719, -0.1497006],
[26.7022331, 0.3145614, 26.7022331, -0.1719482, 26.8556942, -0.1898755],
[27.0091553, 0.2913173, 27.0091553, -0.2370443, 27.1626164, -0.2939869],
[27.3160775, 0.2707732, 27.3160775, -0.4139632, 27.4695386, -0.6910918],
[27.6229997, 0.9436577, 27.6229997, -0.3688420, 27.7764608, -0.2026185],
[27.9299219, 0.9147269, 27.9299219, -0.1095194, 28.0833830, -0.0365244],
[28.2368441, 0.9515446, 28.2368441, 0.0642006, 28.3903052, 0.1596622],
[28.5437664, 0.9296466, 28.5437664, 0.2357042, 28.6972275, 0.3228099],
[28.8506886, 0.9383991, 28.8506886, 0.4390380, 29.0041497, 0.6142345],
[29.1576108, 0.7465231, 29.1576108, 0.6502426, 29.3110719, 0.8357126],
[29.4645330, 0.3472854, 29.4645330, 0.6927159, 29.6179941, 0.6873352],
[29.7714552, 0.2258269, 29.7714552, 0.6147208, 29.9249163, 0.6117702],
[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):
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(
@ -260,7 +59,8 @@ def test_real_self_energy_with_band_indices(si_pbesol):
write_hdf5=False,
frequency_points_at_bands=True,
)
np.testing.assert_allclose(si_pbesol_Delta, delta[0, 0, :], atol=0.01)
# 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):
@ -269,6 +69,117 @@ def test_real_self_energy_with_frequency_points(si_pbesol: Phono3py):
* 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,
-0.00840087,
-0.00960344,
-0.00576660,
-0.00594616,
-0.00840087,
-0.00960344,
-0.01493508,
-0.01639729,
-0.01997820,
-0.02070427,
-0.15511645,
-0.14747203,
-0.14809966,
-0.14230763,
-0.15674925,
-0.15684992,
-0.15983868,
-0.15091767,
-0.15674925,
-0.15684992,
-0.15983868,
-0.15091767,
],
[
-0.01990306,
-0.02077094,
-0.01798066,
-0.01935581,
-0.02158076,
-0.02190634,
-0.02195633,
-0.01882258,
-0.05740055,
-0.05240406,
-0.06252644,
-0.05651015,
-0.13072273,
-0.11929265,
-0.13472599,
-0.13105120,
-0.15191900,
-0.14202698,
-0.14371246,
-0.14168892,
-0.14760248,
-0.13907618,
-0.14275290,
-0.14100562,
],
]
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.init_phph_interaction()
frequency_points = [1.469947, 3.085309, 14.997187, 15.129080]
@ -282,18 +193,14 @@ def test_real_self_energy_with_frequency_points(si_pbesol: Phono3py):
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))
# 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):
@ -303,22 +210,290 @@ def test_real_self_energy_nacl_nac_npoints(nacl_pbe: Phono3py):
* 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 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
)
# 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],
[0.9207667, 0.0238034, 0.9207667, -0.1556618, 1.0742278, -0.1484696],
[1.2276889, 0.0131519, 1.2276889, -0.1515719, 1.3811500, -0.1525501],
[1.5346111, 0.0124129, 1.5346111, -0.1541357, 1.6880722, -0.1566124],
[1.8415333, 0.0165506, 1.8415333, -0.1556995, 1.9949944, -0.1567434],
[2.1484555, 0.0196960, 2.1484555, -0.1556492, 2.3019166, -0.1563791],
[2.4553778, 0.0225276, 2.4553778, -0.1553117, 2.6088389, -0.1562531],
[2.7623000, 0.0269193, 2.7623000, -0.1536529, 2.9157611, -0.1520174],
[3.0692222, 0.0247502, 3.0692222, -0.1530976, 3.2226833, -0.1545476],
[3.3761444, 0.0253261, 3.3761444, -0.1565178, 3.5296055, -0.1617560],
[3.6830666, 0.0360216, 3.6830666, -0.1586552, 3.8365277, -0.1618262],
[3.9899888, 0.0510642, 3.9899888, -0.1523272, 4.1434499, -0.1463620],
[4.2969111, 0.0494221, 4.2969111, -0.1428688, 4.4503722, -0.1388442],
[4.6038333, 0.0474328, 4.6038333, -0.1361604, 4.7572944, -0.1307086],
[4.9107555, 0.0370451, 4.9107555, -0.1345154, 5.0642166, -0.1359581],
[5.2176777, 0.0375397, 5.2176777, -0.1361887, 5.3711388, -0.1356714],
[5.5245999, 0.0345936, 5.5245999, -0.1377925, 5.6780610, -0.1390885],
[5.8315222, 0.0327132, 5.8315222, -0.1425343, 5.9849833, -0.1472336],
[6.1384444, 0.0374592, 6.1384444, -0.1489307, 6.2919055, -0.1562694],
[6.4453666, 0.0555985, 6.4453666, -0.1481088, 6.5988277, -0.1452657],
[6.7522888, 0.0586765, 6.7522888, -0.1409779, 6.9057499, -0.1375563],
[7.0592110, 0.0560470, 7.0592110, -0.1371741, 7.2126721, -0.1365273],
[7.3661333, 0.0550400, 7.3661333, -0.1372613, 7.5195944, -0.1387434],
[7.6730555, 0.0565263, 7.6730555, -0.1399341, 7.8265166, -0.1440837],
[7.9799777, 0.0624044, 7.9799777, -0.1457614, 8.1334388, -0.1582010],
[8.2868999, 0.1009214, 8.2868999, -0.1369623, 8.4403610, -0.1261905],
[8.5938221, 0.1090238, 8.5938221, -0.1087372, 8.7472832, -0.0870371],
[8.9007444, 0.0807253, 8.9007444, -0.0903838, 9.0542055, -0.0835210],
[9.2076666, 0.0658383, 9.2076666, -0.0875470, 9.3611277, -0.0853012],
[9.5145888, 0.0592222, 9.5145888, -0.0867791, 9.6680499, -0.0838069],
[9.8215110, 0.0521071, 9.8215110, -0.0853895, 9.9749721, -0.0820186],
[10.1284332, 0.0419588, 10.1284332, -0.0862870, 10.2818943, -0.0860513],
[10.4353554, 0.0364097, 10.4353554, -0.0892629, 10.5888165, -0.0892228],
[10.7422777, 0.0316193, 10.7422777, -0.0922975, 10.8957388, -0.0926032],
[11.0491999, 0.0274798, 11.0491999, -0.0955924, 11.2026610, -0.0960711],
[11.3561221, 0.0233126, 11.3561221, -0.0994131, 11.5095832, -0.1004271],
[11.6630443, 0.0193464, 11.6630443, -0.1043388, 11.8165054, -0.1065527],
[11.9699665, 0.0174292, 11.9699665, -0.1099907, 12.1234276, -0.1125414],
[12.2768888, 0.0164239, 12.2768888, -0.1157951, 12.4303499, -0.1188108],
[12.5838110, 0.0165308, 12.5838110, -0.1218120, 12.7372721, -0.1253937],
[12.8907332, 0.0184723, 12.8907332, -0.1275275, 13.0441943, -0.1309579],
[13.1976554, 0.0211215, 13.1976554, -0.1325454, 13.3511165, -0.1358727],
[13.5045776, 0.0244407, 13.5045776, -0.1369047, 13.6580387, -0.1400194],
[13.8114999, 0.0282908, 13.8114999, -0.1404205, 13.9649610, -0.1429165],
[14.1184221, 0.0313175, 14.1184221, -0.1435582, 14.2718832, -0.1462968],
[14.4253443, 0.0348581, 14.4253443, -0.1467576, 14.5788054, -0.1495649],
[14.7322665, 0.0386150, 14.7322665, -0.1500845, 14.8857276, -0.1533339],
[15.0391887, 0.0424503, 15.0391887, -0.1550110, 15.1926498, -0.1625259],
[15.3461109, 0.0605747, 15.3461109, -0.1550584, 15.4995720, -0.1547951],
[15.6530332, 0.0727408, 15.6530332, -0.1443151, 15.8064943, -0.1346934],
[15.9599554, 0.0625829, 15.9599554, -0.1346277, 16.1134165, -0.1305515],
[16.2668776, 0.0547742, 16.2668776, -0.1324223, 16.4203387, -0.1301107],
[16.5737998, 0.0459500, 16.5737998, -0.1348752, 16.7272609, -0.1365713],
[16.8807220, 0.0432888, 16.8807220, -0.1397155, 17.0341831, -0.1410867],
[17.1876443, 0.0398939, 17.1876443, -0.1450244, 17.3411054, -0.1471968],
[17.4945665, 0.0351088, 17.4945665, -0.1540409, 17.6480276, -0.1617873],
[17.8014887, 0.0424833, 17.8014887, -0.1635035, 17.9549498, -0.1699703],
[18.1084109, 0.0556675, 18.1084109, -0.1649880, 18.2618720, -0.1626782],
[18.4153331, 0.0532241, 18.4153331, -0.1631579, 18.5687942, -0.1618524],
[18.7222554, 0.0458460, 18.7222554, -0.1672846, 18.8757165, -0.1709147],
[19.0291776, 0.0443150, 19.0291776, -0.1758322, 19.1826387, -0.1808661],
[19.3360998, 0.0454289, 19.3360998, -0.1853104, 19.4895609, -0.1911626],
[19.6430220, 0.0483353, 19.6430220, -0.1952856, 19.7964831, -0.2019949],
[19.9499442, 0.0534799, 19.9499442, -0.2053378, 20.1034053, -0.2123304],
[20.2568665, 0.0599861, 20.2568665, -0.2151359, 20.4103276, -0.2225551],
[20.5637887, 0.0685856, 20.5637887, -0.2240731, 20.7172498, -0.2305568],
[20.8707109, 0.0755373, 20.8707109, -0.2331629, 21.0241720, -0.2419628],
[21.1776331, 0.0906438, 21.1776331, -0.2384003, 21.3310942, -0.2380604],
[21.4845553, 0.0827018, 21.4845553, -0.2476712, 21.6380164, -0.2593767],
[21.7914775, 0.0907297, 21.7914775, -0.2660044, 21.9449386, -0.2793854],
[22.0983998, 0.1025946, 22.0983998, -0.2859215, 22.2518609, -0.3026711],
[22.4053220, 0.1210133, 22.4053220, -0.3084041, 22.5587831, -0.3297156],
[22.7122442, 0.1517773, 22.7122442, -0.3307280, 22.8657053, -0.3535969],
[23.0191664, 0.1918115, 23.0191664, -0.3505455, 23.1726275, -0.3769971],
[23.3260886, 0.2484881, 23.3260886, -0.3659809, 23.4795497, -0.3980745],
[23.6330109, 0.3449959, 23.6330109, -0.3587633, 23.7864720, -0.3785512],
[23.9399331, 0.4932745, 23.9399331, -0.2683146, 24.0933942, -0.1740930],
[24.2468553, 0.4000672, 24.2468553, -0.1714061, 24.4003164, -0.1388929],
[24.5537775, 0.3476085, 24.5537775, -0.1478397, 24.7072386, -0.1341295],
[24.8606997, 0.3031722, 24.8606997, -0.1588841, 25.0141608, -0.1796324],
[25.1676220, 0.3332258, 25.1676220, -0.1716723, 25.3210831, -0.1741573],
[25.4745442, 0.3514490, 25.4745442, -0.1621265, 25.6280053, -0.1553650],
[25.7814664, 0.3514801, 25.7814664, -0.1503569, 25.9349275, -0.1458311],
[26.0883886, 0.3490583, 26.0883886, -0.1437619, 26.2418497, -0.1387991],
[26.3953108, 0.3346525, 26.3953108, -0.1471754, 26.5487719, -0.1497006],
[26.7022331, 0.3145614, 26.7022331, -0.1719482, 26.8556942, -0.1898755],
[27.0091553, 0.2913173, 27.0091553, -0.2370443, 27.1626164, -0.2939869],
[27.3160775, 0.2707732, 27.3160775, -0.4139632, 27.4695386, -0.6910918],
[27.6229997, 0.9436577, 27.6229997, -0.3688420, 27.7764608, -0.2026185],
[27.9299219, 0.9147269, 27.9299219, -0.1095194, 28.0833830, -0.0365244],
[28.2368441, 0.9515446, 28.2368441, 0.0642006, 28.3903052, 0.1596622],
[28.5437664, 0.9296466, 28.5437664, 0.2357042, 28.6972275, 0.3228099],
[28.8506886, 0.9383991, 28.8506886, 0.4390380, 29.0041497, 0.6142345],
[29.1576108, 0.7465231, 29.1576108, 0.6502426, 29.3110719, 0.8357126],
[29.4645330, 0.3472854, 29.4645330, 0.6927159, 29.6179941, 0.6873352],
[29.7714552, 0.2258269, 29.7714552, 0.6147208, 29.9249163, 0.6117702],
[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],
]
vals = np.array(im_part)
i2r = ImagToReal(vals[:, 1], vals[:, 0])
i2r.run()

927
test/phonon3/test_spectral_function.py
View File

@ -1,259 +1,620 @@
"""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."""
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,
],
[
-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,
],
],
]
si_pbesol.mesh_numbers = [9, 9, 9]
if check_bands:
si_pbesol.band_indices = [[4, 5]]
si_pbesol.init_phph_interaction()
sf = SpectralFunction(
si_pbesol.phph_interaction,
@ -266,49 +627,35 @@ def test_SpectralFunction(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.array(ref_shifts)[:, [4, 5], :],
sf.shifts[0, 0],
atol=1e-2,
)
np.testing.assert_allclose(
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]), "],")
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,
)
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
)
def test_SpectralFunction_band_indices(si_pbesol: Phono3py):
"""Spectral function of Si."""
si_pbesol.mesh_numbers = [9, 9, 9]
si_pbesol.band_indices = [[4, 5]]
si_pbesol.init_phph_interaction()
sf = SpectralFunction(
si_pbesol.phph_interaction,
si_pbesol.grid.grg2bzg[[1, 103]],
temperatures=[
300,
],
num_frequency_points=10,
log_level=1,
)
sf.run()
# for line in np.swapaxes(sf.spectral_functions, -2, -1).reshape(-1, 6):
# print(("%.7f, " * 6) % tuple(line))
# raise
np.testing.assert_allclose(
np.reshape(shifts, (-1, 6))[:, [4, 5]],
np.swapaxes(sf.shifts, -2, -1).reshape(-1, 2),
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),
atol=1e-2,
rtol=1e-2,
)
if check_bands:
si_pbesol.band_indices = None

3
test/phonon3/test_triplets.py
View File

@ -5,6 +5,7 @@ from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.symmetry import Symmetry
from phono3py import Phono3py
from phono3py.phonon3.triplets import (
_get_BZ_triplets_at_q,
_get_triplets_reciprocal_mesh_at_q,
@ -68,7 +69,7 @@ def test_get_triplets_at_q_type1(si_pbesol_111):
np.testing.assert_equal(weights, weights_ref)
def test_get_triplets_at_q_type2(si_pbesol_111):
def test_get_triplets_at_q_type2(si_pbesol_111: Phono3py):
"""Test triplets under type2 grid."""
pcell = si_pbesol_111.primitive
psym = si_pbesol_111.primitive_symmetry